# Learning Nonlinear State Space Models with Hamiltonian Sequential Monte
Carlo Sampler

Duo Xu
Department of Electrical and Computer Engineering
Georgian Institute of Technology
Atlanta, GA 30332

###### Abstract

State space models (SSM) have been widely applied for the analysis and
visualization of large sequential datasets. Sequential Monte Carlo (SMC) is a
very popular particle-based method to sample latent states from intractable
posteriors. However, SSM is significantly influenced by the choice of the
proposal. Recently Hamiltonian Monte Carlo (HMC) sampling has shown success in
many practical problems. In this paper, we propose an SMC augmented by HMC
(HSMC) for inference and model learning of nonlinear SSM, which can exempt us
from learning proposals and reduce the model complexity significantly. Based
on the measure preserving property of HMC, the particles directly generated by
transition function can approximate the posterior of latent states arbitrarily
well. In order to better adapt to the local geometry of latent space, the HMC
is conducted on Riemannian manifold defined by a positive definite metric
$\bm{M}(\bm{x})$. In addition, we show that the proposed HSMC method can
improve SSMs realized by both Gaussian Processes (GP) and Neural Network (NN).

##  1 Introduction

System identification (Ljung,, 1998, 2010) is a fundamental ingredient of many
problems, such as model-predictive control (Camacho and Alba,, 2013) and
model-based reinforcement learning (Deisenroth and Rasmussen,, 2011;
Berkenkamp et al.,, 2017), which is to learn system dynamics from practical
data. State-space model (SSM) (Billings,, 2013) is the most popular method to
represent the system with input $\bm{u}_{\\_}t$ and output $\bm{y}_{\\_}t$ as
functions of a latent Markovian state $\bm{x}_{\\_}t$. Specifically, linear
and non-linear Gaussian state space models (GSSM) are most widely used in
practical applications from robotic planning to neural signal processing.
However, despite significant effort in research community over past decades,
efficient learning method for non-linear GSSM is still lacking.

The sequential Monte Carlo (SMC) (Gordon et al.,, 1993) is a classical method
to infer latent state in SSM, which uses weighted particles to approximate the
intractable posterior of the latent states. The proposal distribution, from
which particles are sampled, has significant influence on the approximation
performance. In order to perform model learning and proposal adaptation at the
same time, recent work (Maddison et al.,, 2017; Le et al.,, 2017; Naesseth et
al.,, 2017) combines variational auto-encoder (VAE) (Kingma and Welling,,
2013; Rezende et al.,, 2014) with importance weighted auto-encoder (IWAE)
(Burda et al.,, 2015), and uses SMC as the estimator for marginal observation
likelihood. However, since in nonlinear GSSM the emission and transition
frameworks are realized by deep neural networks, the real posterior of latent
states is intractable to sample, and the proposal distributions in previous
work are always assumed to be Gaussian, different from the true posterior.
Although the approximated log-likelihood is unbiased asymptotically (Maddison
et al.,, 2017), it is not close to the real log-likelihood without large
number of particles.

In order to improve the inference performance of SMC with low sampling and
model complexity, here we propose an SMC sampler augmented by Hamiltonian
dynamics (HSMC). Different from previous SMC methods, we don’t need proposals
here to generate particles of latent states. Using transition function as
initial distribution, we use HMC to sample particles to approximate the
posterior. And since latent space is time-dependent in nonlinear SSM, we
modified the Riemannian Manifold HMC (RMHMC) (Girolami and Calderhead,, 2011),
whose mass matrix is generated by an MLP with latent states as inputs. The
state dynamics in SSM is usually realized by Gaussian Processes (GP) and
neural network. GP is data-efficient and has much less parameters than neural
network, but its expressive capability is not as good as neural network,
especially in high-dimensional problems. Here we show that the proposed method
can improve the learning performance of GP to be comparable as neural network.

##  2 Preliminary

Denote $\bm{u}_{\\_}t,\bm{x}_{\\_}t$ and $\bm{y}_{\\_}t$ as input, latent and
output variables at time $t$ respectively.

###  2.1 Gaussian State Space Model

Gaussian state space model (GSSM) is the most popular method to model the
dynamics of complex sequential data in the latent space (Raiko and Tornio,,
2009). The inference and learning of GSSM are both considered in this work.
The model is defined as

$\displaystyle\bm{x}_{\\_}t$ $\displaystyle\sim$
$\displaystyle\mathcal{N}(\mu_{\\_}{\theta}(\bm{x}_{\\_}{t-1},\bm{u}_{\\_}t),\sigma^{2}_{\\_}{\theta}(\bm{x}_{\\_}{t-1},\bm{u}_{\\_}t))\hskip
25.0pt\text{Transition}$ $\displaystyle\bm{y}_{\\_}t$ $\displaystyle\sim$
$\displaystyle\Pi(f_{\\_}{\theta}(\bm{x}_{\\_}{t-1},\bm{u}_{\\_}t))\hskip
90.0pt\text{Emission}$ (1)

where $\Pi$ is the output distribution parametrized by function $f$, and
latent variables distribute as multivariate Gaussian conditioned on previous
latent variables and input variables. The GSSM defined above includes both
linear and nonlinear GSSM. When functions $f$ and $\mu$ are linear, the model
can be learned by extended Kalman filter (Wan and Nelson,, 1997) and
expectation maximization (Ghahramani and Roweis,, 1999). However, in most
practical problems, the dynamic and emission functions are nonlinear. In this
paper we propose an efficient method to deal with nonlinearties.

###  2.2 Variational Sequential Monte Carlo

SMC performs particle-based approximate inference on a sequence of target
distributions, which is the extension of importance sampling to sequential
data. In the context of SSM, the target distributions are the posterior of
latent variables, i.e.,
$\\{p_{\\_}{\theta}(\bm{x}_{\\_}{1:t}|\bm{y}_{\\_}{1:t})\\}_{\\_}{t=1}^{T}$.
The generative model includes initial distribution of latent variables
$p_{\\_}{\theta}(\bm{x}_{\\_}0)$, transition distribution
$p_{\\_}{\theta}(\bm{x}_{\\_}{t}|\bm{x}_{\\_}{t-1},\bm{u}_{\\_}t)$, and
emission distribution $p_{\\_}{\theta}(\bm{y}_{\\_}t|\bm{x}_{\\_}t)$. The
proposal distribution, defined as
$q_{\\_}{\phi}(\bm{x}_{\\_}t|\bm{x}_{\\_}{t-1},\bm{y}_{\\_}t,\bm{u}_{\\_}t)$,
is the approximate inference on target distributions. Then the generative and
inference models are factorized as

$\displaystyle p_{\\_}{\theta}(\bm{x}_{\\_}{\leq T},\bm{y}_{\\_}{\leq
T}|\bm{u}_{\\_}{\leq T})$ $\displaystyle=$ $\displaystyle
p_{\\_}{\theta}(\bm{x}_{\\_}0)\prod_{\\_}{t=1}^{T}p_{\\_}{\theta}(\bm{x}_{\\_}t|\bm{x}_{\\_}{t-1},\bm{u}_{\\_}t)p_{\\_}{\theta}(\bm{y}_{\\_}t|\bm{x}_{\\_}t)$
$\displaystyle q_{\\_}{\theta}(\bm{x}_{\\_}{\leq T}|\bm{y}_{\\_}{\leq
T},\bm{u}_{\\_}{\leq T})$ $\displaystyle=$ $\displaystyle
p_{\\_}{\theta}(\bm{x}_{\\_}0)\prod_{\\_}{t=1}^{T}q_{\\_}{\theta}(\bm{x}_{\\_}t|\bm{x}_{\\_}{t-1},\bm{y}_{\\_}t,\bm{u}_{\\_}t)$

Recently, a new ELBO objective has been introduced (Le et al.,, 2017; Maddison
et al.,, 2017; Naesseth et al.,, 2017), which is asymptotically unbiased
estimator on log-likelihood

$\text{ELBO}_{\\_}{\text{SMC}}=\mathbb{E}\bigg{[}\sum_{\\_}{t=1}^{T}\log\bigg{(}\frac{1}{K}\sum_{\\_}{k=1}^{K}\omega_{\\_}t^{k}\bigg{)}\bigg{]}$
(2)

where $K$ is the number of particles and $\omega_{\\_}t^{k}$ is the weight of
particle $k$ at time $t$. Each particle is defined by weight
$\omega_{\\_}t^{k}$ and value $\bm{x}_{\\_}t^{p}$. At time $t=0$ each particle
value $\bm{x}_{\\_}t^{p}$ is sampled from initial latent distribution
$p_{\\_}{\theta}(\bm{x}_{\\_}0)$. In this paper, the transition (emission)
function is denoted as $f(\cdot|\cdot)$ ( $g(\cdot|\cdot)$ ). By resampling
from the previous particle set $\\{\bm{x}_{\\_}{t-1}^{k}\\}_{\\_}{k=1}^{K}$,
the weight $\omega_{\\_}t^{k}$ for every particle at each time is defined as
below

$\omega_{\\_}t^{k}=\frac{f_{\\_}{\theta}(\bm{x}_{\\_}t^{k}|\bm{x}_{\\_}{t-1}^{\alpha_{\\_}t^{k}},\bm{u}_{\\_}t,\bm{y}_{\\_}{t-1})g_{\\_}{\theta}(\bm{y}_{\\_}t|\bm{x}_{\\_}t^{k})}{q_{\\_}{\phi}(\bm{x}_{\\_}t^{k}|\bm{x}_{\\_}{t-1}^{\alpha_{\\_}t^{k}},\bm{y}_{\\_}t,\bm{u}_{\\_}t)}$
(3)

where $\bm{x}_{\\_}t^{k}$ is sampled from the proposal
$q_{\\_}{\phi}(\cdot|\bm{x}_{\\_}{t-1}^{\alpha_{\\_}t^{k}},\bm{y}_{\\_}t,\bm{u}_{\\_}t)$,
and the index follows

$\alpha_{\\_}t^{k}\sim\text{Discrete}\bigg{(}\bigg{\\{}\frac{\omega_{\\_}{t-1}^{k}}{\sum_{\\_}{l=1}^{K}\omega_{\\_}{t-1}^{l}}\bigg{\\}}_{\\_}{k=1}^{K}\bigg{)}$

with ancestor index $\alpha_{\\_}1^{k}=k$. Following (Le et al.,, 2017) we
know that

$\mathbb{E}\bigg{[}\prod_{\\_}{t=1}^{T}\frac{1}{K}\sum_{\\_}{k=1}^{K}\omega_{\\_}t^{k}\bigg{]}:=\mathbb{E}\bigg{[}\hat{Z}_{\\_}{\text{SMC}}(\bm{x}_{\\_}{1:T}^{1:K},\alpha_{\\_}{1:T-1}^{1:K})\bigg{]}=p_{\\_}{\theta}(\bm{y}_{\\_}{1:T})$
(4)

showing that product of particle weights is an unbiased estimator to the
marginal likelihood of observations.

###  2.3 Hamiltonian Monte Carlo

In learning transition and emission of SSM, it is key to generate samples from
posterior distribution of latent states given observation data. However, the
real posterior is intractable in nonlinear SSM. In this work, different
previous methods (Le et al.,, 2017; Maddison et al.,, 2017; Naesseth et al.,,
2017), in each time step we don’t learn proposals but use HMC (Neal et al.,,
2011) to directly sample latent states from the posterior, which can reduce
the model complexity significantly. At each time step $t$, we denote the joint
log likelihood of the observations and latent states as

$\mathcal{L}(\bm{x}_{\\_}t)=\log g(\bm{y}_{\\_}t|\bm{x}_{\\_}t)+\log
f(\bm{x}_{\\_}t|\bm{x}_{\\_}{t-1},\bm{y}_{\\_}{t-1},\bm{u}_{\\_}t)$ (5)

summing log probability of emission and transition function together. Then we
add to it a term involving "momenta" variables $\bm{p}_{\\_}t$, to obtain the
Hamiltonian energy function,

$H(\bm{x}_{\\_}t,\bm{p}_{\\_}t)=-\mathcal{L}(\bm{x}_{\\_}t)+\frac{1}{2}\bm{p}_{\\_}t^{T}M^{-1}\bm{p}_{\\_}t$
(6)

This quantity can be interpreted as in physical terminology as the sum of the
potential energy $\mathcal{L}(\bm{x}_{\\_}t)$ and the kinetic energy
$\frac{1}{2}\bm{p}_{\\_}t^{T}M^{-1}\bm{p}_{\\_}t$, where $M$ acts as the
canonical mass matrix. The joint distribution of $\bm{x}_{\\_}t$ and
$\bm{p}_{\\_}t$ is then defined as
$p(\bm{x}_{\\_}t,\bm{p}_{\\_}t)\propto\exp(-H(\bm{x}_{\\_}t,\bm{p}_{\\_}t))$.

Denote the time derivatives with the dot notation, i.e.,
$\dot{\bm{x}_{\\_}t}=d\bm{x}_{\\_}t/d\tau$, where $\tau$ is the refined time
in $[t,t+1)$. The Hamiltonian equations of motion governing the dynamics of
this system can be written as

$\dot{\bm{x}_{\\_}t}=\frac{\partial}{\partial\bm{p}_{\\_}t}H(\bm{x}_{\\_}t,\bm{p}_{\\_}t),\hskip
20.0pt\dot{\bm{p}_{\\_}t}=-\frac{\partial}{\partial\bm{x}_{\\_}t}H(\bm{x}_{\\_}t,\bm{p}_{\\_}t)$

Obviously these equations are time-reversible, and the dynamics conserve the
total energy. These continuous-time equations can be discretized to give
"leapfrog" algorithms which are used for Monte Carlo simulations along with
Metropolis-Hastings correction steps (Neal et al.,, 2011).

##  3 Model Learning Method

In this section, the proposed Hamiltonian Sequential Monte Carlo (HSMC) is
thoroughly described. We first extend Riemann manifold HMC to recurrent
setting, and then formulate the learning objective function based on
variational inference. Finally, based on HSMC, we will introduce a new meta-
learning method for nonstationary state space model.

###  3.1 Riemann Manifold HMC

The sampling quality of HMC is heavily influenced by the choice of mass matrix
(Girolami and Calderhead,, 2011). In order to automatically select mass matrix
across different time steps, we adopt Riemann manifold HMC (Girolami and
Calderhead,, 2011) to incorporate local geometric properties of latent states.
Here, at time step $t$, we parametrize the mass matrix as a function of latent
states, i.e., $\bm{M}(\bm{x}_{\\_}t)$, and it is a positive definite metric
tensor defining the Riemann manifold on which we are sampling. Defining the
kinetic energy in terms of the mass matrix, we can get the Hamiltonian energy
function as below

$H_{\\_}{\text{rm}}(\bm{x}_{\\_}t,\bm{p}_{\\_}t)=-\mathcal{L}(\bm{x}_{\\_}t)+\frac{1}{2}\bm{p}_{\\_}t^{T}\bm{M}(\bm{x}_{\\_}t)^{-1}\bm{p}_{\\_}t+\frac{1}{2}\log\\{(2\pi)^{D_{\\_}x}|\bm{M}(\bm{p}_{\\_}t)|\\}$
(7)

where $D_{\\_}x$ is the dimension of latent state. The desired marginal
density of $\bm{x}_{\\_}t$ can be obtained by integrating out the momenta
$\bm{p}_{\\_}t$. The equations of motion at every time step $t$ are as below

$\displaystyle\dot{\bm{x}_{\\_}t}$ $\displaystyle=$
$\displaystyle\bm{M}(\bm{x}_{\\_}t)^{-1}\bm{p}_{\\_}t$ (8)
$\displaystyle\dot{\bm{p}_{\\_}t}$ $\displaystyle=$
$\displaystyle\nabla\mathcal{L}(\bm{x}_{\\_}t)-\frac{1}{2}\text{tr}\big{(}\bm{M}(\bm{x}_{\\_}t)^{-1}\nabla\bm{M}(\bm{x}_{\\_}t)\big{)}$
(9)
$\displaystyle+\frac{1}{2}\bm{p}_{\\_}t^{T}\bm{M}(\bm{x}_{\\_}t)^{-1}\nabla\bm{M}(\bm{x}_{\\_}t)\bm{M}(\bm{x}_{\\_}t)^{-1}\bm{p}_{\\_}t$

where the last equation is denoted as a function $\hat{U}(\cdot)$ of latent
state. To discretize this system of equations, we use the generalized leapfrog
algorithm, where a first order symplectic integrator is composed with its
adjoint; the resultant second order integrator can be shown to be both time-
reversible and symplectic (Leimkuhler and Reich,, 2004). The Riemann manifold
HMC operation is formulated as a function shown as below.

Algorithm 1 Riemann Manifold HMC function:
RMHMC($\bm{x},\bm{p},\bm{M},\mathcal{L},S,\epsilon$;$\phi$)

0: $\bm{M}(\cdot)$: mass matrix function parameterized by $\phi$;
$\mathcal{L}(\cdot)$: joint log likelihood function

0: $S$: number of HMC steps; $\epsilon$: step size ($S\epsilon<1$)

Define
$\hat{U}(\bm{x},\bm{p})=\nabla\mathcal{L}(\bm{x})-\frac{1}{2}\text{tr}\big{(}\bm{M}(\bm{x})^{-1}\nabla\bm{M}(\bm{x})\big{)}+\frac{1}{2}\bm{p}^{T}\bm{M}(\bm{x})^{-1}\nabla\bm{M}(\bm{x})\bm{M}(\bm{x})^{-1}\bm{p}$

for $s\leftarrow 1$ to $S$ do

$\tilde{\bm{p}}\leftarrow\bm{p}-\epsilon/2\odot\hat{U}(\bm{x},\bm{p})$

$\bm{x}\leftarrow\bm{x}+\epsilon\odot(\bm{M}(\bm{x})^{-1}\tilde{\bm{p}})$

$\bm{p}\leftarrow\tilde{\bm{p}}-\epsilon/2\odot\hat{U}(\bm{x},\tilde{\bm{p}})$

end for

return augmented latent state $\bm{x}$ and momenta variables $\bm{p}$

In implementation, the mass matrix is realized by a low-rank matrix, i.e.,
$\bm{M}(\bm{x})=\bm{u}(\bm{x})+\bm{\nu}(\bm{x})\bm{\nu}(\bm{x})^{T}$, where
$\bm{u}(\bm{x}),\bm{\nu}(\bm{x})\in\mathbb{R}^{D_{\\_}x\times 1}$ are MLPs
with current latent states as inputs, and their parameters are denoted as
$\phi$. Then we denote mass matrix as $\bm{M}_{\\_}{\phi}$.

###  3.2 Hamiltonian SMC

Here we describe the proposed method, Hamiltonian Sequential Monte Carlo
(HSMC). At each time step, the particles of latent states are directly sampled
from the transition function conditioned on previous particles
$\tilde{\bm{x}}_{\\_}{t-1}^{\alpha_{\\_}{t-1}^{k}}$, current input variables
$\bm{u}_{\\_}t$ and previous observations $\bm{y}_{\\_}{t-1}$, which plays the
role of proposal in previous work (Le et al.,, 2017; Maddison et al.,, 2017;
Naesseth et al.,, 2017). The momenta variables are introduced necessarily,
sampled from Gaussian distribution $\mathcal{N}(0,\bm{M}_{\\_}{\phi})$. And
each particle consists of a tuple of latent state and momenta variable. Then
we use the Riemann Manifold Hamiltonian Monte Carlo (RMHMC) to transform the
sampled particles to follow the posterior of latent states given current
observations. For each particle $k$, denote
$(\bm{x}_{\\_}t^{0,k},\bm{p}_{\\_}t^{0,k})$ and
$(\bm{x}_{\\_}t^{S,k},\bm{p}_{\\_}t^{S,k})$ as initial particle and
transformed particle with $S$-step RMHMC. Due to the measure preserving
property of HMC (Neal et al.,, 2011), the initial and transformed particles
have the same density value, i.e.,
$q^{0}_{\\_}t(\bm{x}_{\\_}t^{0,k},\bm{p}_{\\_}t^{0,k})=q^{S}_{\\_}t(\bm{x}_{\\_}t^{S,k},\bm{p}_{\\_}t^{S,k})$,
even though they follow different distributions. Then the weight for $k$-th
particle can be defined as

$\displaystyle\omega_{\\_}t^{k}$ $\displaystyle=$
$\displaystyle\frac{g_{\\_}{\theta}(\bm{y}_{\\_}t|\bm{x}_{\\_}t^{S,k})f_{\\_}{\theta}(\bm{x}^{S,k}_{\\_}t|\tilde{\bm{x}}_{\\_}t^{\alpha^{k}_{\\_}{t-1}},\bm{u}_{\\_}t,\bm{y}_{\\_}{t-1})\mathcal{N}(\bm{p}_{\\_}t^{S,k}|0,\bm{M}_{\\_}{\phi}(\bm{x}_{\\_}t^{S,k}))}{q^{S}_{\\_}t(\bm{x}_{\\_}t^{S,k},\bm{p}_{\\_}t^{S,k})}$
(10) $\displaystyle=$
$\displaystyle\frac{g_{\\_}{\theta}(\bm{y}_{\\_}t|\bm{x}^{S,k}_{\\_}t)f_{\\_}{\theta}(\bm{x}^{S,k}_{\\_}t|\tilde{\bm{x}}_{\\_}t^{\alpha^{k}_{\\_}{t-1}},\bm{u}_{\\_}t,\bm{y}_{\\_}{t-1})\mathcal{N}(\bm{p}^{S,k}_{\\_}t|0,\bm{M}_{\\_}{\phi}(\bm{x}_{\\_}t^{S,k}))}{f_{\\_}{\theta}(\bm{x}^{0,k}_{\\_}t|\tilde{\bm{x}}_{\\_}t^{\alpha^{k}_{\\_}{t-1}},\bm{u}_{\\_}t,\bm{y}_{\\_}{t-1})\mathcal{N}(\bm{p}^{0,k}_{\\_}t|0,\bm{M}(\bm{x}_{\\_}t^{0,k}))}$

where the second equation comes from the measure preserving property of HMC.
Therefore, the proposed algorithm is summarized as below.

Algorithm 2 Hamiltonian Sequential Monte Carlo (HSMC)

0: observation data $\bm{y}_{\\_}{1:T}$, control variables $\bm{u}_{\\_}{1:T}$

0: model parameters $\theta$, mass matrix $\mathcal{M}_{\\_}{\phi}$ with
parameters $\phi$.

Sample initial particle values $\bm{x}_{\\_}1^{0,k}\sim
p_{\\_}{1,\theta}(\cdot|\bm{u}_{\\_}1)$.

Compute particle weights
$\omega_{\\_}1^{k}={g_{\\_}{\theta}(\bm{y}_{\\_}1|\bm{x}_{\\_}1^{0,k})}$.

Initialize particle set
$\tilde{\bm{x}}_{\\_}1^{k}\leftarrow\bm{x}_{\\_}1^{0,k}$

for $t=2,\ldots,T$ do

Sample ancestor index
$\alpha^{k}_{\\_}{t-1}\sim\text{Discrete}(\cdot|\omega_{\\_}{t-1}^{1},\ldots,\omega_{\\_}{t-1}^{K})$.

Sample particle value from transition function $\bm{x}_{\\_}t^{0,k}\sim
f_{\\_}{\theta}(\cdot|\bm{u}_{\\_}t,\bm{y}_{\\_}{t-1},\tilde{\bm{x}}_{\\_}{t-1}^{\alpha_{\\_}{t-1}^{k}})$.

Sample momenta variables from Gaussian distribution
$\bm{p}_{\\_}t^{k}\sim\mathcal{N}(0,\bm{M}_{\\_}{\phi}(\bm{x}_{\\_}t^{0,k}))$.

Process each particle value by RMHMC operation
$\bm{x}_{\\_}t^{S,k},\bm{p}_{\\_}t^{S,k}=\text{RMHMC}(\bm{x}_{\\_}t^{0,k},\bm{p}^{0,k}_{\\_}t,\bm{M}_{\\_}{\phi},\mathcal{L},S,\epsilon;\phi)$
(11)

Update particle set
$\tilde{\bm{x}}_{\\_}t^{k}\leftarrow\\{\tilde{\bm{x}}^{k}_{\\_}t,\bm{x}_{\\_}t^{S,k}\\}$

Compute particle weights $\omega_{\\_}t^{k}$ in (10).

end for

Compute marginal likelihood
$\hat{Z}_{\\_}{\text{HSMC}}=\prod_{\\_}{t=1}^{T}\frac{1}{K}\sum_{\\_}{k=1}^{K}\omega_{\\_}t^{k}$

return particles $\bm{x}_{\\_}{1:T}^{1:K}$, weights $\omega_{\\_}{1:T}^{1:K}$
and likelihood estimate $\hat{Z}_{\\_}{\text{HMSC}}$.

###  3.3 Objective Function

In this work, we choose the objective function to be the ELBO defined as SMC
marginal likelihood estimator in (2). However, we can show that incorporation
of HMC can make the ELBO objective arbitrarily tight. Defining
$\hat{Z}_{\\_}{\text{HSMC}}:=\prod_{\\_}{t=1}^{T}\frac{1}{K}\sum_{\\_}{k=1}^{K}\omega_{\\_}t^{k}$,
the ELBO objective can be formulated as

$\displaystyle\text{ELBO}_{\\_}{\text{HSMC}}(\theta,\phi,S,K,\epsilon,\bm{u}_{\\_}{1:T},\bm{y}_{\\_}{1:T})$
(12) $\displaystyle=\int
Q_{\\_}{\text{HSMC}}(\bm{x}^{S,1:K}_{\\_}{1:T},\bm{p}^{S,1:K}_{\\_}{2:T},\alpha_{\\_}{1:T-1}^{1:K};\bm{u}_{\\_}{1:T},\bm{y}_{\\_}{1:T})\bigg{(}\log\hat{Z}_{\\_}{\text{HSMC}}\bigg{)}d\bm{x}_{\\_}{1:T}^{S,1:K}d\bm{p}^{S,1:K}_{\\_}{2:T}d\alpha_{\\_}{1:T-1}^{1:K}$

where $Q_{\\_}{\text{HSMC}}$ is formulated as

$\displaystyle
Q_{\\_}{\text{HSMC}}(\bm{x}_{\\_}{S,1:T}^{1:K},\bm{p}_{\\_}{2:T}^{S,1:K},\alpha_{\\_}{1:T-1}^{1:K};\bm{u}_{\\_}{1:T},\bm{y}_{\\_}{1:T})$
$\displaystyle=\bigg{(}\prod_{\\_}{k=1}^{K}p_{\\_}{1}(\bm{x}_{\\_}1^{S,k}|\bm{u}_{\\_}1)\bigg{)}\bigg{(}\prod_{\\_}{t=2}^{T}\prod_{\\_}{k=1}^{K}q^{S}_{\\_}{t,\phi}(\bm{x}_{\\_}t^{S,k},\bm{p}_{\\_}t^{S,k}\big{|}\bm{u}_{\\_}t,\bm{y}_{\\_}{t-1},\tilde{\bm{x}}_{\\_}{t-1}^{\alpha^{k}_{\\_}{t-1}})\cdot\text{Discrete}(\alpha_{\\_}{t-1}^{k}\big{|}\omega_{\\_}{t-1}^{1:K})\bigg{)}$

where $\bm{x}_{\\_}t^{S,k},\bm{p}_{\\_}t^{S,k}$ are obtained as (11) and
$q^{S}_{\\_}{t,\phi}$ is the distribution of particle and momenta variables
after Hamiltonian dynamics.

We adopt stochastic gradient descent (Hoffman et al.,, 2013) to learn optimal
model and mass matrix parameters $\theta,\phi$. The gradient of objective
function can be derived as below, where observed and control data
$\bm{y}_{\\_}{1:T},\bm{u}_{\\_}{1:T}$ are omitted here,

$\displaystyle\nabla_{\\_}{\theta,\phi}\text{ELBO}_{\\_}{\text{HSMC}}$
$\displaystyle=$ $\displaystyle\nabla_{\\_}{\theta,\phi}\int
Q_{\\_}{\text{HSMC}}(\bm{x}^{S,1:K}_{\\_}{1:T},\bm{p}^{S,1:K}_{\\_}{2:T},\alpha_{\\_}{1:T-1}^{1:K})\bigg{(}\log\hat{Z}_{\\_}{\text{HSMC}}\bigg{)}d\bm{x}_{\\_}{1:T}^{S,1:K}d\bm{p}^{S,1:K}_{\\_}{2:T}d\alpha_{\\_}{1:T-1}^{1:K}$
$\displaystyle=$
$\displaystyle\int\nabla_{\\_}{\theta,\phi}Q_{\\_}{\text{HSMC}}(\bm{x}^{S,1:K}_{\\_}{1:T},\bm{p}^{S,1:K}_{\\_}{2:T},\alpha_{\\_}{1:T-1}^{1:K})\bigg{(}\log\hat{Z}_{\\_}{\text{HSMC}}\bigg{)}$
$\displaystyle+Q_{\\_}{\text{HSMC}}(\bm{x}^{S,1:K}_{\\_}{1:T},\bm{p}^{S,1:K}_{\\_}{2:T},\alpha_{\\_}{1:T-1}^{1:K})\nabla_{\\_}{\theta,\phi}\log\hat{Z}_{\\_}{\text{HSMC}}d\bm{x}_{\\_}{1:T}^{S,1:K}d\bm{p}^{S,1:K}_{\\_}{2:T}d\alpha_{\\_}{1:T-1}^{1:K}$
$\displaystyle=$ $\displaystyle\int
Q_{\\_}{\text{HSMC}}(\bm{x}^{S,1:K}_{\\_}{1:T},\bm{p}^{S,1:K}_{\\_}{2:T},\alpha_{\\_}{1:T-1}^{1:K})\bigg{[}\nabla_{\\_}{\theta,\phi}\log
Q_{\\_}{\text{HSMC}}(\bm{x}^{S,1:K}_{\\_}{1:T},\bm{p}^{S,1:K}_{\\_}{2:T},\alpha_{\\_}{1:T-1}^{1:K})\log\hat{Z}_{\\_}{\text{HSMC}}$
$\displaystyle+\nabla_{\\_}{\theta,\phi}\log\hat{Z}_{\\_}{\text{HSMC}}\bigg{]}d\bm{x}^{S,1:K}_{\\_}{2:T}d\bm{p}^{S,1:K}_{\\_}{1:T}d\alpha_{\\_}{1:T-1}^{1:K}$

In order to reduce the gradient variance, we ignore the first term in the
squared bracket above.

###  3.4 Theoretical Analysis

Based on the property of importance sampling (Murphy,, 2012), we can easily
show that at each time step $t$ the weight expression (10) is an unbiased
estimator of marginal likelihood likelihood of observations $\bm{y}_{\\_}t$
conditioned on $\bm{u}_{\\_}t$ and $\bm{y}_{\\_}{t-1}$. In this section, we
show that as step number $S$ increasing with particle number $K$ fixed, our
learning objective (12) can be arbitrarily close or converge to the marginal
log likelihood of observations. As (Le et al.,, 2017; Maddison et al.,, 2017)
we can assume the state space model has independent structure
$p_{\\_}{\theta}(\bm{x}_{\\_}{1:t-1}|\bm{y}_{\\_}{1:t})=p_{\\_}{\theta}(\bm{x}_{\\_}{1:t-1}|\bm{y}_{\\_}{1:t-1})$
for $t=2,\ldots,T$. It is a reasonable assumption since in online learning and
many practical applications, the future observations can not be known in
advance. Define conditional marginal likelihood of observations for each time
$t$ and particle $k$ as

$\int
f_{\\_}{\theta}(\bm{x}|\tilde{\bm{x}}_{\\_}{t-1}^{\alpha_{\\_}t^{k}},\bm{u}_{\\_}t,\bm{y}_{\\_}{t-1})g_{\\_}{\theta}(\bm{y}_{\\_}t|\bm{x})\mathcal{N}(\bm{p}|0,M(\bm{x}))d\bm{x}d\bm{p}:=p_{\\_}{\theta}(\bm{y}_{\\_}t|\tilde{\bm{x}}_{\\_}{t-1}^{\alpha_{\\_}t^{k}},\bm{u}_{\\_}t,\bm{y}_{\\_}{t-1})$

Then, if we can show the convergence of ELBO at each time $t$ and particle
$k$, the convergence of ELBO objective across all time in (12) can be shown.

According to the Hamilton energy defined in (6), the HMC adopted in algorithm
2 is ergodic with invariant distribution as

$p_{\\_}{\theta}(\bm{x},\bm{p}\big{|}\bm{u}_{\\_}t,\bm{y}_{\\_}{t-1},\tilde{\bm{x}}_{\\_}{t-1}^{\alpha^{k}_{\\_}{t-1}})\propto
f_{\\_}{\theta}(\bm{x}|\tilde{\bm{x}}_{\\_}{t-1}^{\alpha_{\\_}t^{k}},\bm{u}_{\\_}t,\bm{y}_{\\_}{t-1})g_{\\_}{\theta}(\bm{y}_{\\_}t|\bm{x})\mathcal{N}(\bm{p}|0,M(\bm{x}))$
(13)

which is the posterior of joint probabilities of emission, transition and
momenta distributions, or alternatively,

$p_{\\_}{\theta}(\bm{x},\bm{p}\big{|}\bm{u}_{\\_}t,\bm{y}_{\\_}{t-1},\tilde{\bm{x}}_{\\_}{t-1}^{\alpha^{k}_{\\_}{t-1}})=\frac{f_{\\_}{\theta}(\bm{x}|\tilde{\bm{x}}_{\\_}{t-1}^{\alpha_{\\_}t^{k}},\bm{u}_{\\_}t,\bm{y}_{\\_}{t-1})g_{\\_}{\theta}(\bm{y}_{\\_}t|\bm{x})\mathcal{N}(\bm{p}|0,M(\bm{x}))}{p_{\\_}{\theta}(\bm{y}_{\\_}t|\tilde{\bm{x}}_{\\_}t^{\alpha_{\\_}{t-1}^{k}},\bm{u}_{\\_}t,\bm{y}_{\\_}{t-1})}$
(14)

Based on properties of HMC (Neal et al.,, 2011), the particle distribution
will tend to invariant distribution in total variation, with the increase of
step number,i.e., for all $t$ and $k$,

$\lim_{\\_}{S\to\infty}\|q^{S}_{\\_}{t,\phi}(\bm{x},\bm{p})-p_{\\_}{\theta}(\bm{x},\bm{p}|\bm{u}_{\\_}t,\bm{y}_{\\_}{t-1},\tilde{\bm{x}}_{\\_}{t-1}^{\alpha^{k}_{\\_}{t-1}})\|_{\\_}{\text{TV}}=0$
(15)

Then based on the weight definition (10), we have, for each $t$ and $k$,

$\displaystyle\lim_{\\_}{S\to\infty}\int
q^{S}_{\\_}{t,\phi}(\bm{x},\bm{p})\log\big{(}\omega^{k}_{\\_}t\big{)}d\bm{x}d\bm{p}$
(16) $\displaystyle=$ $\displaystyle\int
q^{\infty}_{\\_}{t,\phi}(\bm{x},\bm{p})\log\bigg{(}\frac{f_{\\_}{\theta}(\bm{x}|\tilde{\bm{x}}_{\\_}{t-1}^{\alpha_{\\_}t^{k}},\bm{u}_{\\_}t,\bm{y}_{\\_}{t-1})g_{\\_}{\theta}(\bm{y}_{\\_}t|\bm{x})\mathcal{N}(\bm{p}|0,M(\bm{x}))}{q^{\infty}_{\\_}{t,\phi}(\bm{x},\bm{p})}\bigg{)}d\bm{x}d\bm{p}$
$\displaystyle=$ $\displaystyle\log
p_{\\_}{\theta}(\bm{y}_{\\_}t|\tilde{\bm{x}}_{\\_}t^{\alpha_{\\_}{t-1}^{k}},\bm{u}_{\\_}t,\bm{y}_{\\_}{t-1})$

where the second equality is due to (14) and (15). Due to the concavity of log
function, the ELBO objective (12) at each time $t$ can be lower bounded as

$\displaystyle\int\prod_{\\_}{k=1}^{K}q^{S}_{\\_}{t,\phi}(\bm{x}^{S,k}_{\\_}t,\bm{p}_{\\_}t^{S,k})\log\bigg{[}\frac{1}{K}\sum_{\\_}{k=1}^{k}\omega_{\\_}t^{k}\bigg{]}d\bm{x}_{\\_}t^{S,1:K}d\bm{p}_{\\_}t^{S,1:K}$
(17) $\displaystyle\geq$
$\displaystyle\int\prod_{\\_}{k=1}^{K}q^{S}_{\\_}{t,\phi}(\bm{x}^{S,k}_{\\_}t,\bm{p}_{\\_}t^{S,k})\bigg{[}\frac{1}{K}\sum_{\\_}{k=1}^{K}\log\omega^{k}_{\\_}t\bigg{]}d\bm{x}_{\\_}t^{S,1:K}d\bm{p}_{\\_}t^{S,1:K}$
$\displaystyle=$ $\displaystyle\int
q^{S}_{\\_}{t,\phi}(\bm{x},\bm{p})\log(\omega^{k}_{\\_}t)d\bm{x}d\bm{p}$

Combining (16) and (17) yields that at each time $t$ the ELBO objective
converges to conditional log likelihood of observations. Due to the
independent structure of latent state space, we can show that the overall
objective (12) converges to the log likelihood of observations
$p(\bm{y}_{\\_}{1:T})$ with the increase of step number $S$. The empirical
study tells us that the performance is good enough when $S$ is just around 10.

##  4 Improving Gaussian Process State Space Model

Gaussian Process State Space Models (GP-SSM) are a popular class of stochastic
SSMs (Frigola et al.,, 2013, 2014; Eleftheriadis et al.,, 2017; Doerr et al.,,
2018; Ialongo et al.,, 2018). In SSM, at each time step the system is taken to
evolve as a Markov chain by the transition function, mapping a latent state to
the next. By placing a Gaussian Process (GP) prior on the transition function,
we can obtain the Gaussian process state-space model (GP-SSM), which is a
fully Bayesian non-parametric treatment on modeling problem. It has many
advantages: 1) better uncertainty estimates based on the posterior of the
transition function; 2) avoiding overfitting with little amount of data; 3)
handling large amount of data without model saturation.

###  4.1 Preliminary

Same as standard SSM, we model the sequence of observations
$\bm{Y}=\\{\bm{y}_{\\_}t\\}_{\\_}{t=1}^{T}$ by a corresponding sequence of
latent states $\bm{X}=\\{\bm{x}_{\\_}t\\}_{\\_}{t=1}^{T}$, and
$\bm{x}_{\\_}t\in\mathbb{R}^{D_{\\_}x},\bm{y}_{\\_}t\in\mathbb{R}^{D_{\\_}y}$.
Here the state transition is assumed to be governed by a nonparametric
stochastic function $\bm{f}\in\mathbb{R}^{D_{\\_}x}$ following a GP prior.
Specifically, we have

$\bm{x}_{\\_}{t+1}\sim\mathcal{N}(\bm{f}(\bm{x}_{\\_}t),\bm{Q}),\hskip
15.0ptf_{\\_}d\sim\mathcal{G}\mathcal{P}(0,k_{\\_}d(\cdot,\cdot)),\hskip
15.0ptd=1,\ldots,D_{\\_}x$

where $\bm{Q}$ is the variance matrix and every $d$-th latent dimension has
its own GP function $f_{\\_}d$. The emission function is still parametric same
as (1). In order to reduce the computation complexity in learning GP, we adopt
the induced-inputs method (Snelson and Ghahramani,, 2006) and variational
sparse GP (Titsias,, 2009) in model learning, which achieves success in many
practical problems. For every GP function $f_{\\_}d$, we first introduce $P$
inducing GP targets $\bm{z}_{\\_}d=[z^{1}_{\\_}d,\ldots,z^{P}_{\\_}d]$ at
inducing GP inputs $\zeta_{\\_}d=[\zeta^{1}_{\\_}d,\ldots,\zeta^{P}_{\\_}d]$,
which are jointly Gaussian with the transition function $f_{\\_}d$. Then for
every latent dimension, the true GP predictive distribution can be
approximated by using the set inducing inputs and outputs as below, with
notation $d$ omitted,

$p(f^{*}_{\\_}d|\bm{x}^{*},\bm{f},\bm{X})\approx
p(f^{*}|\bm{x}^{*},\bm{z},\bm{\zeta}),\hskip
15.0ptp(\bm{z})=\mathcal{N}(\bm{z}|0,\bm{K}_{\\_}{\zeta,\zeta})$ (18)

where the covariance matrix $\bm{K}_{\\_}{\zeta,\zeta}$ with entries
$K_{\\_}{ij}=k(\zeta_{\\_}i,\zeta_{\\_}j)$. Following (1) we show the joint
distribution of GP-SSM for completeness,

$p(\bm{y}_{\\_}{1:T},\bm{x}_{\\_}{1:T},\bm{f}_{\\_}{2:T},\bm{z})=\bigg{[}\prod_{\\_}{t=1}^{T}g_{\\_}{\theta}(\bm{y}_{\\_}t|\bm{x}_{\\_}t)\bigg{]}p(\bm{x}_{\\_}1)p(\bm{z})\bigg{[}\prod_{\\_}{t=2}^{T}p(\bm{x}_{\\_}t|\bm{f}_{\\_}t)p(\bm{f}_{\\_}t|\bm{x}_{\\_}{t-1},\bm{u}_{\\_}t,\bm{y}_{\\_}{t-1},\bm{z})\bigg{]}$
(19)

where
$p(\bm{f}_{\\_}t|\hat{\bm{x}}_{\\_}{t-1},\bm{z})=\prod_{\\_}{d=1}^{D_{\\_}x}p(f_{\\_}{t,d}|\bm{x}_{\\_}{t-1},\bm{u}_{\\_}t,\bm{y}_{\\_}{t-1},\bm{z}_{\\_}d)$
and $\bm{z}:=[\bm{z}_{\\_}1,\ldots,\bm{z}_{\\_}{D_{\\_}x}]$. Here
$p(\bm{x}_{\\_}t|\bm{f}_{\\_}t)=\mathcal{N}(\bm{x}_{\\_}t|\bm{f}_{\\_}t,\text{diag}(\sigma^{2}_{\\_}{x,1},\ldots,\sigma^{2}_{\\_}{x,D_{\\_}x})$
and $p(\bm{x}_{\\_}1),p(\bm{z})$ are assumed to be Gaussian.

###  4.2 Motivation

As far as we know, all previous work on GP-SSM (Frigola et al.,, 2013, 2014;
Frigola-Alcade,, 2015; Eleftheriadis et al.,, 2017; Doerr et al.,, 2018;
Ialongo et al.,, 2018) assume the emission function to be a linear mapping
between latent state $\bm{x}_{\\_}t$ and the mean of observation
$\bm{y}_{\\_}t|\bm{x}_{\\_}t$. However, in many practical applications the
emission involves nonlinear function of $\bm{x}_{\\_}t$ (Gultekin and
Paisley,, 2017). In this case the posterior of latent states is non-Gaussian
and possibly multi-modal. So, the Gaussian distribution of transition function
cannot approximate the real posterior well enough.

In previous work (Frigola et al.,, 2013, 2014; Ialongo et al.,, 2018) the
variational distributions need to optimize parameter matrices
$\bm{A}_{\\_}t,\bm{b}_{\\_}t$ and $\bm{S}_{\\_}t$, which increases model and
learning complexity. Although authors in (Doerr et al.,, 2018) directly use
transition prior as variational approximation for latent posterior without
extra parameters, it cannot exploit information contained in the current
observations other than by adapting $q(\bm{u})$, which cannot handle cases
with high observation noise and long sequence length (Ialongo et al.,, 2018).

###  4.3 Variational Inference for GP-SSM with HSMC

Here we use HSMC to solve problems mentioned above. We first design the
variational distribution. Following the structure of real latent posterior,
the latent states and transition function are not factorized, i.e.,
$q(\bm{X},f)=q(\bm{X}|f)q(f)$. By introducing induced outputs $\bm{z}$, the
variational distribution can be defined as

$q(\bm{x}_{\\_}{1:T},\bm{f}_{\\_}{2:T},\bm{z})=q(\bm{x}_{\\_}1)\bigg{[}\prod_{\\_}{t=2}^{T}p(\bm{x}_{\\_}t|\bm{f}_{\\_}t)\prod_{\\_}{d=1}^{D_{\\_}x}p(f_{\\_}{x,d}|\hat{\bm{x}}_{\\_}{t-1},\bm{z}_{\\_}d)\bigg{]}\cdot\bigg{[}\prod_{\\_}{d=1}^{D_{\\_}x}q(\bm{z}_{\\_}d)\bigg{]}$
(20)

where
$\hat{\bm{x}}_{\\_}{t-1}:=(\bm{x}_{\\_}{t-1},\bm{u}_{\\_}t,\bm{y}_{\\_}{t-1})$,
and $q(\bm{x}_{\\_}1)$ and $q(\bm{z}_{\\_}d)$ are assumed to be Gaussian,
i.e., $q(\bm{z}_{\\_}d)=\mathcal{N}(\bm{\mu}_{\\_}d,\bm{\Sigma}_{\\_}d)$ where
$\bm{\mu}_{\\_}d\in\mathbb{R}^{P}$ and
$\bm{\Sigma}_{\\_}d\in\mathbb{R}^{P\times P}$. We can further simplify the
variational inference by integrating out induced-outputs $\bm{z}$ (Hensman et
al.,, 2013). Then following (19) the variational transition function can be
derived as

$q(\bm{x}_{\\_}t|\hat{\bm{x}}_{\\_}{t-1})=\int\bigg{[}p(\bm{x}_{\\_}t|\bm{f}_{\\_}t)p(\bm{f}_{\\_}t|\hat{\bm{x}}_{\\_}{t-1},\bm{z})\prod_{\\_}{d=1}^{D_{\\_}x}q(\bm{z}_{\\_}d)\bigg{]}d\bm{f}_{\\_}td\bm{z}=\mathcal{N}(\tilde{\bm{\mu}},\tilde{{\Sigma}})$

where
$\tilde{\bm{\mu}}=[\tilde{\mu}_{\\_}1,\ldots,\tilde{\mu}_{\\_}{D_{\\_}x}]$ and
$\tilde{{\Sigma}}=\text{diag}(\tilde{\sigma}_{\\_}1^{2},\ldots,\tilde{\sigma}_{\\_}{D_{\\_}x}^{2})$,
for $d=1,\ldots,D_{\\_}x$,

$\displaystyle\tilde{\mu}_{\\_}d$ $\displaystyle=$ $\displaystyle
k_{\\_}{\hat{\bm{x}}_{\\_}{t-1},\bm{\zeta}_{\\_}d}K^{-1}_{\\_}{\bm{\zeta}_{\\_}d,\bm{\zeta}_{\\_}d}\bm{\mu}_{\\_}d$
$\displaystyle\tilde{\sigma}^{2}_{\\_}d$ $\displaystyle=$ $\displaystyle
k_{\\_}{\hat{\bm{x}}_{\\_}{t-1},\hat{\bm{x}}_{\\_}{t-1}}-k_{\\_}{\hat{\bm{x}}_{\\_}{t-1},\bm{\zeta}_{\\_}d}K^{-1}_{\\_}{\bm{\zeta}_{\\_}d,\bm{\zeta}_{\\_}d}(K_{\\_}{\bm{\zeta}_{\\_}d,\bm{\zeta}_{\\_}d}-\Sigma_{\\_}d)K^{-1}_{\\_}{\bm{\zeta}_{\\_}d,\bm{\zeta}_{\\_}d}k_{\\_}{\hat{\bm{x}}_{\\_}{t-1},\bm{\zeta}_{\\_}d}^{T}+\sigma^{2}_{\\_}{x,d},\hskip
15.0pt$

In this work, we use $K$-particle HSMC to approximate the intractable
posterior. Denote
$\hat{\bm{x}}_{\\_}{t-1}^{k}=(\bm{x}_{\\_}{t-1}^{\alpha^{k}_{\\_}{t-1}},\bm{u}_{\\_}t,\bm{y}_{\\_}{t-1})$
where $\alpha_{\\_}{t-1}^{k}$ is the re-sampling index. At each time step $t$,
we first sample multiple particles of latent states $\bm{x}_{\\_}t^{0,k}$ and
momenta $\bm{p}_{\\_}t^{0,k}$ from the transition function
$q(\bm{x}_{\\_}t|\hat{\bm{x}}^{k}_{\\_}{t-1})$ and Gaussian
$\mathcal{N}(0,\bm{M}(\bm{x}_{\\_}t^{0,k}))$ respectively. Then by using
$S$-step RMHMC operation we transform the sampled particles to approximate the
posterior as

$p_{\\_}{\theta}(\bm{x},\bm{p}|\hat{\bm{x}}_{\\_}{t-1})\propto
g_{\\_}{\theta}(\bm{y}_{\\_}t|\bm{x})q(\bm{x}|\hat{\bm{x}}_{\\_}{t-1})\mathcal{N}(\bm{p}|0,\bm{M}(\bm{x}))$

Due to the measure preserving and reversibility property of HMC, we define the
particle weights as

$\omega_{\\_}t^{k}=\frac{g_{\\_}{\theta}(\bm{y}_{\\_}t|\bm{x}_{\\_}t^{S,k})q(\bm{x}_{\\_}t^{S,k}|\hat{\bm{x}}_{\\_}{t-1}^{k})\mathcal{N}(\bm{p}_{\\_}t^{S,k}|0,\bm{M}(\bm{x}^{S,k}_{\\_}t))}{q(\bm{x}_{\\_}t^{0,k}|\hat{\bm{x}}_{\\_}{t-1}^{k})\mathcal{N}(\bm{p}_{\\_}t^{0,k}|0,\bm{M}(\bm{x}_{\\_}t^{0,k}))}$
(21)

According to the joint and variational distribution in (19)(20), the ELBO can
be computed as

$\displaystyle\text{ELBO}_{\\_}{\text{GP-SSM}}$ $\displaystyle=$
$\displaystyle\sum_{\\_}{t=1}^{T}\int\prod_{\\_}{k=1}^{K}q(\bm{x}_{\\_}t^{S,k}|\hat{\bm{x}}_{\\_}{t-1}^{k})\log\bigg{(}\frac{1}{K}\sum_{\\_}{k=1}^{K}\omega_{\\_}t^{k}\bigg{)}d\bm{x}_{\\_}t^{0,1:K}$
$\displaystyle-\sum_{\\_}{d=1}^{D_{\\_}x}\text{KL}\big{(}q(\bm{z}_{\\_}d)\|p(\bm{z}_{\\_}d))\big{)}$

where $p(\bm{z}_{\\_}d)$ is defined in (18). In experimental study, we show
that without increasing model complexity, HSMC can improve learning
performance on GP-SSM with nonlinear emission.

##  5 Experiments

In this section, we conduct experiments to show the benefits of proposed
algorithms. First, we use synthetic data to verify the advantages of HSMC over
conventional variational SMC, where both transition and emission are realized
by neural networks. Then based on real-world dataset, we show that HSMC can
better learn GP-SSM when emission function is nonlinear.

###  5.1 Synthetic Data

The synthetic data is generated by a nonlinear state-space model

$\bm{x}_{\\_}{t+1}=\bm{A}\bm{x}_{\\_}t+\bm{\epsilon}_{\\_}t,\hskip
20.0pt\bm{y}_{\\_}t=g_{\\_}{\theta}(\bm{x}_{\\_}t)+\bm{\xi}_{\\_}t$

where $\bm{x}_{\\_}t\in\mathbb{R}^{10},\bm{y}_{\\_}t\in\mathbb{R}^{30}$ and
$g_{\\_}{\theta}$ is realized by two-layer neural network with 20 hidden
neurons. The non-linearities in $g_{\\_}{\theta}$ are ReLU and Sigmoid. And
the noise $\bm{\epsilon}_{\\_}t,\bm{\xi}_{\\_}t$ follow independent Gaussian
with variance of $0.2$. Every sequence has length of 100. Both training and
testing datasets only contain observations $\bm{y}_{\\_}t$.

In this experiment, the base models are variational recurrent neural network
(VRNN) (Chung et al.,, 2015) and stochastic recurrent neural network (SRNN)
(Fraccaro et al.,, 2016), where the emission is realized by deep neural
network and transition is implemented by gated recurrent unit (GRU) (Chung et
al.,, 2014). The comparison is between the base models (VRNN or SRNN)
augmented by FIVO (Maddison et al.,, 2017) and HSMC. And generation network in
every compared model is the same. However, in base model with HSMC, the
proposal neural network is omitted. The performance metric is the log-
likelihood per step. Here $K$ is the number of particles and $S$ is the number
of step in RMHMC. Performance of FIVO is not related with $S$.

Table 1: Performance Comparison in Synthetic Data | $K=10,S=5$ | $K=10,S=10$ | $K=5,S=5$ | $K=5,S=10$
---|---|---|---|---
VRNN-FIVO | 2.56 | 2.56 | 2.12 | 2.12
VRNN-HSMC | 2.78 | 2.86 | 2.37 | 2.54
SRNN-FIVO | 2.64 | 2.64 | 2.31 | 2.31
SRNN-HSMC | 2.80 | 2.86 | 2.49 | 2.53

###  5.2 Bike-sharing Demand Data

In this experiment, we use bike-sharing record data in New York city
(dataset), from July 2014 to July 2017. In every transaction record, there are
trip duration, bike check out/in time, names of start and end stations, and
customer information such as age and gender. We first aggregate transaction
records to bike demands at each station, ignoring customer information. Then
we further remove stations existing for less than two years. And stations with
less than one bike used per hour are also deleted. Finally, there are only 269
stations left in the dataset.

Here we compare the orginal GP-SSM (Doerr et al.,, 2018) with multiple
particles (GP-SSM-SMC) and that augmented by HSMC (GP-SSM-HSMC). Both models
have same latent dimension and number of inducing points. The performance
metric is log likelihood. The results are shown below. We also the results of
VRNN-FIVO for comparison.

Table 2: Performance Comparison in Bike-sharing Demand Data | $K=10,S=5$ | $K=10,S=10$ | $K=5,S=5$ | $K=5,S=10$
---|---|---|---|---
GP-SSM-SMC | -653.2 | -653.2 | -668.7 | -668.7
GP-SSM-HSMC | -607.1 | -591.3 | -621.4 | -610.9
VRNN-FIVO | -589.2 | -589.2 | -599.1 | -599.1

We find that HSMC can help diminish the performance gap between GP model and
neural network model, even though GP has much less parameters than neural
networks.

## References

  * Berkenkamp et al., (2017) Berkenkamp, F., Turchetta, M., Schoellig, A., and Krause, A. (2017).  Safe model-based reinforcement learning with stability guarantees.  In Advances in Neural Information Processing Systems, pages 908–918.
  * Billings, (2013) Billings, S. A. (2013).  Nonlinear system identification: NARMAX methods in the time, frequency, and spatio-temporal domains.  John Wiley & Sons.
  * Burda et al., (2015) Burda, Y., Grosse, R., and Salakhutdinov, R. (2015).  Importance weighted autoencoders.  arXiv preprint arXiv:1509.00519.
  * Camacho and Alba, (2013) Camacho, E. F. and Alba, C. B. (2013).  Model predictive control.  Springer Science & Business Media.
  * Chung et al., (2014) Chung, J., Gulcehre, C., Cho, K., and Bengio, Y. (2014).  Empirical evaluation of gated recurrent neural networks on sequence modeling.  arXiv preprint arXiv:1412.3555.
  * Chung et al., (2015) Chung, J., Kastner, K., Dinh, L., Goel, K., Courville, A. C., and Bengio, Y. (2015).  A recurrent latent variable model for sequential data.  In Advances in neural information processing systems, pages 2980–2988.
  * Deisenroth and Rasmussen, (2011) Deisenroth, M. and Rasmussen, C. E. (2011).  Pilco: A model-based and data-efficient approach to policy search.  In Proceedings of the 28th International Conference on machine learning (ICML-11), pages 465–472.
  * Doerr et al., (2018) Doerr, A., Daniel, C., Schiegg, M., Nguyen-Tuong, D., Schaal, S., Toussaint, M., and Trimpe, S. (2018).  Probabilistic recurrent state-space models.  arXiv preprint arXiv:1801.10395.
  * Eleftheriadis et al., (2017) Eleftheriadis, S., Nicholson, T., Deisenroth, M., and Hensman, J. (2017).  Identification of gaussian process state space models.  In Advances in Neural Information Processing Systems, pages 5309–5319.
  * Fraccaro et al., (2016) Fraccaro, M., Sønderby, S. K., Paquet, U., and Winther, O. (2016).  Sequential neural models with stochastic layers.  In Advances in neural information processing systems, pages 2199–2207.
  * Frigola et al., (2014) Frigola, R., Chen, Y., and Rasmussen, C. E. (2014).  Variational gaussian process state-space models.  In Advances in Neural Information Processing Systems, pages 3680–3688.
  * Frigola et al., (2013) Frigola, R., Lindsten, F., Schön, T. B., and Rasmussen, C. E. (2013).  Bayesian inference and learning in gaussian process state-space models with particle mcmc.  In Advances in Neural Information Processing Systems, pages 3156–3164.
  * Frigola-Alcade, (2015) Frigola-Alcade, R. (2015).  Bayesian time series learning with gaussian processes.  Uni-versity of Cambridge.
  * Ghahramani and Roweis, (1999) Ghahramani, Z. and Roweis, S. T. (1999).  Learning nonlinear dynamical systems using an em algorithm.  In Advances in neural information processing systems, pages 431–437.
  * Girolami and Calderhead, (2011) Girolami, M. and Calderhead, B. (2011).  Riemann manifold langevin and hamiltonian monte carlo methods.  Journal of the Royal Statistical Society: Series B (Statistical Methodology), 73(2):123–214.
  * Gordon et al., (1993) Gordon, N. J., Salmond, D. J., and Smith, A. F. (1993).  Novel approach to nonlinear/non-gaussian bayesian state estimation.  In IEE Proceedings F-radar and signal processing, volume 140, pages 107–113. IET.
  * Gultekin and Paisley, (2017) Gultekin, S. and Paisley, J. (2017).  Nonlinear kalman filtering with divergence minimization.  IEEE Transactions on Signal Processing, 65(23):6319–6331.
  * Hensman et al., (2013) Hensman, J., Fusi, N., and Lawrence, N. D. (2013).  Gaussian processes for big data.  arXiv preprint arXiv:1309.6835.
  * Hoffman et al., (2013) Hoffman, M. D., Blei, D. M., Wang, C., and Paisley, J. (2013).  Stochastic variational inference.  The Journal of Machine Learning Research, 14(1):1303–1347.
  * Ialongo et al., (2018) Ialongo, A. D., van der Wilk, M., Hensman, J., and Rasmussen, C. E. (2018).  Non-factorised variational inference in dynamical systems.  arXiv preprint arXiv:1812.06067.
  * Kingma and Welling, (2013) Kingma, D. P. and Welling, M. (2013).  Auto-encoding variational bayes.  arXiv preprint arXiv:1312.6114.
  * Le et al., (2017) Le, T. A., Igl, M., Rainforth, T., Jin, T., and Wood, F. (2017).  Auto-encoding sequential monte carlo.  arXiv preprint arXiv:1705.10306.
  * Leimkuhler and Reich, (2004) Leimkuhler, B. and Reich, S. (2004).  Simulating hamiltonian dynamics, volume 14.  Cambridge university press.
  * Ljung, (1998) Ljung, L. (1998).  System identification.  In Signal analysis and prediction, pages 163–173. Springer.
  * Ljung, (2010) Ljung, L. (2010).  Perspectives on system identification.  Annual Reviews in Control, 34(1):1–12.
  * Maddison et al., (2017) Maddison, C. J., Lawson, J., Tucker, G., Heess, N., Norouzi, M., Mnih, A., Doucet, A., and Teh, Y. (2017).  Filtering variational objectives.  In Advances in Neural Information Processing Systems, pages 6573–6583.
  * Murphy, (2012) Murphy, K. P. (2012).  Machine learning: A probabilistic perspective. adaptive computation and machine learning.
  * Naesseth et al., (2017) Naesseth, C. A., Linderman, S. W., Ranganath, R., and Blei, D. M. (2017).  Variational sequential monte carlo.  arXiv preprint arXiv:1705.11140.
  * Neal et al., (2011) Neal, R. M. et al. (2011).  Mcmc using hamiltonian dynamics.  Handbook of Markov Chain Monte Carlo, 2(11):2.
  * Raiko and Tornio, (2009) Raiko, T. and Tornio, M. (2009).  Variational bayesian learning of nonlinear hidden state-space models for model predictive control.  Neurocomputing, 72(16-18):3704–3712.
  * Rezende et al., (2014) Rezende, D. J., Mohamed, S., and Wierstra, D. (2014).  Stochastic backpropagation and approximate inference in deep generative models.  In International Conference on Machine Learning, pages 1278–1286.
  * Snelson and Ghahramani, (2006) Snelson, E. and Ghahramani, Z. (2006).  Sparse gaussian processes using pseudo-inputs.  In Advances in neural information processing systems, pages 1257–1264.
  * Titsias, (2009) Titsias, M. (2009).  Variational learning of inducing variables in sparse gaussian processes.  In Artificial Intelligence and Statistics, pages 567–574.
  * Wan and Nelson, (1997) Wan, E. A. and Nelson, A. T. (1997).  Dual kalman filtering methods for nonlinear prediction, smoothing and estimation.  In Advances in neural information processing systems, pages 793–799.


# Particle-field duality in QFT measurements

Maria Papageorgiou<EMAIL_ADDRESS>Department of Physics,
University of Patras, 26500 Greece  Department of Applied Mathematics,
University of Waterloo, Waterloo, Ontario, N2L 3G1, Canada  Institute for
Quantum Computing, University of Waterloo, Waterloo, Ontario, N2L 3G1, Canada
Jose de Ramón<EMAIL_ADDRESS>Departamento de Fisica, Universidad de Burgos
Pza. Misael Bañuelos s.n., 09001-Burgos, Spain  Charis Anastopoulos
<EMAIL_ADDRESS>Department of Physics, University of Patras, 26500 Greece

###### Abstract

Pointlike systems coupled to quantum fields are often employed as toy models
for measurements in quantum field theory. In this paper, we identify the field
observables recorded by such models. We show that in models that work in the
strong coupling regime, the apparatus is correlated with smeared field
amplitudes, while in models that work in weak coupling the apparatus records
particle aspects of the field, such as the existence of a particle-like time
of arrival and resonant absorption. Then, we develop an improved field-
detector interaction model, adapting the formalism of Quantum Brownian motion,
that is exactly solvable. This model confirms the association of field and
particle properties in the strong and weak coupling regimes, respectively.
Further, it can also describe the intermediate regime, in which the field-
particle characteristics ‘merge’. In contrast to standard perturbation
techniques, this model also recovers the relativistic Breit-Wigner resonant
behavior in the weak coupling regime. The modulation of field-particle-duality
by a single tunable parameter is a novel feature that is, in principle,
experimentally accessible.

##  I Introduction

Quantum field theory (QFT) has proved immensely successful in describing the
three fundamental interactions (weak, strong, and electromagnetic). Despite
the severe challenges in its rigorous mathematical formulation, theorists have
developed pragmatic tools for calculations that work really well, and provide
unambiguous predictions to a large class of experiments. There is little
reason to expect a failure of QFT at high energies, as long as one does not
attempt to deal with gravity.

For this reason, it is really astonishing that the development of a consistent
measurement theory for local measurements in QFT is an ongoing project until
today. This is not apparent in high-energy applications that focus on
scattering experiments that can be treated via S-matrix theory (see discussion
in [1]). The latter conceptualizes scattering as a process with a single
measurement event in the far future, so there is no reason to employ the
state-update rule. This rule is employed, for example, in measurements of
Bell-type correlations, sequential measurements, and post-selected
measurements. The problem is that the usual state-update rule conflicts with
relativistic causality—see, [3, 4, 2] and references therein. For this reason,
even axiomatizations of QFT [7, 11] refrained from including a state-update
rule into their axioms (see discussion in [5]).

The state-update rule is essentially the rule for joint probabilities for
multiple measurement events, and no probabilistic theory can work in its
absence. Many experiments in quantum optics require the use of such
probabilities, for example, to describe photon bunching and anti-bunching.
These came to be described by photo-detection theories, the most prominent of
which is that of Glauber [13].

Glauber’s theory was arguably the first QFT measurement theory, and it proved
remarkably successful. However, its scope is limited. It works only for
photons, and it involves approximations that may not be appropriate in set-ups
of current interest, for example, quantum experiments in space that involve
long-distance propagation of photons [14, 15]. These approximations can
introduce friction with relativistic causality. More generally, the extension
of quantum information theory to the relativistic regime, in order to take
into account the effects of gravity and motion, arguably requires and
consistent a practicable theory of QFT measurements [4, 16].

A consistent and practicable measurement theory should fulfill some
requirements. A tentative and non-exhaustive list is the following:

  1. 1.

Causality: Probabilities are compatible with relativistic causality, that is,
there is no signaling outside the lightcone.

  2. 2.

Covariance: The probability assignment does not depend on the choice of
coordinate system. However, it may depend on the state of motion of measuring
apparatuses.

  3. 3.

Localization: The theory must provide a pragmatic definition of spatially
localized systems that is compatible with experimental practice. This would
lead to the notion of multi-partite systems and to the construction of the
associated probabilities.

  4. 4.

Non-relativistic limit: The probability assignment recovers the standard non-
relativistic measurement theory in an appropriate limit.

  5. 5.

Practicality: It should provide concrete predictions for experiments that
cannot be modeled by S-matrix theory, for example, in quantum optics.

A theory fulfilling all these requirements is still missing. Progress has been
made on several fronts, but a recurring feature is that when considering
relativistic measurement models the apparatus performing the measurement plays
a more prominent role than in the non-relativistic counterpart.

There are two paradigms that address these challenges, with varying degrees of
success. First, histories-based approaches and, second, local generalizations
of scattering processes.

An example of a histories-based approach is the Quantum Temporal Probabilities
(QTP) approach [22, 23, 24, 4]. This is based on the decoherent histories
approach to quantum measurement that emphasizes emerging classicality in the
measuring apparatus [25, 26, 27, 28, 29]. QTP emphasizes the spacetime aspects
of QFT probabilities, namely, that every measurement record is localized in
spacetime, and that physical probabilities must correspond to genuine
densities in spacetime.

The second approach essentially formulates a QFT version of von Neumann’s
measurement theory [17, 6] that involves switching on an interaction between
the measured system and a probe in a compact region of spacetime. In this
paradigm, sequential measurements are envisioned as sequential implementations
of interactions with independent probes, all of which interact locally with
the field. Therefore, these approaches deal with a generalization of particle
physics-like scattering for modeling sequential measurements. When the probe
is modeled as a relativistic quantum field, as in [18, 19, 20, 21], and such
that the interaction between apparatus and field satisfies some locality
conditions, covariance and some versions of causality are guaranteed.

Another approach to QFT measurements that falls into the scattering-like
paradigm, is the detector model approach. It is a pragmatic approach that is
based on modeling the coupling of the quantum field to pointlike (or extended)
quantum-mechanical systems that follow prescribed spacetime trajectories in
flat or curved spacetime [31]. Standard measurement theory is applied to the
non-relativistic detector, and the induced state-update ruled for the field
can be derived [30]. Frictions with causality and covariance arise [32, 33,
34, 35] and can be handled by specifying the regime of validity of the model.

Detector models are a simple and useful tool for exploring many issues
relevant to information and measurements in QFT [36]. The main example of such
models is the Unruh-DeWitt particle detector [37, 38], originally employed in
the discussion of the Unruh effect. Particle detector models, as their name
suggests, were originally intended to tackle the challenges associated with
the notion of particle in relativistic quantum field theory (QFT), that is, to
deal with the complications that a naive notion of particle introduces in
scenarios involving non-inertial observers or curved spacetimes. Certainly,
many authors argue that QFT does not admit a particle ontology [8, 9, 10]. In
addition, there is controversy around field ontology (see [12] and references
therein). In this paper, we are agnostic about this issue. We make no claim
that particle detector models are fundamental, rather we view them as useful
tools for extracting QFT phenomenology in certain regimes (see discussion in
[10]).

This paper follows the detector model approach to QFT measurements. Our
original motivation was to answer the question: What do detectors detect? Or
equivalently: Which field observables are measured by pointlike detectors? To
this end, we considered two models, one that works well in the strong coupling
regime and one that is commonly used in the weak coupling regime. The first
model is the von Neumann model of measurement [17] with an apparatus that
consists of one localized degree of freedom that is strongly coupled to a
scalar field and it behaves as a pointer variable. The second model is the
Unruh-DeWitt detector coupled to the scalar field, and treated to lowest order
in perturbation theory. For both models, we analyze in detail the limit in
which the procedure of switching the detector on and off is irrelevant, that
is, when the detector and the field interact freely without external
influences. We found that in the von Neumann model, the pointer variable is
correlated with spacetime averaged field amplitude, and we analyzed the exact
sense in which in the Unruh-DeWitt model, the detector captures the particle
aspects of the field at the weak-coupling limit.

Our results suggest that the coupling strength is a key factor in modulating
field-particle duality in field-detector interactions. This property,
pertaining to a fundamental aspect of QFT, could be measured in any set-up
that allows for a tunable interaction strength of a localized system with a
quantum field, for example, in circuit quantum electrodynamics [41], cavity
optomechanics [42] or in analogue gravity experiments with Bose-Einstein
condensates [43, 44].

For this reason, we constructed a new QFT measurement model that enables us to
analyze both strong and weak coupling, but also the interpolating regimes.
This model is an adaptation of the quantum Brownian motion (QBM) models, that
have been extensively studied in the theory of open quantum systems [46, 47,
48]. QBM models describe a harmonic oscillator in contact with a heat bath,
itself modeled by a large number of harmonic oscillators. With some
modifications that model can describe the measurement of a quantum field by an
oscillator detector.

The key point is that our model can be solved exactly, so it goes well beyond
the results obtained in the lowest order in perturbation theory in the Unruh-
DeWitt model, or through the drastic approximations of von Neumann’s model. We
show that, indeed, the detector captures the field aspects at strong coupling
and particle aspects, such as resonance and time-of-arrival behavior, at weak
coupling. In the regime of intermediate couplings, field and particle
characteristics ”merge”. It is an advantage of solvable models that they can
be used to characterize intermediate regimes, in which neither the field nor
the particle concepts can be used exclusively for interpreting the detector’s
response.

The structure of this paper is the following. In Sec. 2, we analyze field
measurements in von Neumann’s model, and in Sec. 3., we explain in what sense
the Unruh-DeWitt model acts as a particle detector. In Sec. 4, we develop our
QBM-based measurement model, and we analyze its properties. In Sec. 5, we
summarize and discuss our results.

##  II Field measurements in von Neumann’s model

The simplest quantum measurement model involves the recording of an observable
$\hat{A}$ by a single pointer variable $\hat{X}$ with continuous spectrum. Let
${\cal H}_{S}$ be the Hilbert space of the microscopic system that is being
measured and ${\cal H}_{A}$ be the Hilbert space of the apparatus. Let the
Hamiltonian of the system be $\hat{H}_{S}$ and the Hamiltonian of the
apparatus $\hat{H}_{A}$. In an ideal apparatus, the pointer variable only
changes during the interaction with the microscopic system. Once this
interaction is over, the pointer variable $\hat{X}$ is frozen. This means that
$[\hat{H}_{A},\hat{X}]=0$. The essence of this analysis remains unchanged if
we take $\hat{H}_{A}=0$, or if we assume that $\hat{H}_{A}$ is negligible in
comparison to the interaction Hamiltonian.

In measurement models, it is common to assume that the interaction between
system and apparatus is switched on for a finite time. We consider an
interaction Hamiltonian of the form

$\hat{H}_{int}=f(t)\hat{A}\otimes\hat{P},$ (1)

where $\hat{P}$ is the ‘conjugate momentum’ of the pointer variable $\hat{X}$,
i.e., $[\hat{X},\hat{P}]=i$. The switching function $f(t)$ vanishes outside
the time interval $[0,T]$.

Let $|\Omega\rangle$ be the initial state of the apparatus, and
$|\psi_{0}\rangle$ the initial state of the measured system. It is
straightforward to calculate that the probability density $p(X)$ for the
pointer variable at times $t>T$

$\displaystyle
p(X)=\langle\psi_{0}|\hat{D}^{\dagger}_{X}\hat{D}_{X}|\psi_{0}\rangle,$ (2)

where

$\displaystyle\hat{D}_{X}=\int\frac{\mathrm{d}k}{\sqrt{2\pi}}e^{ikX}\langle
k|\Omega\rangle{\cal T}e^{-ik\int\mathrm{d}tf(t)\hat{A}(t)}$ (3)

is the measurement operator defined in terms of the eigenvectors of $\hat{P}$
denoted as $|k\rangle$ . In Eq. (3), ${\cal T}$ stands for the time-ordered
expansion of the exponential and $\hat{A}(t)$ is the Heisenberg-picture
operator $e^{i\hat{H}_{S}t}\hat{A}e^{-i\hat{H}_{S}t}$. It can be shown that
the operator

$\displaystyle\hat{\Pi}_{X}=\hat{D}^{\dagger}_{X}\hat{D}_{X}$ (4)

defines a Positive-Operator-Valued measure for the time-averaged values of the
observable $\hat{A}$, weighted by the function $f(t)$ [45].

Here we apply this model to the measurement of a scalar field $\hat{\phi}$. We
take
$\hat{A}=\hat{\phi}(g)=\int\mathrm{d}^{3}x\hat{\phi}({\bm{x}})g({\bm{x}})$ for
some smearing function $g({\bm{x}})$ that models the spatial extension of the
interaction. We work in a reference frame where the detector is located at the
coordinate origin, so $g({\bm{x}})$ is peaked around ${\bm{x}}=0$.

For a free field of mass $m$,

$\displaystyle\hat{H}_{S}=\frac{1}{2}\int\mathrm{d}^{3}x\left(\hat{\pi}^{2}+(\nabla\hat{\phi})^{2}+m^{2}\hat{\phi}^{2}\right).$
(5)

If we denote
$\hat{\phi}_{t}(g)=e^{i\hat{H}_{S}t}\hat{\phi}(g)e^{-i\hat{H}_{S}t}$ The field
commutation relations are

$\displaystyle[\hat{\phi}_{t_{1}}(g),\hat{\phi}_{t_{2}}(g)]$
$\displaystyle=\mathrm{i}\int\mathrm{d}^{3}x\mathrm{d}^{3}x^{\prime}\Delta(t_{1}-t_{2},{\bm{x}}-{\bm{x}}^{\prime})g({\bm{x}})g({\bm{x}}^{\prime}),$
(6)

where $\epsilon_{\bm{p}}=\sqrt{{\bm{p}}^{2}+m^{2}}$, and

$\displaystyle\Delta(t,{\bf
x})=\int\frac{\mathrm{d}^{3}p}{(2\pi)^{3}\epsilon_{\bm{p}}}\sin({\bm{p}}\cdot{\bm{x}}-\epsilon_{\bm{p}}t)$
(7)

is the commutator distribution for the Klein-Gordon field.

Consider a pointer variable is locally coupled to the field through the
interaction Hamiltonian

$\displaystyle\hat{H}_{int}=f(t)\hat{\phi}(g)\otimes\hat{P}.$ (8)

We define the spacetime smearing function $F(t,{\bm{x}})=f(t)g({\bm{x}})$, and
we write the associated smeared Heisenberg-picture field
$\hat{\phi}(F)=\int\mathrm{d}t\mathrm{d}\bm{x}F(t,{\bm{x}})\hat{\phi}(t,{\bm{x}})$.
Then, we find that the probability density of the pointer variable after the
interaction with the field is given by equation (2) with

$\displaystyle\hat{D}_{X}=w[X-\hat{\phi}(F)],$ (9)

in terms of

$\displaystyle
w(x)=\int\frac{\mathrm{d}k}{\sqrt{2\pi}}e^{ikx-\frac{i}{2}k^{2}\delta_{F}}\langle
k|\Omega\rangle,$ (10)

where

$\displaystyle\delta_{F}=\int\mathrm{d}t$
$\displaystyle\mathrm{d}^{3}x\mathrm{d}t^{\prime}\mathrm{d}^{3}x^{\prime}\theta(t-t^{\prime})$
$\displaystyle\times
F(t,{\bm{x}})\Delta(t-t^{\prime},{\bm{x}}-{\bm{x}}^{\prime})F(t^{\prime},{\bm{x}}^{\prime}).$
(11)

To derive (10) we have used the the Magnus expansion for the time-ordered
exponential that appears in (3), namely,

$\displaystyle{\cal
T}e^{-ik\int\mathrm{d}tf(t)\hat{\phi}_{t}(g)}=\exp\left[-ik\hat{\phi}(F)-\frac{i}{2}k^{2}\delta_{F}\right].$
(12)

It is straightforward to check that $\int\mathrm{d}x|w(x)|^{2}=1$, hence, that
Eq. ((4)) defines a Positive-Operator-Valued-Measure (POVM). This POVM
corresponds to the measurement of time-averaged smeared-field amplitudes.
Indeed, we will associate the expectation value and variance of the pointer
variable after the interaction to the expectation value and variance of the
smeared field amplitudes.

To this end, we choose an initial state $|\Omega\rangle$ for the pointer
variables, such that $\langle\Omega|\hat{X}|\Omega\rangle=0$ and
$\langle\Omega|\hat{P}|\Omega\rangle=0$, and no initial correlation between
position and momentum. Then, $\int\mathrm{d}xx|w(x)|^{2}=0$, and the variance
${\cal N}^{2}$ of $|w(x)|^{2}$ is

$\displaystyle{\cal N}^{2}=\sigma_{X}^{2}+\sigma_{P}^{2}\delta_{F}^{2},$ (13)

where $\sigma^{2}_{X}$ and $\sigma^{2}_{P}$ are the position and momentum
variances for $|\Omega\rangle$, respectively.

It follows that

$\displaystyle\langle\hat{X}\rangle=\langle\psi_{0}|\hat{\phi}(F)|\psi_{0}\rangle,$
(14) $\displaystyle(\Delta X)^{2}=[\Delta\hat{\phi}(F)]^{2}+{\cal N}^{2},$
(15)

where $\psi_{0}$ the initial state of the field. This implies that the pointer
variable $\hat{X}$ is correlated with the smeared field $\hat{\phi}(F)$. The
correlation is not perfect, as it is characterized by a noise ${\cal N}$ that
is intrinsic to the detector. By the uncertainty relation
$\sigma_{X}\sigma_{P}\geq\frac{1}{2}$, we derive a lower bound for ${\cal N}$,

$\displaystyle{\cal N}\geq\sqrt{|\delta_{F}|}.$ (16)

This bound does not depend on the initial state of the apparatus, but only on
the localization area of the field, and the strength of the coupling as
encoded in $\delta_{F}$.

The signal-to-noise ratio for this system is $|\langle\hat{X}\rangle|/{\cal
N}$. Hence, the necessary condition for a well-defined measurement is that

$\displaystyle|\langle\hat{\phi}(F)\rangle|>>\sqrt{|\delta_{F}|}$ (17)

For a state with a fixed number of particles, e.g., a wavepacket, it holds
that $\langle\hat{\phi}(F)\rangle=0$ and as a result there is no signal. This
implies that this model cannot capture the particle aspects of the quantum
field.

To estimate the strength of the noise in a field amplitude measurement,
consider a Gaussian smearing function

$\displaystyle F(t,{\bf
x})=\frac{\lambda}{4\pi^{2}s_{X}^{3}s_{T}}e^{-\frac{{\bf
x}^{2}}{2s_{X}^{2}}-\frac{t^{2}}{2s_{T}^{2}}},$ (18)

with spatial width $s_{X}$ and temporal width $s_{T}$; $\lambda$ plays the
role of the coupling constant. The function $F$ becomes proportional to a
delta function for $s_{X},s_{T}\rightarrow 0$, while the limit of a constant
interaction corresponds to $\sigma_{T},\lambda\rightarrow\infty$, with
$\lambda/\sigma_{T}$ constant. We straightforwardly calculate

$\displaystyle\delta_{F}=\frac{2\lambda^{2}}{\pi^{7/2}}\int_{0}^{\infty}\frac{p^{2}\mathrm{d}p}{\epsilon_{p}}e^{-s_{X}^{2}p^{2}}D(s_{T}\epsilon_{p}),$
(19)

where $D(x):=e^{-x^{2}}\int_{0}^{x}dte^{t^{2}}$ is the Dawson function. For
$m\neq 0$, and $ms_{T}>>1$, we can use the asymptotic expression for the
Dawson function $D(x)\simeq\frac{1}{2x}$, to obtain
$\delta_{F}\sim\sigma_{x}^{-3}\sigma_{T}^{-1}$, i.e., $\delta_{F}$ is
inversely proportional to the spacetime volume of the interaction region. For
$m=0$,

$\displaystyle\delta_{F}=\frac{\lambda^{2}}{2\pi^{3}}\frac{s_{T}}{s_{X}(s_{X}^{2}+s_{T}^{2})}.$
(20)

##  III The perturbative Unruh-DeWitt detector

Next, we consider the case where the field-apparatus interaction term is much
weaker compared to the Hamiltonian $\hat{H}_{A}$ of the apparatus. In this
case, we can evaluate the detection probability using perturbation theory. The
Unruh-DeWitt detectors have been used extensively as measurement models in
this regime. The weak coupling approximation is physically relevant, as it
applies, for example to photodetection and to neutrino detection.

Let us denote by $|\epsilon\rangle$, the eigenstates of the Hamiltonian
$\hat{H}_{A}$, $\hat{H}_{A}|\epsilon\rangle=\epsilon|\epsilon\rangle$. We take
the initial state of the apparatus $|\Omega\rangle$ to coincide with the
ground state of $\hat{H}_{A}$, and we choose the energy scale so that
$\hat{H}_{A}|\Omega\rangle=0$. We assume an interaction Hamiltonian of the
form $f(t)\hat{\mu}\hat{\phi}(g)$, as in Sec. 2, where $\hat{\mu}$ is an
operator on the Hilbert space ${\cal H}_{A}$ of the apparatus. The switching
function $f(t)$ vanishes for $t$ outside $[0,T]$.

To leading order in the interaction, the excitation probability for the
apparatus is

$\displaystyle
P=\sum_{\epsilon>0}|\mu(\epsilon)|^{2}\int\mathrm{d}t\mathrm{d}^{3}x\int\mathrm{d}t^{\prime}\mathrm{d}^{3}x^{\prime}F(t,{\bm{x}})F(t,{\bm{x}}^{\prime})$
$\displaystyle\times
e^{-i\epsilon(t-t^{\prime})}\langle\psi_{0}|\hat{\phi}(t^{\prime},{\bm{x}}^{\prime})\hat{\phi}(t,{\bm{x}})|\psi_{0}\rangle,\;\;$
(21)

where $\mu(\epsilon)=\langle\epsilon|\hat{\mu}|\Omega\rangle$, and again
$F(t,{\bm{x}})=f(t)g({\bm{x}})$.

Consider a field state $|\psi_{0}\rangle$ with a definite number of particles.
Then, $P$ in Eq. (21) is a sum of two terms, $P=P_{0}+P_{1}$. The first term

$\displaystyle
P_{0}=\sum_{\epsilon>0}|\mu(\epsilon)|^{2}\int\frac{\mathrm{d}\bm{k}}{2\epsilon_{\bm{k}}}|\tilde{F}(\epsilon+\epsilon_{\bm{k}},{\bm{k}})|^{2}$
(22)

is independent of the initial state of the quantum field, and it is due to the
excitation of the vacuum from the switching of the interaction. In Eq. (22),
$\tilde{F}(\omega,{\bm{k}})$ is the Fourier transform of $F(t,{\bm{x}})$, and
$\mathrm{d}{\bm{k}}=\mathrm{d}^{3}k/(2\pi)^{3/2}$. $P_{0}$ functions as a
noise term in the measurement scheme; it becomes stronger for smearing
functions that are strongly localized in spacetime. If $\tau$ is the effective
duration of the interaction, the noise is suppressed if $\tau\epsilon_{1}>>1$,
where $\epsilon_{1}$ is the energy of the first excited state in the detector.

The second term $P_{1}$ depends on the one-particle reduced density matrix of
the field
$\rho_{0}({\bm{k}},{\bm{k}}^{\prime})=\langle\psi_{0}|\hat{a}^{\dagger}_{\bm{k}}\hat{a}_{\bm{k}}|\psi_{0}\rangle$.
For simplicity, we consider a pure
$\rho_{0}({\bm{k}},{\bm{k}}^{\prime})=\psi_{0}({\bm{k}})\psi_{0}^{*}({\bm{k}}^{\prime})$,
to obtain

$\displaystyle P_{1}$
$\displaystyle=\sum_{\epsilon>0}|\mu(\epsilon)|^{2}\left|\int\frac{\mathrm{d}\bm{k}}{\sqrt{2\epsilon_{\bm{k}}}}\tilde{F}^{*}(\epsilon_{\bm{k}}+\epsilon,\bm{k})\psi_{0}({\bm{k}})\right|^{2}$
(23)
$\displaystyle+\sum_{\epsilon>0}|\mu(\epsilon)|^{2}\left|\int\frac{\mathrm{d}\bm{k}}{\sqrt{2\epsilon_{\bm{k}}}}\tilde{F}^{*}(\epsilon_{\bm{k}}-\epsilon,\bm{k})\psi_{0}(\bm{k})\right|^{2}.$

is the ”signal” as it depends on the initial field state, e.g, a wavepacket of
the form

$\ket{\psi_{0}}=\int\mathrm{d}\bm{k}\psi_{0}(\bm{k})\hat{a}^{\dagger}(\bm{k})\ket{0}.$
(24)

A proper measurement requires a large signal-to-noise ratio, hence,
$\epsilon_{1}\tau>>1$. In this regime, the first term in the right hand side
of Eq. (23) is negligible.

This model functions as a particle detector for sufficiently large duration of
the interaction. To show this, we write the excitation probability (23) as a
spacetime integral:

$\displaystyle
P_{1}=\sum_{\epsilon>0}|\mu(\epsilon)|^{2}\left|\int\mathrm{d}t\mathrm{d}^{3}xe^{i\epsilon
t}F(t,\bm{x})\psi_{\textsc{nw}}(t,\bm{x})\right|^{2}$
$\displaystyle+\sum_{\epsilon>0}|\mu(\epsilon)|^{2}\left|\int\mathrm{d}t\mathrm{d}^{3}xe^{-i\epsilon
t}F(t,\bm{x})\psi_{\textsc{nw}}(t,\bm{x})\right|^{2}$

where we defined the Newton-Wigner wavefunction [39]

$\psi_{\textsc{NW}}(t,\bm{x})=\int\frac{\mathrm{d}\bm{k}}{\sqrt{2\epsilon_{\bm{k}}}}\psi_{0}(\bm{k})e^{i\bm{k}\bm{x}-i\epsilon_{\bm{k}}t}.$
(25)

Assume an almost monochromatic initial state with momentum ${\bm{p}}_{0}$,
such that $\psi_{\textsc{NW}}(0,{\bm{x}})$ is concentrated around a point
${\bm{x}}_{0}$. Then, $\psi_{NW}(t,{\bm{x}})$ is concentrated in a world tube
$W$ that surrounds surrounding the classical path
${\bm{x}}(t)={\bm{x}}_{0}+{\bm{p}}_{0}t$. Let $C$ stand for the region of
support of the function $F$. The probability density $P_{1}$ is appreciably
different from zero only if $W\cap C\neq\emptyset$. In a semi-classical
language, this means that the bundle of classical particle paths associated
with the initial conditions cross into the interaction region. This conforms
to our classical intuitions about the behavior of a particle detector.

Realistic field-detector couplings are typically isotropic, that is, the
smearing function $g({\bm{x}})$ is a function only of $|{\bm{x}}|$. In this
case, $\tilde{F}$ depends only on $k=|{\bm{k}}|$, and the excitation
probability $P_{1}$ involves an integral $\int d{\bf n}\psi(k,{\bf n})$, where
${\bf n}={\bf k}/|{\bf k}|$. Writing
$\psi({\bm{k}})=\sum_{\ell,m_{\ell}}\psi_{\ell,m_{\ell}}(k)Y_{\ell,m_{\ell}}({\bf
n})$, in terms of spherical harmonics $Y_{\ell,m_{\ell}}({\bf n})$, we find
that only the zero angular momentum $\psi_{0},0(k)$ component survives. This
means that propagation is essentially one-dimensional, along the axis that
connects the source and the detector, again in accordance with our classical
intuition. The fact that the detector ‘clicks’ in response to zero angular
momentum spherical waves (and not plane waves) can be used to explain the
counter-intuitive behavior that was reported in [40], that in certain
spacetime dimensions the detector’s excitation probability decreases in the
monochromatic limit.

In the limit of an interaction that is adiabatically switched on, that is, for
$f(t)=1$, we have $\tilde{F}(\epsilon,{\bf k})=2\pi\delta(\epsilon)g(k)$, and

$\displaystyle
P_{1}=\sum_{\epsilon>0}\frac{\epsilon|\mu(\epsilon)g(\sqrt{\epsilon^{2}-m^{2}})|^{2}}{8\pi^{2}(\epsilon^{2}-m^{2})}|\psi_{0,0}(\sqrt{\epsilon^{2}-m^{2}})|^{2}.$

If the energy levels of the detector are well separated, and the initial state
is approximately monochromatic, then the detection probability is
characterized by resonances with respect to the incoming energy. However,
these resonances are not of the standard form (Breit-Wigner for scattering, or
Wigner-Weisskopf for photodetection), but they mirror the profile of the
wavefunction $\psi_{0,0}$. This behavior is arguably unphysical, or at least,
it does not account for the observed resonances. However, this is to be
expected, because resonances cannot be consistently identified from the
leading-order terms in perturbation theory [49, 50]; the derivation of the
standard form requires at least a partial resummation of a perturbative
series—see, for example, chapter 7 in [51], or chapter 18 in [52].

##  IV Quantum Brownian motion as a QFT measurement model

In Secs. 2 and 3, we saw that a field-particle interaction with an UDW-type
coupling corresponds (i) to a measurement of field properties if the apparatus
Hamiltonian is negligible in compared to the interaction term, and (ii) to a
measurement of particle properties if the interaction is weak compared to the
apparatus Hamiltonian. In this section, we present a model that incorporates
both cases for different parameter regimes. Furthermore, this model does not
require the trick of switching on the field-apparatus interaction.

To this end, we consider a Hamiltonian
$\hat{H}=\hat{H}_{S}+\hat{H}_{A}+\hat{H}_{I}$, where $\hat{H}_{S}$ is the
field Hamiltonian (5), the apparatus Hamiltonian

$\displaystyle\hat{H}_{A}=\frac{1}{2M}\hat{P}^{2}+\frac{1}{2}M\omega_{0}^{2}\hat{X}^{2}$
(26)

describes an harmonic oscillator of mass $M$ and frequency $\omega_{0}$, and
the interaction Hamiltonian is $\hat{H}_{I}=\lambda\hat{\phi}(g)\hat{X}$.

The total Hamiltonian is of the Quantum Brownian motion (QBM) form

$\displaystyle\hat{H}=\frac{\hat{P}^{2}}{2M}+\frac{1}{2}M\omega_{0}^{2}\hat{X}^{2}$
$\displaystyle+\sum_{i}\left(\frac{\hat{p}_{i}^{2}}{2m_{i}}+\frac{1}{2}m_{i}\omega_{i}^{2}\hat{q}_{i}^{2}\right)$
$\displaystyle+\hat{X}\sum_{i}c_{i}\hat{q}_{i},$ (27)

i.e., it describes the interaction of a distinguished oscillator with a bath
of oscillators labeled by $i$.

QBM models admit exact solutions for their time evolution [48, 54, 55, 56,
47]. Here, we follow the presentation of [52]. We first define the dissipation
kernel

$\displaystyle\gamma(t):=\sum_{i}\frac{c_{i}^{2}}{2m_{i}\omega_{i}^{2}}\cos(\omega_{i}t),$
(28)

and then we identify the function $u(t)$ that solves the integrodifferential
equation

$\displaystyle\ddot{u}(t)+\bar{\omega}^{2}u(t)+\frac{2}{M}\int_{0}^{t}dt^{\prime}\gamma(t-t^{\prime})\dot{u}(t^{\prime})=0,$
(29)

subject to initial conditions $u(0)=0$ and $\dot{u}(0)=1$. We have defined the
regularized frequency $\bar{\omega}$ by
$\bar{\omega}^{2}=\omega_{0}^{2}-2\gamma(0)/M$.

Then, we can solve for $\hat{X}(t)$ and $\hat{P}(t)$ in the Heisenberg
picture,

$\displaystyle\hat{X}(t)$ $\displaystyle=$
$\displaystyle\dot{u}(t)\hat{X}(0)+\frac{1}{M}u(t)\hat{P}(0)$ (30)
$\displaystyle-$
$\displaystyle\frac{1}{M}\int_{0}^{t}dt^{\prime}u(t-t^{\prime})\sum_{i}c_{i}\hat{q}_{i}^{0}(t^{\prime}),$
$\displaystyle\hat{P}(t)$ $\displaystyle=$ $\displaystyle
M\frac{d}{dt}\hat{X}(t).$ (31)

where $\hat{q}_{i}^{0}(t)$ stands for the evolution of the bath harmonic
oscillators in the absence of interactions.

A QBM model is specified by the distribution of frequencies $\omega_{i}$ and
couplings $c_{i}$ among the environment oscillators. These determine the
dissipation kernel uniquely, and then, the function $u(t)$, by the solution of
Eq. (29). In the context of the field-detector interaction, the bath
oscillators are labeled by the momenta ${\bf k}$, $\omega_{i}$ corresponds to
$\epsilon_{\bf k}$, the oscillator masses $m_{i}=1$, and the coupling
constants $c_{i}$ correspond to $\lambda\tilde{g}({\bf k})$, where $\tilde{g}$
is the Fourier transform of the smearing function $g({\bf x})$. The
dissipation kernel is then,

$\displaystyle\gamma(t)=\frac{\lambda^{2}}{4\pi^{2}}\int_{0}^{\infty}\frac{dkk^{2}|\tilde{g}(k)^{2}|}{\epsilon_{k}^{2}}\cos(\epsilon_{k}t),$
(32)

while the Heisenberg-picture operators $\hat{X}(t)$ and $\hat{P}(t)$ are given
by

$\displaystyle\hat{X}(t)=$
$\displaystyle\dot{u}(t)\hat{X}(0)+\frac{1}{M}u(t)\hat{P}(0)$
$\displaystyle-\frac{\lambda}{M}\int_{0}^{t}dt^{\prime}u(t-t^{\prime})\hat{\phi}_{g}(t^{\prime}),$
(33) $\displaystyle\hat{P}(t)=$ $\displaystyle
M\ddot{u}(t)\hat{X}(0)+\dot{u}(t)\hat{P}(0)$
$\displaystyle-\lambda\int_{0}^{t}dt^{\prime}\dot{u}(t-t^{\prime})\hat{\phi}_{g}(t^{\prime}).$
(34)

Here, $\hat{\phi}_{g}(t)$ is the Heisenberg-evolution of the smeared field
$\hat{\phi}(g)$.

We assume that the apparatus oscillator is prepared in an initial state with
$\langle\hat{X}\rangle=0$, $\langle\hat{P}\rangle=0$ and $c_{XP}=0$. Then, we
obtain the following equations for the mean and the variances of the
detector’s position and momentum.

$\displaystyle\langle\hat{X}(t)\rangle=-\frac{\lambda}{M}\int_{0}^{t}dt^{\prime}u(t-t^{\prime})\langle\hat{\phi}_{g}(t^{\prime})\rangle,$
(35)
$\displaystyle\langle\hat{P}(t)\rangle=-\lambda\int_{0}^{t}dt^{\prime}\dot{u}(t-t^{\prime})\langle\hat{\phi}_{g}(t^{\prime})\rangle,$
(36)
$\displaystyle\langle\hat{X}^{2}(t)\rangle=\sigma_{X}^{2}\dot{u}^{2}(t)+\frac{\sigma_{P}^{2}}{M^{2}}{u}^{2}(t)$
(37)
$\displaystyle+\frac{\lambda^{2}}{M^{2}}\int_{0}^{t}dt^{\prime}\int_{0}^{t}dt^{\prime\prime}u(t-t^{\prime})u(t-t^{\prime\prime})\langle\hat{\phi}_{g}(t^{\prime})\hat{\phi}_{g}(t^{\prime\prime})\rangle,$
$\displaystyle\langle\hat{P}^{2}(t)\rangle=M^{2}\sigma_{X}^{2}\ddot{u}^{2}(t)+\sigma_{P}^{2}\dot{u}^{2}(t)$
(38)
$\displaystyle+\lambda^{2}\int_{0}^{t}dt^{\prime}\int_{0}^{t}dt^{\prime\prime}\dot{u}(t-t^{\prime})\dot{u}(t-t^{\prime\prime})\langle\hat{\phi}_{g}(t^{\prime})\hat{\phi}_{g}(t^{\prime\prime})\rangle.$

To proceed further, we need to evaluate the function $u(t)$. We will analyze
the case $m=0$. For concreteness, we take $|g(k)|^{2}=e^{-k/\Lambda}$, where
$\Lambda$ is an ultra-violet cut-off; $\Lambda^{-1}$ is the size of the
interaction region. Then,

$\displaystyle\gamma(t)=\frac{\lambda^{2}}{4\pi^{2}}\frac{\Lambda}{1+\Lambda^{2}t^{2}}.$
(39)

For $\Lambda>>\bar{\omega}$, we can approximate $\gamma(t)$ in Eq. with
$\frac{\lambda^{2}}{4\pi}\delta(t)$. Then, Eq. (29) becomes

$\displaystyle\ddot{u}+2\Gamma\dot{u}+\bar{\omega}^{2}u=0,$ (40)

where $\Gamma=\frac{\lambda^{2}}{16\pi M}$. This is the equation of a damped
harmonic oscillator, with solution

$\displaystyle u(t)=\left\\{\begin{array}[]{ccc}\omega^{-1}\sin(\omega
t)e^{-\Gamma
t},&\omega=\sqrt{\bar{\omega}^{2}-\Gamma^{2}},&\Gamma<\bar{\omega}\\\
w^{-1}\sinh(wt)e^{-\Gamma
t},&w=\sqrt{\Gamma^{2}-\bar{\omega}^{2}},&\Gamma>\bar{\omega}\end{array}\right.$
(43)

At the strong-coupling limit, $\bar{\omega}/\Gamma<<1$, and
$u(t)\simeq\frac{1}{2\Gamma}(e^{-\frac{\bar{\omega}^{2}}{2\Gamma}t}-e^{-2\Gamma
t})$, i.e., $u(t)$ interpolates between 0 at $t=0$ and $\frac{1}{2\Gamma}$ at
times $t$, such that $\Gamma^{-1}<<t<<\Gamma\bar{\omega}^{-2}$. Therefore, we
can approximate $u(t)\simeq\frac{1}{2\Gamma}\theta(t)$. Then, Eq. (36)
becomes,

$\displaystyle\langle\hat{P}(t)\rangle=-\frac{\lambda}{4\Gamma}\langle\hat{\phi}_{g}(t)\rangle,$
(44)

that is, the detector’s momentum observable records the field’s expectation
value. Eq. (44 holds for any field pulse that (i) arrives at the detector at
times $t>>\Gamma^{-1}$, and (ii) is of duration $\tau$, such that
$\bar{\omega}^{2}\tau/\Gamma<<1$. In this regime, $(\Delta
P)^{2}(t)\simeq\frac{\lambda^{2}}{16\Gamma^{2}}(\Delta\phi_{g})^{2}(t)$.
Hence, the oscillator acts as an antenna: it records faithfully ”classical”
field states, that is, field states with small quantum fluctuations.

For an initial state with a fixed number of particles,
$\langle\hat{X}(t)\rangle=0$ and $\langle\hat{P}(t)\rangle=0$. Then, the
natural quantity for a detection signal is the energy of the harmonic
oscillator,

$\displaystyle\langle\hat{H}(t)\rangle=\frac{(\Delta\hat{P})^{2}(t)}{2M}+\frac{M\bar{\omega}^{2}(\Delta\hat{X})^{2}(t)}{2}.$
(45)

The expressions simplify if we assume an initial state that is localized at
${\bf x}={\bf L}$, such that $|{\bf L}|\Gamma>>1$. Then, the signal from
incoming particles appears at times $t>|{\bf L}|$, when many transient noise
terms have vanished. We obtain (for $m=0$),

$\displaystyle\langle\hat{H}(t)\rangle=h(t)+{\cal N},$ (46)

where

$\displaystyle{\cal
N}=\frac{2\Gamma}{\pi}\int_{0}^{\infty}\frac{dkk|g(k)|^{2}(\bar{\omega}^{2}+k^{2})}{(\bar{\omega}^{2}-k^{2})^{2}+4\Gamma^{2}k^{2}}$
(47)

is the vacuum noise, and

$\displaystyle h$ $\displaystyle(t)$ (48)
$\displaystyle=8\pi\Gamma\left|\int\frac{d{\bf k}}{\sqrt{2\epsilon_{\bf
k}}}\psi_{0}({\bf k})g({\bf k})\frac{\bar{\omega}^{2}+\epsilon_{\bf
k}^{2}}{\bar{\omega}^{2}+(\Gamma-i\epsilon_{\bf k})^{2}}e^{-i\epsilon_{\bf
k}t}\right|^{2}$

is the detection signal111 The noise term diverges logarithmically as the size
of the detector vanishes, that is, $\tilde{g}\to 1$. Note that we have used
the pointlike renormalized dynamics for solving the integrodifferential
equation, but we kept the smearing in the source term in order to calculate
the noise. Although this might seem artificial, it can be shown that this
noise term represents the leading order behavior of the noise in a power
expansion involving a UV cutoff.. We assume that the initial single-particle
state is almost monochromatic at ${\bf k}={\bf k}_{0}$. Then,

$\displaystyle
h(t)=\frac{8\pi\Gamma|g(\epsilon_{0})|^{2}(\epsilon_{0}^{2}+\bar{\omega}^{2})}{(\bar{\omega}-\epsilon_{0}^{2})^{2}+4\Gamma^{2}\epsilon_{0}^{2}}|\psi_{NW}(t,0)|^{2},$
(49)

where $\epsilon_{0}=|{\bf k}_{0}|$.

We see that the detection signal appears after the time $t=|{\bf L}|$ that the
Newton-Wigner wave packet arrives at ${\bf x}=0$, as in the perturbative
Unruh-DeWitt model of Sec. 3. Note that no switching function has been
employed. Except for noise terms, the field interacts with the detector only
when the wavepacket has reached the locus of the detector.

We also note that for $\Gamma<<\bar{\omega}$, the signal is characterized by a
resonance that is independent of the shape of the initial state. The resonance
is described by the relativistic Breit-Wigner distribution [53].

We conclude that the harmonic oscillator in this model functions as a particle
detector that records energy at small $\Gamma$, and as an antenna that records
the field’s amplitude at large $\Gamma$. The model is well defined also in the
regime of intermediate $\Gamma$, thus interpolating between the field regime
(also described by von Neumann’s model) and the particle regime (also
described by the perturbative Unruh-DeWitt detector).

##  V Conclusions

In this article, we explored the principle of field-particle duality using
simple models for QFT measurements. We saw that von Neumann’s theory describes
the measurement of classical field amplitude in the regime of strong field-
detector coupling. We also saw that the perturbative Unruh-DeWitt model indeed
describes particle detection in the regime of weak field-detector coupling,
with the caveat that it does not provide a description of resonances
reproducing observed experimental results. Indeed, resonances of the Breit-
Wigner type require absorption effects that can only be captured
nonperturbatively, or after a partial resummation of the perturbative series.

Then, we presented a new measurement model that is based on QBM. This model is
exactly solvable. It enables us to show that in the long-time limit, field
amplitude measurements occur at strong coupling and particle measurements at
weak coupling. It also accounts for the regime of intermediate couplings that
corresponds to the transition from field to particle behavior. The fact that
this transition is modulated by a single parameter is a novel aspect of field-
particle duality, which we expect to be experimentally accessible in set-ups
that allow for tunable coupling of a quantum oscillator with a quantum field.

Besides experimental relevance, the model presented here is a substantial
improvement over existing models for QFT measurements that are based on the
coupling of fields to pointlike systems. Being exactly solvable, it provides a
consistent analysis of the strong coupling regime, and it also provides a
consistent description of resonances. Furthermore, the model remains exactly
solvable if we introduce multiple detectors [55, 56], for example, in order to
analyze causality in relation to field measurements, relativistic quantum
communication [57], or even measurements in an extended spatial region.

Acknowledgements: MP acknowledges support of the ID$\\#$ 62312 grant from the
John Templeton Foundation, as part of the ‘The Quantum Information Structure
of Spacetime’ Project (QISS). JRR’s work was supported by the Regional
Government of Castilla y León (Junta de Castilla y León) and by the Ministry
of Science and Innovation MICIN and the European Union NextGenerationEU (PRTR
C17.I1). JRR would like to thank the Univerisity of Patras and Charis
Anastopoulos’ group for the hospitality during the course of this research.

## References

  * [1] A. S. Blum, The state is not abolished, it withers away: How quantum field theory became a theory of scattering, Studies in History and Philosophy of Modern Physics 60:46-80 (2017).
  * [2] M. Papageorgiou and D. Fraser, Eliminating the ‘impossible’: Recent progress on local measurement theory for quantum field theory, http://philsci-archive.pitt.edu/22322/ (2023).
  * [3] C. Anastopoulos and N. Savvidou, Quantum Information in Relativity: The Challenge of QFT Measurements, Entropy 24, 4 (2022).
  * [4] C. Anastopoulos, B.L. Hu and K. Savvidou, Quantum field theory based quantum information: Measurements and correlations, Ann. Phys. 450, 169239 (2023).
  * [5] J. Earman and G. Valente, Relativistic Causality in Algebraic Quantum Field Theory, International Studies in the Philosophy of Science 28 (1):1-48.
  * [6] P. Busch, P. J. Lahti and P. Mittelstaedt, The Quantum Theory of Measurement, Lecture Notes in Physics Monographs, volume 2 (1996).
  * [7] R. F. Streater and A. S. Wightman, PCT, Spin and statistics, and all that (Princeton University Press, 2000).
  * [8] D. Malament In Defense of Dogma: Why There Cannot Be a Relativistic Quantum Mechanical Theory of (Localizable) Particles, In R. Clifton (ed.), Perspectives on Quantum Reality. Kluwer Academic Publishers (1996).
  * [9] D. Fraser, The fate of ’particles’ in quantum field theories with interactions, Studies in History and Philosophy of Science Part B: Studies in History and Philosophy of Modern Physics 39 (4):841-859 (2008).
  * [10] H. Halvorson and R. Clifton, No Place for Particles in Relativistic Quantum Theories?, Philosophy of Science, Vol. 69, No. 1 (2002).
  * [11] R. Haag, Local Quantum Physics: Fields, Particles, Algebras (Springer, 1996).
  * [12] M. Kuhlmann, Quantum Field Theory, The Stanford Encyclopedia of Philosophy, https://plato.stanford.edu/archives/sum2023/entries/quantum-field-theory/ (2023).
  * [13] R. J. Glauber, The Quantum Theory of Optical Coherence, Phys. Rev. 130, 2529 (1963).
  * [14] D. Rideout et al, Fundamental Quantum Optics Experiments Conceivable with Satellites – Reaching Relativistic Distances and Velocities, Class. Quantum Grav. 29, 224011 (2012).
  * [15] M. Mohageg et al, The Deep Space Quantum Link: Prospective Fundamental Physics Experiments using Long-Baseline Quantum Optics, EJP Quantum Technology, 9, 25 (2022).
  * [16] C. Anastopoulos, B.L. Hu and K. Savvidou, Towards a Field-Theory based Relativistic Quantum Information J. Phys.: Conf. Ser. 2533, 012004 (2023).
  * [17] J. von Neumann, Mathematical Foundations of Quantum Mechanics (Princeton University Press, Princeton, 1955).
  * [18] K. E. Hellwig and K. Kraus, Formal description of measurements in local quantum field theory, Phys. Rev. D 1, 566 (1970).
  * [19] K. Okamura and M. Ozawa, Measurement theory in local quantum physics, J. Math. Phys. 57, 015209 (2015).
  * [20] C. J. Fewster and R. Verch, Quantum fields and local measurements, Comm. Math. Phys. 378, 851 (2020).
  * [21] C. J. Fewster, I. Jubb, and M. H. Ruep, Asymptotic measurement schemes for every observable of a quantum field theory, Ann. Henri Poincaré 24, 1137–1184 (2023)
  * [22] C. Anastopoulos and N. Savvidou, Time-of-Arrival Probabilities for General Particle Detectors, Phys. Rev. A86, 012111 (2012).
  * [23] C. Anastopoulos and N. Savvidou, Time-of-Arrival Correlations, Phys. Rev. A95, 032105 (2017).
  * [24] C. Anastopoulos and N. Savvidou, Time of arrival and Localization of Relativistic Particles, J. Math. Phys. 60, 0323301 (2019).
  * [25] R. B. Griffiths, Consistent Quantum Theory (Cambridge University Press, Cambridge 2003).
  * [26] R. Omnés, The Interpretation of Quantum Mechanics, (Princeton University Press, Princeton 1994).
  * [27] R. Omnés, Understanding Quantum Mechanics (Princeton University Press, Princeton 1999).
  * [28] M. Gell-Mann and J. B. Hartle, Classical Equations for Quantum Systems, Phys. Rev. D47, 3345 (1993).
  * [29] J.B. Hartle, Spacetime Quantum Mechanics and the Quantum Mechanics of Spacetime in ‘Gravitation and Quantizations’, in the Proceedings of the 1992 Les Houches Summer School, ed. by B. Julia and J. Zinn- Justin, Les Houches Summer School Proceedings, Vol. LVII, (North Holland, Amsterdam, 1995); [gr-qc/9304006].
  * [30] J. Polo-Gómez, L. J. Garay, L. J. and E. Martín-Martínez, A detector-based measurement theory for quantum field theory, Phys. Rev. D 105, 065003 (2022).
  * [31] T. R. Perche, Localized nonrelativistic quantum systems in curved spacetimes: A general characterization of particle detector models, Phys. Rev. D106, 025018 (2022).
  * [32] E. Martin-Martinez, T. R. Perche and B. de S. L. Torres, Broken covariance of particle detector models in relativistic quantum information, Phys. Rev. D 103, 025007 (2021).
  * [33] E. Martin-Martinez, Causality issues of particle detector models in QFT and quantum optics, Phys. Rev. D 92, 104019 (2015).
  * [34] J. de Ramon, M. Papageorgiou and E. Martin-Martinez, Relativistic causality in particle detector models: Faster-than-light signaling and impossible measurements , Phys. Rev. D 103, 085002 (2021).
  * [35] J. de Ramon, M. Papageorgiou and E. Martin-Martinez, Causality and signalling in noncompact detector-field interactions, Phys. Rev. D 108, 045015 (2023).
  * [36] B. L. Hu, S-Y Lin, J. Louko, Relativistic Quantum Information in Detectors–Field Interactions, Class. Quantum Grav. 29, 224005 (2012).
  * [37] W. G. Unruh, Notes on Black Hole Evaporation, Phys. Rev. D14, 870 (1976).
  * [38] B. S. DeWitt, Quantum Gravity: the New Synthesis in General Relativity: An Einstein Centenary Survey, ed. by S. W. Hawking and W. Israel (Cambridge University Press, Cambridge, 1979), p. 680\.
  * [39] T. D. Newton and E. P. Wigner Localized states for elementary systems, Reviews of Modern Physics, Volume 21, Number 3 (1949).
  * [40] E. Tjoa, I. Lopez- Gutierrez, A. Sachs, E. Martin-Martinez, What makes a particle detector click, Phys. Rev. D 103, 125021 (2021).
  * [41] A. Blais, A. L. Grimsmo, S. M. Girvin, and A. Wallraff, Circuit quantum electrodynamics, Rev. Mod. Phys. 93, 025005 (2021).
  * [42] M. Aspelmeyer, T. J. Kippenberg, and F. Marquardt (edts.) Cavity Optomechanics (Springer: Berlin, Heidelberg, 2014).
  * [43] C. Barceló, S Liberati, and Matt Visser, Analogue gravity from Bose-Einstein condensates, Class. Quantum Grav. 18, 1137 (2001).
  * [44] C. Viermann, M. Sparn, N. Liebster, et al, Quantum field simulator for dynamics in curved spacetime Nature 611, 260 (2022).
  * [45] C. Anastopoulos and N. Savvidou, Quantum probabilities for time-extended alternatives, J. Math. Phys. 48, 032106 (2008).
  * [46] A. O. Caldeira, and A. J. Leggett, Path Integral Approach to Quantum Brownian Motion, Physica A121, 587 (1983).
  * [47] W. G. Unruh, and W. H. Zurek, Reduction of a wave packet in quantum Brownian motion Phys. Rev. D 40, 1071 (1989).
  * [48] B. L. Hu, J. P. Paz, and Y. Zhang, Quantum Brownian motion in a general environment: Exact master equation with nonlocal dissipation and colored noise, Phys. Rev. D46, 2843 (1992).
  * [49] L. S. Rodberg, and R. M. Thaler, Introduction to the Quantum Theory of Scattering (New York: Academic Press, 1967).
  * [50] C. Cohen-Tannoudji, B. Diu, and F. Laloë, Quantum Mechanics, vol. 2 (New York: Wiley, 2019).
  * [51] M. Peskin and D. Schroeder,An Introduction to Quantum Field Theory (Reading, MA: Addison Wesley, 1995).
  * [52] C. Anastopoulos,Quantum Theory: A Foundational Approach (Cambridge: Cambridge University Press, 2023).
  * [53] G. Breit, and E. Wigner, Capture of slow neutrons, Phys. Rev.Vol. 49, Iss. 7 — April 1936.
  * [54] J. J. Halliwell, and T. Yu, Alternative derivation of the Hu-Paz-Zhang master equation of quantum Brownian motion Phys. Rev. D53, 2012 (1995)
  * [55] C. Anastopoulos, S. Kechribaris and D. Mylonas, Generalized uncertainty relations and entanglement dynamics in quantum Brownian motion models, Phys. Rev. A82, 042119 (2010).
  * [56] C. H. Fleming, A. Roura and B. L. Hu, Exact analytical solutions to the master equation of quantum Brownian motion for a general environment, Ann. Phys. 326, 1207 (2010).
  * [57] A. Lapponi, D. Moustos, D. E. Bruschi, and S. Mancini, Relativistic quantum communication between harmonic oscillator detectors, Phys. Rev. D 107, 125010 (2023).


On an Invariance Problem for Parameterized Concurrent Systems

Marius Bozga Lucas Bueri Radu Iosif1
Univ. Grenoble Alpes, CNRS, Grenoble INP, VERIMAG, 38000, France

[1]
E. Ahrens, M. Bozga, R. Iosif, and J. Katoen.
Local reasoning about parameterized reconfigurable distributed
CoRR, abs/2107.05253, 2021.

[2]
Y. Bar-Hillel, M. Perles, and E. Shamir.
On formal properties of simple phrase structure grammars.
Sprachtypologie und Universalienforschung, 14:143–172, 1961.

[3]
R. Bloem, S. Jacobs, A. Khalimov, I. Konnov, S. Rubin, H. Veith, and J. Widder.
Decidability of Parameterized Verification.
Synthesis Lectures on Distributed Computing Theory. Morgan &
Claypool Publishers, 2015.

[4]
M. Bozga, L. Bueri, and R. Iosif.
Decision problems in a logic for reasoning about reconfigurable
distributed systems.
arXiv 2202.09637, cs.LO, 2022.

[5]
M. Bozga, J. Esparza, R. Iosif, J. Sifakis, and C. Welzel.
Structural invariants for the verification of systems with
parameterized architectures.
In Tools and Algorithms for the Construction and Analysis of
Systems - 26th International Conference, TACAS 2020, volume 12078 of LNCS, pages 228–246. Springer, 2020.

[6]
M. Bozga and R. Iosif.
Specification and safety verification of parametric hierarchical
distributed systems.
In Formal Aspects of Component Software - 17th International
Conference, FACS 2021, Virtual Event, October 28-29, 2021, Proceedings,
volume 13077 of Lecture Notes in Computer Science, pages 95–114.
Springer, 2021.

[7]
A. Bucchiarone and J. P. Galeotti.
Dynamic software architectures verification using dynalloy.
Electron. Commun. Eur. Assoc. Softw. Sci. Technol., 10, 2008.

[8]
C. Calcagno, P. W. O'Hearn, and H. Yang.
Local action and abstract separation logic.
In 22nd IEEE Symposium on Logic in Computer Science (LICS
2007), 10-12 July 2007, Wroclaw, Poland, Proceedings, pages 366–378. IEEE
Computer Society, 2007.

[9]
B. Cook, C. Haase, J. Ouaknine, M. J. Parkinson, and J. Worrell.
Tractable reasoning in a fragment of separation logic.
In CONCUR, volume 6901 of Lecture Notes in Computer
Science, pages 235–249. Springer, 2011.

[10]
J. Dormoy, O. Kouchnarenko, and A. Lanoix.
Using temporal logic for dynamic reconfigurations of components.
In L. S. Barbosa and M. Lumpe, editors, Formal Aspects of
Component Software - 7th International Workshop, FACS 2010, volume 6921 of
Lecture Notes in Computer Science, pages 200–217. Springer, 2010.

[11]
A. El-Hokayem, M. Bozga, and J. Sifakis.
A temporal configuration logic for dynamic reconfigurable systems.
In C. Hung, J. Hong, A. Bechini, and E. Song, editors, SAC
'21: The 36th ACM/SIGAPP Symposium on Applied Computing, Virtual Event,
Republic of Korea, March 22-26, 2021, pages 1419–1428. ACM, 2021.

[12]
K. Foerster and S. Schmid.
Survey of reconfigurable data center networks: Enablers, algorithms,
SIGACT News, 50(2):62–79, 2019.

[13]
D. Hirsch, P. Inverardi, and U. Montanari.
Graph grammars and constraint solving for software architecture
In Proceedings of the Third International Workshop on Software
Architecture, ISAW '98, page 69–72, New York, NY, USA, 1998. Association
for Computing Machinery.

[14]
Y. Kesten, A. Pnueli, E. Shahar, and L. D. Zuck.
Network invariants in action.
In CONCUR 2002 - Concurrency Theory, 13th International
Conference, volume 2421 of LNCS, pages 101–115. Springer, 2002.

[15]
A. Lanoix, J. Dormoy, and O. Kouchnarenko.
Combining proof and model-checking to validate reconfigurable
Electron. Notes Theor. Comput. Sci., 279(2):43–57, 2011.

[16]
D. Le Metayer.
Describing software architecture styles using graph grammars.
IEEE Transactions on Software Engineering, 24(7):521–533,

[17]
D. Lesens, N. Halbwachs, and P. Raymond.
Automatic verification of parameterized linear networks of processes.
In The 24th ACM SIGPLAN-SIGACT Symposium on Principles of
Programming Languages, pages 346–357. ACM Press, 1997.

[18]
M. Noormohammadpour and C. S. Raghavendra.
Datacenter traffic control: Understanding techniques and tradeoffs.
IEEE Commun. Surv. Tutorials, 20(2):1492–1525, 2018.

[19]
J. C. Reynolds.
Separation logic: A logic for shared mutable data structures.
In 17th IEEE Symposium on Logic in Computer Science (LICS
2002), 22-25 July 2002, Copenhagen, Denmark, Proceedings, pages 55–74.
IEEE Computer Society, 2002.

[20]
Z. Shtadler and O. Grumberg.
Network grammars, communication behaviors and automatic verification.
In J. Sifakis, editor, Automatic Verification Methods for Finite
State Systems, International Workshop, volume 407 of LNCS, pages
151–165. Springer, 1989.

[21]
P. Wolper and V. Lovinfosse.
Verifying properties of large sets of processes with network
In Automatic Verification Methods for Finite State Systems,
International Workshop, volume 407 of LNCS, pages 68–80. Springer,