# Bisimulation Relations Between Automata,
Stochastic Differential Equations and Petri Nets

Mariken H.C. Everdij Henk A.P. Blom National Aerospace Laboratory NLR
Amsterdam, Netherlands<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Two formal stochastic models are said to be bisimilar if their solutions as a
stochastic process are probabilistically equivalent. Bisimilarity between two
stochastic model formalisms means that the strengths of one stochastic model
formalism can be used by the other stochastic model formalism. The aim of this
paper is to explain bisimilarity relations between stochastic hybrid automata,
stochastic differential equations on hybrid space and stochastic hybrid Petri
nets. These bisimilarity relations make it possible to combine the formal
verification power of automata with the analysis power of stochastic
differential equations and the compositional specification power of Petri
nets. The relations and their combined strengths are illustrated for an air
traffic example.

##  1 Introduction

Two formal stochastic models are said to be _bisimilar_ if their solutions as
a stochastic process (i.e. their _executions_) are probabilistically
equivalent [9, 24]. Bisimilarity relations between formal stochastic models
are very useful to study since they allow one stochastic model to take
advantage of the strengths of the other stochastic model. The aim of this
paper is to show bisimulation relations between three different stochastic
modelling formalisms: stochastic hybrid automata, stochastic differential
equations on hybrid space, and stochastic hybrid Petri nets. These
bisimulation relations make it possible to combine the formal verification
power of automata with the analysis power of stochastic differential equations
and the compositional specification power of Petri nets.

For the stochastic automata formalism, we take the _general stochastic hybrid
system_ (GSHS) theoretical setting developed by [8]. A GSHS is a hybrid
automaton defined on a hybrid state space. This hybrid state space consists of
a countable set of discrete modes, and per discrete mode a Euclidean subset.
Per discrete mode, a stochastic differential equation (SDE) is defined. Two
additional GSHS elements are a jump rate function and a GSHS transition
measure. The execution of these elements provides a stochastic process that
follows the solution of the SDE connected to the initial discrete mode. After
a time period, defined by the jump rate function, the process state may
spontaneously jump to another mode, defined by the GSHS transition measure. A
jump may also be forced if the process state hits the boundary of the
Euclidean subset. The GSHS execution is referred to as _general stochastic
hybrid process_ (GSHP). One of the main strengths of the automata formalism is
the availability of formal verification tools.

For the hybrid stochastic differential equations formalism we take the _hybrid
stochastic differential equations_ (HSDE) theoretical setting developed in a
series of complementary studies [2, 3, 21, 22]. A HSDE consists of a sequence
of SDEs on a hybrid state space, driven by a Poisson random measure. When the
Poisson random measure generates a multivariate point, a spontaneous jump
occurs. A jump may also be forced if the process state hits the boundary of a
Euclidean subset. The HSDE solution process is referred to as _general
stochastic hybrid process_ (GSHP). In [17] it is shown that whereas the GSHS
formalism is at some points more general than HSDE (for GSHS the dimension of
the Euclidean subset may depend on the discrete mode; for HSDE this dimension
is fixed), HSDE has the advantage of an established semi-martingale property
and includes the coverage of jump-linear systems.

For the stochastic hybrid Petri nets formalism, we take _stochastically and
dynamically coloured Petri nets_ (SDCPN) developed in a series of studies by
[14, 15, 16]. A Petri net has _places_ (circles), which model possible
discrete states or conditions, and which may contain one or more _tokens_
(dots), modelling which of these states are current. The places are connected
by _transitions_ (squares), which model state switches by removing input
tokens and producing output tokens along _arcs_ (arrows). In SDCPN, the tokens
have Euclidean-valued colours that follow SDEs. Some of the transitions remove
and produce tokens spontaneously, other transitions are forced and occur when
the colours of their input tokens reach the boundary of a Euclidean subset.
The collection of token colours in all places forms a _general stochastic
hybrid process_ (GSHP). The specific strength of SDCPN is their compositional
specification power, which makes available a hierarchical modelling approach
that separates local modelling issues from global modelling issues. This is
illustrated for a large distributed example in air traffic management [18],
which covers many distributed agents each of which interacts in a dynamic way
with the others. Other typical Petri net features are concurrency and
synchronisation mechanism, hierarchical and modular construction, and natural
expression of causal dependencies, in combination with graphical and
equational representation.

The aim of this paper is to illustrate the relations between SDCPN, GSHP, HSDE
and GSHS which show that SDCPN, GSHS and HSDE are bisimilar. This means that
if we take the elements of any one of these formalisms, we can construct the
elements of another formalism in such a way that their associated GSHPs are
probabilistically equivalent. Fig. 1 shows the relations between the
formalisms, and the key tools available for each of them.

Figure 1: Relationship between SDCPN, GSHS, GSHP and HSDE, and their key
properties and advantages. The [B] arrow is established in [2]. The [BL] arrow
is established in [8]. The [E1] arrows are established in [16]. The [E2]
arrows are established in [17].

With these relations, the properties and advantages of the various approaches
come within reach of each other. The compositional specification power of
SDCPN makes it relatively easy to develop a model for a complex system with
multiple interactions. Subsequently, in the analysis stage three alternative
approaches can be taken. The first is direct execution of SDCPN and evaluation
through e.g. Monte Carlo simulation. The second is mapping the SDCPN into a
GSHS and evaluating its execution, with the advantages of connection to formal
methods in automata theory and to optimal control theory [7]. The third is
mapping the SDCPN into HSDE and evaluating its solution, with the advantages
of stochastic analysis for semi-martingales [11, 12]. With the GSHP resulting
from any of these three means, properties become available such as convergence
of discretisation, existence of limits, existence of event probabilities,
strong Markov properties, and reachability analysis [8, 10, 13].

The organisation of this paper is as follows. Section 2 defines SDCPN and the
related SDCPN process. Section 3 presents an example SDCPN model for a simple
but illustrative air traffic situation. Section 4 defines GSHS and illustrates
how the example SDCPN can be mapped to a bisimilar GSHS. Section 5 defines
HSDE and illustrates how the example SDCPN can be mapped to a bisimilar HSDE.
Section 6 gives conclusions.

##  2 SDCPN

This section outlines _stochastically and dynamically coloured Petri net_
(SDCPN). For a more formal definition, we refer to [17].

#### Definition 2.1

(Stochastically and dynamically coloured Petri net.) _An SDCPN is a collection
of elements $({\cal P}$, ${\cal T}$, ${\cal A}$, ${\cal N}$, ${\cal S}$,
${\cal C}$, ${\cal I}$, ${\cal V}$, ${\cal W}$, ${\cal G}$, ${\cal D}$, ${\cal
F})$, together with an SDCPN execution prescription which makes use of a
sequence $\\{U_{i};i=0,1,\ldots\\}$ of independent uniform $U[0,1]$ random
variables, of independent sequences of mutually independent standard Brownian
motions $\\{B_{t}^{i,P};i=1,2,\ldots\\}$ of appropriate dimensions, one
sequence for each place $P$, and of five rules R0–R4 that solve enabling
conflicts._

###  2.1 SDCPN elements

The SDCPN elements (${\cal P}$, ${\cal T}$, ${\cal A}$, ${\cal N}$, ${\cal
S}$, ${\cal C}$, ${\cal I}$, ${\cal V}$, ${\cal W}$, ${\cal G}$, ${\cal D}$,
${\cal F}$) are defined as follows:

  * •

${\cal P}$ is a finite set of places.

  * •

${\cal T}$ is a finite set of transitions which consists of 1) a set ${\cal
T}_{G}$ of guard transitions, 2) a set ${\cal T}_{D}$ of delay transitions and
3) a set ${\cal T}_{I}$ of immediate transitions.

  * •

${\cal A}$ is a finite set of arcs which consists of 1) a set ${\cal A}_{O}$
of ordinary arcs, 2) a set ${\cal A}_{E}$ of enabling arcs and 3) a set ${\cal
A}_{I}$ of inhibitor arcs.

  * •

${\cal N}:{\cal A}\to{\cal P}\times{\cal T}\cup{\cal T}\times{\cal P}$ is a
node function which maps each arc $A\in{\cal A}$ to a pair of ordered nodes
${\cal N}(A)$, where a node is a place or a transition.

  * •

${\cal S}\subset\\{\mathbb{R}^{0},\mathbb{R}^{1},\mathbb{R}^{2},\ldots\\}$ is
a finite set of colour types, with $\mathbb{R}^{0}\triangleq\emptyset$.

  * •

${\cal C}:{\cal P}\to{\cal S}$ is a colour type function which maps each place
$P\in{\cal P}$ to a specific colour type. Each token in $P$ is to have a
colour in ${\cal C}(P)$. If ${\cal C}(P)=\mathbb{R}^{0}$ then a token in $P$
has no colour.

  * •

${\cal I}$ is a probability measure, which defines the initial marking of the
net: for each place it defines a number $\geq 0$ of tokens initially in it and
it defines their initial colours.

  * •

${\cal V}=\\{{\cal V}_{P};P\in{\cal P},{\cal C}(P)\neq\mathbb{R}^{0}\\}$ is a
set of token colour functions. For each place $P\in{\cal P}$ for which ${\cal
C}(P)\neq\mathbb{R}^{0}$, it contains a function ${\cal V}_{P}:{\cal
C}(P)\to{\cal C}(P)$ that defines the drift coefficient of a differential
equation for the colour of a token in place $P$.

  * •

${\cal W}=\\{{\cal W}_{P};P\in{\cal P},{\cal C}(P)\neq\mathbb{R}^{0}\\}$ is a
set of token colour matrix functions. For each place $P\in{\cal P}$ for which
${\cal C}(P)\neq\mathbb{R}^{0}$, it contains a measurable mapping ${\cal
W}_{P}:{\cal C}(P)\to\mathbb{R}^{n(P)\times h(P)}$ that defines the diffusion
coefficient of a stochastic differential equation for the colour of a token in
place $P$, where $h:{\cal P}\rightarrow\mathbb{N}$ and $n:{\cal
P}\rightarrow\mathbb{N}$ is such that ${\cal C}(P)=\mathbb{R}^{n(P)}$. It is
assumed that ${\cal W}_{P}$ and ${\cal V}_{P}$ satisfy conditions that ensure
a probabilistically unique solution of each stochastic differential equation.

  * •

${\cal G}=\\{{\cal G}_{T};T\in{\cal T}_{G}\\}$ is a set of transition guards.
For each $T\in{\cal T}_{G}$, it contains a transition guard ${\cal G}_{T}$,
which is an open Euclidean subset with boundary $\partial{\cal G}_{T}$.

  * •

${\cal D}=\\{{\cal D}_{T};T\in{\cal T}_{D}\\}$ is a set of transition delay
rates. For each $T\in{\cal T}_{D}$, it contains a locally integrable
transition delay rate ${\cal D}_{T}$.

  * •

${\cal F}=\\{{\cal F}_{T};T\in{\cal T}\\}$ is a set of firing measures. For
each $T\in{\cal T}$, it contains a firing measure ${\cal F}_{T}$, which
generates the number and colours of the tokens produced when transition $T$
fires, given the value of the vector that collects all input tokens: For each
output arc, zero or one token is produced. For each fixed $H$, ${\cal
F}_{T}(H;\cdot)$ is measurable. For any $c$, ${\cal F}_{T}(\cdot;c)$ is a
probability measure.

For the places, transitions and arcs, the graphical notation is as in Figure
2.

Figure 2: Graphical notation for places, transitions and arcs in an SDCPN

###  2.2 SDCPN execution

The execution of an SDCPN provides a series of increasing stopping times,
$\\{\tau_{i};i=0,1,\ldots\\}$, $\tau_{0}=0$, with for
$t\in(\tau_{k},\tau_{k+1})$ a fixed number of tokens per place and per token a
colour which is the solution of a stochastic differential equation.

#### Initiation.

The probability measure ${\cal I}$ characterises an initial marking at
$\tau_{0}$, i.e. it gives each place $P\in{\cal P}$ zero or more tokens and
gives each token in $P$ a colour in ${\cal C}(P)$, i.e. a Euclidean-valued
vector.

#### Token colour evolution.

For each token in each place $P$ for which ${\cal C}(P)\neq\mathbb{R}^{0}$: if
the colour of this token is equal to $C_{0}^{P}$ at time $t=\tau_{0}$, and if
this token is still in this place at time $t>\tau_{0}$, then the colour
$C_{t}^{P}$ of this token equals the probabilistically unique solution of the
stochastic differential equation $dC_{t}^{P}={\cal V}_{P}(C_{t}^{P})dt+{\cal
W}_{P}(C_{t}^{P})dB_{t}^{i,P}$ with initial condition
$C_{\tau_{0}}^{P}=C_{0}^{P}$, and with $\\{B_{t}^{i,P}\\}$ an
$h(P)$-dimensional standard Brownian motion. Each token in a place for which
${\cal C}(P)=\mathbb{R}^{0}$ remains without colour.

#### Transition enabling.

A transition $T$ is pre-enabled if it has at least one token per incoming
ordinary and enabling arc in each of its input places and has no token in
places to which it is connected by an inhibitor arc. For each transition $T$
that is pre-enabled at $\tau_{0}$, consider one token per ordinary and
enabling arc in its input places and write $C_{t}^{T}$, $t\geq\tau_{0}$, as
the column vector containing the colours of these tokens; $C_{t}^{T}$ evolves
through time according to its corresponding token colour functions. If this
vector is not unique (i.e., if one input place contains several tokens per
arc), all possible such vectors are executed in parallel. A transition $T$ is
enabled if it is pre-enabled and a second requirement holds true. For
$T\in{\cal T}_{I}$, the second requirement automatically holds true at the
time of pre-enabling. For $T\in{\cal T}_{G}$, the second requirement holds
true when $C_{t}^{T}\in\partial{\cal G}_{T}$. For $T\in{\cal T}_{D}$, the
second requirement holds true at $t=\tau_{0}+\sigma_{1}^{T}$, where
$\sigma_{1}^{T}$ is generated from a probability distribution function
$D_{T}(t-\tau_{0})=1-\mbox{exp}(-\int_{\tau_{0}}^{t}{\cal
D}_{T}(C_{s}^{T})ds)$. A Uniform random variable $U_{i}$ is used to determine
this $\sigma_{1}^{T}$. In the case of competing enablings, the following rules
apply:

  * R0

The firing of an immediate transition has priority over the firing of a guard
or a delay transition.

  * R1

If one transition becomes enabled by two or more sets of input tokens at
exactly the same time, and the firing of any one set will not disable one or
more other sets, then it will fire these sets of tokens independently, at the
same time.

  * R2

If one transition becomes enabled by two or more sets of input tokens at
exactly the same time, and the firing of any one set disables one or more
other sets, then the set that is fired is selected randomly, with the same
probability for each set.

  * R3

If two or more transitions become enabled at exactly the same time and the
firing of any one transition will not disable the other transitions, then they
will fire at the same time.

  * R4

If two or more transitions become enabled at exactly the same time and the
firing of any one transition disables some other transitions, then each
combination of transitions that can fire independently without leaving enabled
transitions gets the same probability of firing.

#### Transition firing.

If $T$ is enabled, suppose this occurs at time $\tau_{1}$ and in a particular
vector of token colours $C_{\tau_{1}}^{T}$, it removes one token per ordinary
input arc corresponding with $C_{\tau_{1}}^{T}$ from each of its input places
(i.e. tokens are not removed along enabling arcs). Next, $T$ produces zero or
one token along each output arc: If $(e_{\tau_{1}}^{T},a_{\tau_{1}}^{T})$ is a
random hybrid vector generated from probability measure ${\cal
F}_{T}(\cdot;C_{\tau_{1}}^{T})$ (by making use of a Uniform random variable
$U_{i}$), then vector $e_{\tau_{1}}^{T}$ is a vector of zeros and ones, where
the $i$th vector element corresponds with the $i$th outgoing arc of transition
$T$. An output place gets a token iff it is connected to an arc that
corresponds with a vector element 1. Moreover, $a_{\tau_{1}}^{T}$ specifies
the colours of the produced tokens.

#### Execution from first transition firing onwards.

At $t=\tau_{1}$, zero or more transitions are pre-enabled (if this number is
zero, no transitions will fire anymore). If these include immediate
transitions, then these are fired without delay, but with use of rules R0–R4.
If after this, still immediate transitions are enabled, then these are also
fired, and so forth, until no more immediate transitions are enabled. Next,
the SDCPN is executed in the same way as described above for the situation
from $\tau_{0}$ onwards.

###  2.3 SDCPN stochastic process

The marking of the SDCPN is given by the numbers of tokens in the places and
the associated colour values of these tokens and can be mapped to a
probabilistically unique SDCPN stochastic process $\\{M_{t},C_{t}\\}$ as
follows: For any $t\geq\tau_{0}$, let a token distribution be characterised by
the vector $M_{t}^{\prime}=(M_{1,t}^{\prime},\ldots,M_{|{\cal
P}|,t}^{\prime})$, where $M_{i,t}^{\prime}\in\mathbb{N}$ denotes the number of
tokens in place $P_{i}$ at time $t$ and $1,\ldots,|{\cal P}|$ refers to a
unique ordering of places adopted for SDCPN. At times
$t\in(\tau_{k-1},\tau_{k})$ when no transition fires, the token distribution
is unique and we define $M_{t}=M_{t}^{\prime}$. The associated colours of
these tokens are gathered in a column vector $C_{t}$ which first contains all
colours of tokens in place $P_{1}$, next (i.e. below it) all colours of tokens
in place $P_{2}$, etc, until place $P_{|{\cal P}|}$. If at time $t=\tau_{k}$
one or more transitions fire, then the SDCPN discrete process state at time
${\tau_{k}}$ is defined by $M_{\tau_{k}}=$ the token distribution that occurs
after all transitions that fire at time ${\tau_{k}}$ have been fired. The
associated colours of these tokens are gathered in a column vector
$C_{\tau_{k}}$ in the same way as described above. This construction ensures
that the process $\\{M_{t},C_{t}\\}$ has limits from the left and is
continuous from the right, i.e., it satisfies the càdlàg property.

##  3 Air traffic example and its SDCPN model

To illustrate the advantages of SDCPN when modelling a complex system,
consider a simplified model of the evolution of an aircraft in one sector of
airspace. The deviation of this aircraft from its intended path is affected by
its engine system and its navigation system. Each of these aircraft systems
can be in either Working (functioning properly) or Not working (operating in
some failure mode). Both systems switch between these modes independently and
with exponentially distributed sojourn times, with finite rates $\delta_{3}$
(engine repaired), $\delta_{4}$ (engine fails), $\delta_{5}$ (navigation
repaired) and $\delta_{6}$ (navigation fails), respectively. If both systems
are Working, the aircraft evolves in Nominal mode and the position $Y_{t}$ and
velocity $S_{t}$ of the aircraft are determined by $dX_{t}={\cal
V}_{1}(X_{t})dt+{\cal W}_{1}dW_{t}$, where $X_{t}=(Y_{t},S_{t})^{\prime}$. If
either one, or both, of the systems is Not working, the aircraft evolves in
Non-nominal mode and the position and velocity of the aircraft are determined
by $dX_{t}={\cal V}_{2}(X_{t})dt+{\cal W}_{2}dW_{t}$. The factors ${\cal
W}_{1}$ and ${\cal W}_{2}$ are determined by wind fluctuations. Initially, the
aircraft has position $Y_{0}$ and velocity $S_{0}$, while both its systems are
Working. The evaluation of this process may be stopped when the aircraft has
Landed, i.e. its vertical position and velocity are equal to zero.

An SDCPN graph for this example is developed in two stages. In the first
stage, the agents of the operation are modelled separately, by one local SDCPN
each, see Fig. 3a. In the next stage, the interactions between the agents are
modelled, thus connecting the local SDCPN, Fig. 3b.

Figure 3: SDCPN graph for the aircraft evolution example

Fig. 3b shows the SDCPN graph for this example, where,

  * •

$P_{1}$ denotes aircraft evolution Nominal, i.e. evolution is according to
${\cal V}_{1}$ and ${\cal W}_{1}$.

  * •

$P_{2}$ denotes aircraft evolution Non-nominal, i.e. evolution is according to
${\cal V}_{2}$ and ${\cal W}_{2}$.

  * •

$P_{3}$ and $P_{4}$ denote engine system Not working and Working,
respectively.

  * •

$P_{5}$ and $P_{6}$ denote navigation system Not working and Working,
respectively.

  * •

$P_{7}$ denotes the aircraft has landed.

  * •

$T_{1a}$ and $T_{1b}$ denote a transition of aircraft evolution from Nominal
to Non-nominal, due to engine system or navigation system Not working,
respectively.

  * •

$T_{2}$ denotes a transition of aircraft evolution from Non-nominal to
Nominal, due to engine system and navigation system both Working again.

  * •

$T_{3}$ through $T_{6}$ denote transitions between Working and Not working of
the engine and navigation systems.

  * •

$T_{7}$ and $T_{8}$ denote transitions of the aircraft landing.

The graph in Fig. 3b completely defines SDCPN elements ${\cal P}$, ${\cal T}$,
${\cal A}$ and ${\cal N}$, where ${\cal T}_{G}=\\{T_{7},T_{8}\\}$, ${\cal
T}_{D}=\\{T_{3},T_{4},T_{5},T_{6}\\}$ and ${\cal
T}_{I}=\\{T_{1a},T_{1b},T_{2}\\}$. The other SDCPN elements are specified
below:

  * ${\cal S}$:

Two colour types are defined; ${\cal S}=\\{\mathbb{R}^{0},\mathbb{R}^{6}\\}$.

  * ${\cal C}$:

${\cal C}(P_{1})={\cal C}(P_{2})={\cal C}(P_{7})=\mathbb{R}^{6}$, i.e. tokens
in $P_{1}$, $P_{2}$ and $P_{7}$ have colours in $\mathbb{R}^{6}$; the colour
components model the 3-dimensional position and 3-dimensional velocity of the
aircraft. ${\cal C}(P_{3})={\cal C}(P_{4})={\cal C}(P_{5})={\cal
C}(P_{6})=\mathbb{R}^{0}\triangleq\emptyset$.

  * ${\cal I}$:

Place $P_{1}$ initially has a token with colour
$X_{0}=(Y_{0},S_{0})^{\prime}$, with $Y_{0}\in\mathbb{R}^{2}\times(0,\infty)$
and $S_{0}\in\mathbb{R}^{3}\setminus\mbox{Col}\\{0,0,0\\}$. Places $P_{4}$ and
$P_{6}$ initially each have a token without colour.

  * ${\cal V}$,

${\cal W}$: The token colour functions for places $P_{1}$, $P_{2}$ and $P_{7}$
are determined by $({\cal V}_{1},{\cal W}_{1})$, $({\cal V}_{2},{\cal
W}_{2})$, and $({\cal V}_{7},{\cal W}_{7})$, respectively, where $({\cal
V}_{7},{\cal W}_{7})=(0,0)$. For places $P_{3}$ – $P_{6}$ there is no token
colour function.

  * ${\cal G}$:

Transitions $T_{7}$ and $T_{8}$ have a guard defined by ${\cal
G}_{T_{7}}={\cal
G}_{T_{8}}=\mathbb{R}^{2}\times(0,\infty)\times\mathbb{R}^{2}\times(0,\infty)$.

  * ${\cal D}$:

The jump rates for transitions $T_{3}$, $T_{4}$, $T_{5}$ and $T_{6}$ are
${\cal D}_{T_{3}}(\cdot)=\delta_{3}$, ${\cal D}_{T_{4}}(\cdot)=\delta_{4}$,
${\cal D}_{T_{5}}(\cdot)=\delta_{5}$ and ${\cal D}_{T_{6}}(\cdot)=\delta_{6}$.

  * ${\cal F}$:

Each transition has a unique output place, to which it fires a token with a
colour (if applicable) equal to the colour of the token removed.

##  4 From SDCPN to GSHS

Following [8], this section first presents a definition of _general stochastic
hybrid system_ (GSHS) and its execution. In [16] it has been proven that under
a few conditions, SDCPN and GSHS are bisimilar. In Subsection 4.2 this is
illustrated by showing how the SDCPN example of the previous section can be
mapped to a bisimilar GSHS.

###  4.1 GSHS definition

#### Definition 4.1

(General stochastic hybrid system) _A GSHS is an automaton ( ${\bf K}$, $d$,
${\cal X}$, $f$, $g$, _Init_ , $\lambda$, $Q$), where_

  * •

${\bf K}$_is a countable set._

  * •

$d:{\bf K}\rightarrow\mathbb{N}$_maps each_ $\theta\in{\bf K}$ to a natural
number.

  * •

${\cal X}:{\bf K}\rightarrow\\{E_{\theta};\theta\in{\bf K}\\}$_maps each_
$\theta\in{\bf K}$ to an open subset $E_{\theta}$ of $\mathbb{R}^{d(\theta)}$.
With this, the hybrid state space is given by
$E\triangleq\\{\\{\theta\\}\times E_{\theta};\theta\in{\bf K}\\}$.

  * •

$f:E\to\\{\mathbb{R}^{d(\theta)};\theta\in{\bf K}\\}$_is a vector field._

  * •

$g:E\to\\{\mathbb{R}^{d(\theta)\times h};\theta\in{\bf K}\\}$_is a matrix
field, with_ $h\in\mathbb{N}$.

  * •

_Init_ $:{\cal B}(E)\rightarrow[0,1]$_is an initial probability measure, with_
${\cal B}(E)$ the Borel $\sigma$-algebra on $E$.

  * •

$\lambda:E\to\mathbb{R}^{+}$_is a jump rate function._

  * •

$Q:\mathcal{B}(E)\times(E\cup\partial E)\rightarrow[0,1]$_is a GSHS transition
measure, where_ $\partial E\triangleq\\{\\{\theta\\}\times\partial
E_{\theta};\theta\in{\bf K}\\}$ is the boundary of $E$, in which $\partial
E_{\theta}$ is the boundary of $E_{\theta}$.

#### Definition 4.2

(GSHS execution) _A stochastic process $\\{\theta_{t},X_{t}\\}$ is called a
_GSHS execution_ if there exists a sequence of stopping times
$0=\tau_{0}<\tau_{1}<\tau_{2}\cdots$ such that for each $k\in\mathbb{N}$:_

  * •

$(\theta_{0},X_{0})$_is an_ $E$-valued random variable extracted according to
probability measure _Init_.

  * •

_For_ $t\in[\tau_{k},\tau_{k+1})$, $\theta_{t}=\theta_{\tau_{k}}$ and
$X_{t}=X_{t}^{k}$, where for $t\geq\tau_{k}$, $X_{t}^{k}$ is a solution of the
stochastic differential equation
$dX_{t}^{k}=f(\theta_{\tau_{k}},X_{t}^{k})dt+g(\theta_{\tau_{k}},X_{t}^{k})dB_{t}^{\theta_{\tau_{k}}}$
with initial condition $X_{\tau_{k}}^{k}=X_{\tau_{k}}$, and where
$\\{B_{t}^{\theta}\\}$ is $h$-dimensional standard Brownian motion for each
$\theta\in{\bf K}$.

  * •

$\tau_{k+1}=\tau_{k}+\sigma_{k}$_, where_ $\sigma_{k}$ is chosen according to
a survivor function given by $F(t)=$
${\bf 1}_{(t<\tau^{*})}\exp(-\int_{0}^{t}\lambda(\theta,X_{s}^{k})ds)$. Here,
$\tau^{*}=\inf\\{t>\tau_{k}\mid X_{t}^{k}\in\partial E_{\theta_{\tau_{k}}}\\}$
and ${\bf 1}$ is indicator function.

  * •

_The probability distribution of_ $(\theta_{\tau_{k+1}},X_{\tau_{k+1}})$, i.e.
the hybrid state right after the jump, is governed by the law
$Q(\cdot;(\theta_{\tau_{k}},X_{\tau_{k+1}-}))$.

[8] show that under assumptions G1-G4 below, a GSHS execution is a strong
Markov Process and has the càdlàg property (right continuous with left hand
limits).

  * G1

$f(\theta,\cdot)$ and $g(\theta,\cdot)$ are Lipschitz continuous and bounded.
This yields that for each initial state $(\theta,x)$ at initial time $\tau$
there exists a pathwise unique solution $X_{t}$ to
$dX_{t}=f(\theta,X_{t})dt+g(\theta,X_{t})dB_{t}$, where $\\{B_{t}\\}$ is
$h$-dimensional standard Brownian motion.

  * G2

$\lambda:E\to\mathbb{R}^{+}$ is a measurable function such that for all
$\xi\in E$, there is $\epsilon(\xi)>0$ such that
$t\to\lambda(\theta_{t},X_{t})$ is integrable on $[0,\epsilon(\xi))$.

  * G3

For each fixed $A\in{\cal B}(E)$, the map $\xi\to Q(A;\xi)$ is measurable and
for any $(\theta,x)\in E\cup\partial E$, $Q(\cdot;\theta,x)$ is a probability
measure.

  * G4

If $N_{t}=\sum_{k}{\bf 1}_{(t\geq\tau_{k})}$, then it is assumed that for
every starting point $(\theta,x)$ and for all $t\in\mathbb{R}^{+}$,
$\mathbb{E}N_{t}<\infty$. This means, there will be a finite number of jumps
in finite time.

###  4.2 A bisimilar GSHS for the example SDCPN

Next we transform the SDCPN example model of Section 3 into a bisimilar GSHS.
The first step is to construct the state space ${\bf K}$ for the GSHS discrete
process $\\{\theta_{t}\\}$. This is done by identifying the SDCPN reachability
graph. Nodes in the reachability graph provide the number of tokens in each of
the SDCPN places. Arrows connect these nodes as they represent transitions
firing. The SDCPN of Fig. 3b has seven places hence the reachability graph for
this example has elements that are vectors of length 7. These nodes, excluding
the nodes that enable immediate transitions, form the GSHS discrete state
space.

Figure 4: Reachability graph for the SDCPN of Fig. 3b. The nodes in bold type
face correspond with the elements of the GSHS discrete state space ${\bf K}$.

The reachability graph is shown in Fig. 4, with nodes that form the GSHS
discrete state space in Bold typeface, i.e. ${\bf
K}=\\{V_{1},\ldots,V_{8}\\}$, with $V_{1}=(1,0,0,1,0,1,0)$,
$V_{2}=(0,1,1,0,0,1,0)$, $V_{3}=(0,1,1,0,1,0,0)$, $V_{4}=(0,1,0,1,1,0,0)$,
$V_{5}=(0,0,0,1,0,1,1)$, $V_{6}=(0,0,1,0,0,1,1)$, $V_{7}=(0,0,1,0,1,0,1)$,
$V_{8}=(0,0,0,$ $1,1,0,1)$. Since initially there is a token in places
$P_{1}$, $P_{4}$ and $P_{6}$, the GSHS initial mode equals
$\theta_{0}=V_{1}=(1,0,0,1,0,1,0)$. The GSHS initial continuous state value
equals the vector containing the initial colours of all initial tokens. Since
the initial colour of the token in Place $P_{1}$ equals $X_{0}$, and the
tokens in places $P_{4}$ and $P_{6}$ have no colour, the GSHS initial
continuous state value equals Col$\\{X_{0},\emptyset,\emptyset\\}=X_{0}$. The
GSHS drift coefficient $f$ is given by $f(\theta,\cdot)={\cal V}_{1}(\cdot)$
for $\theta=V_{1}$, $f(\theta,\cdot)={\cal V}_{2}(\cdot)$ for
$\theta\in\\{V_{2},V_{3},V_{4}\\}$, and $f(\theta,\cdot)=0$ otherwise. For the
diffusion coefficient, $g(\theta,\cdot)={\cal W}_{1}$ for $\theta=V_{1}$,
$g(\theta,\cdot)={\cal W}_{2}$ for $\theta\in\\{V_{2},V_{3},V_{4}\\}$, and
$g(\theta,\cdot)=0$ otherwise. The hybrid state space is given by
$E=\\{\\{\theta\\}\times E_{\theta};\theta\in\mathbb{M}\\}$, where for
$\theta\in\\{V_{1},V_{2},V_{3},V_{4}\\}$:
$E_{\theta}=\mathbb{R}^{2}\times(0,\infty)\times\mathbb{R}^{2}\times(0,\infty)$
and for $\theta\in\\{V_{5},V_{6},V_{7},V_{8}\\}$: $E_{\theta}=\mathbb{R}^{6}$.
Always two delay transitions are pre-enabled: either $T_{3}$ or $T_{4}$ and
either $T_{5}$ or $T_{6}$. This yields
$\lambda(V_{1},\cdot)=\lambda(V_{5},\cdot)=\delta_{4}+\delta_{6}$,
$\lambda(V_{2},\cdot)=\lambda(V_{6},\cdot)=\delta_{3}+\delta_{6}$,
$\lambda(V_{3},\cdot)=\lambda(V_{7},\cdot)=\delta_{3}+\delta_{5}$,
$\lambda(V_{4},\cdot)=\lambda(V_{8},\cdot)=\delta_{4}+\delta_{5}$. For the
determination of GSHS transition measure $Q$, we make use of the reachability
graph, the sets ${\cal D}$, ${\cal G}$ and ${\cal F}$ and the rules R0–R4. In
Table 1, $Q(\theta^{\prime},x^{\prime};\theta,x)=p$ denotes that if
$(\theta,x)$ is the value of the GSHS state before the hybrid jump, then, with
probability $p$, $(\theta^{\prime},x^{\prime})$ is the value of the GSHS state
immediately after the jump.

Table 1: Example GSHS transition measure for size of jump For $x\notin\partial E_{V_{1}}$: |  $Q(V_{2},x;V_{1},x)=\frac{\delta_{4}}{\delta_{4}+\delta_{6}}$, | $Q(V_{4},x;V_{1},x)=\frac{\delta_{6}}{\delta_{4}+\delta_{6}}$
---|---|---
For $x\in\partial E_{V_{1}}$: | $Q(V_{5},x;V_{1},x)=1$ |
For $x\notin\partial E_{V_{2}}$: |  $Q(V_{3},x;V_{2},x)=\frac{\delta_{6}}{\delta_{3}+\delta_{6}}$, | $Q(V_{1},x;V_{2},x)=\frac{\delta_{3}}{\delta_{3}+\delta_{6}}$
For $x\in\partial E_{V_{2}}$: | $Q(V_{6},x;V_{2},x)=1$ |
For $x\notin\partial E_{V_{3}}$: |  $Q(V_{4},x;V_{3},x)=\frac{\delta_{3}}{\delta_{3}+\delta_{5}}$, | $Q(V_{2},x;V_{3},x)=\frac{\delta_{5}}{\delta_{3}+\delta_{5}}$
For $x\in\partial E_{V_{3}}$: | $Q(V_{7},x;V_{3},x)=1$ |
For $x\notin\partial E_{V_{4}}$: |  $Q(V_{3},x;V_{4},x)=\frac{\delta_{4}}{\delta_{4}+\delta_{5}}$, | $Q(V_{1},x;V_{4},x)=\frac{\delta_{5}}{\delta_{4}+\delta_{5}}$
For $x\in\partial E_{V_{4}}$: | $Q(V_{8},x;V_{4},x)=1$ |
For all $x$: |  $Q(V_{6},x;V_{5},x)=\frac{\delta_{4}}{\delta_{4}+\delta_{6}}$, | $Q(V_{8},x;V_{5},x)=\frac{\delta_{6}}{\delta_{4}+\delta_{6}}$
For all $x$: |  $Q(V_{7},x;V_{6},x)=\frac{\delta_{6}}{\delta_{3}+\delta_{6}}$, | $Q(V_{5},x;V_{6},x)=\frac{\delta_{3}}{\delta_{3}+\delta_{6}}$
For all $x$: |  $Q(V_{8},x;V_{7},x)=\frac{\delta_{3}}{\delta_{3}+\delta_{5}}$, | $Q(V_{6},x;V_{7},x)=\frac{\delta_{5}}{\delta_{3}+\delta_{5}}$
For all $x$: |  $Q(V_{7},x;V_{8},x)=\frac{\delta_{4}}{\delta_{4}+\delta_{5}}$, | $Q(V_{5},x;V_{8},x)=\frac{\delta_{5}}{\delta_{4}+\delta_{5}}$

With this, the SDCPN of the aircraft evolution example is uniquely mapped to
an GSHS. It can be shown that the SDCPN execution and the execution of the
resulting GSHS are probabilistically equivalent, i.e. the SDCPN and the GSHS
are bisimilar. Thanks to this bisimilarity we can now use the automata
framework to analyse the GSHP that is defined by the SDCPN model for the
example.

##  5 From SDCPN to HSDE

Following [2] and [3], this section first presents a definition of _hybrid
stochastic differential equation_ (HSDE) and gives conditions under which the
HSDE has a pathwise unique solution. This pathwise unique solution is referred
to as _HSDE solution process_ or GSHP. The basic advantage of using HSDE in
defining a GSHP over using GSHS is that with the HSDE approach the spontaneous
jump mechanism is explicitly built on an underlying stochastic basis, whereas
in GSHS the execution itself imposes an underlying stochastic basis. In [17]
it has been proven that under a few conditions, SDCPN and HSDE are bisimilar.
In Subsection 5.2 this is illustrated by showing how the SDCPN example of the
previous section can be mapped to a bisimilar HSDE.

###  5.1 HSDE definition

For the HSDE setting we start with a complete stochastic basis
$(\Omega,\Im,\mathbb{F},\mathbb{P},\mathbb{T})$, in which a complete
probability space $(\Omega,\Im,\mathbb{P})$ is equipped with a right-
continuous filtration $\mathbb{F}=\\{\Im_{t}\\}$ on the positive time line
$\mathbb{T}=\mathbb{R}^{+}$. This stochastic basis is endowed with a
probability measure $\mu_{\theta_{0},X_{0}}$ for the initial state, an
independent $h$-dimensional standard Wiener process $\\{W_{t}\\}$ and an
independent homogeneous Poisson random measure $p_{P}(dt,dz)$ on
$\mathbb{T}\times\mathbb{R}^{d+1}$.

#### Definition 5.1

(Hybrid stochastic differential equation) _An HSDE on stochastic basis
$(\Omega,\Im,\mathbb{F},\mathbb{P},\mathbb{T})$, is defined as a set of
equations (1)-(8) in which a collection of elements ($\mathbb{M}$, $E$, $f$,
$g$, $\mu_{\theta_{0},X_{0}}$, $\Lambda$, $\psi$, $\rho$, $\mu$, $p_{P}$,
$\\{W_{t}\\}$) appear._
The elements ($\mathbb{M}$, $E$, $f$, $g$, $\mu_{\theta_{0},X_{0}}$,
$\Lambda$, $\psi$, $\rho$, $\mu$, $p_{P}$, $\\{W_{t}\\}$) are defined as
follows:

  * •

$\mathbb{M}=\\{\vartheta_{1},\ldots,\vartheta_{N}\\}$ is a finite set,
$N\in\mathbb{N}$, $1\leq N<\infty$.

  * •

$E=\\{\\{\theta\\}\times E_{\theta};\theta\in\mathbb{M}\\}$ is the hybrid
state space, where for each $\theta\in\mathbb{M}$, $E_{\theta}$ is an open
subset of $\mathbb{R}^{n}$ with boundary $\partial E_{\theta}$. The boundary
of $E$ is $\partial E=\\{\\{\theta\\}\times\partial
E_{\theta};\theta\in\mathbb{M}\\}$.

  * •

$f:\mathbb{M}\times\mathbb{R}^{n}\to\mathbb{R}^{n}$ is a measurable mapping.

  * •

$g:\mathbb{M}\times\mathbb{R}^{n}\to\mathbb{R}^{n\times h}$ is a measurable
mapping.

  * •

$\mu_{\theta_{0},X_{0}}:\Omega\times{\cal B}(E)\to[0,1]$ is a probability
measure for the initial random variables $\theta_{0}$, $X_{0}$, which are
defined on the stochastic basis; $\mu_{\theta_{0},X_{0}}$ is assumed to be
invertible.

  * •

$\Lambda:\mathbb{M}\times\mathbb{R}^{n}\to[0,\infty)$ is a measurable mapping.

  * •

$\psi:\mathbb{M}\times\mathbb{M}\times\mathbb{R}^{n}\times\mathbb{R}^{d}\to\mathbb{R}^{n}$
is a measurable mapping such that $x+\psi(\vartheta,\theta,x,\underline{z})\in
E_{\vartheta}$ for all $x\in E_{\theta}$, $\underline{z}\in\mathbb{R}^{d}$,
and $\vartheta,\theta\in\mathbb{M}$.

  * •

$\rho:\mathbb{M}\times\mathbb{M}\times\mathbb{R}^{n}\to[0,\infty)$ is a
measurable mapping such that $\sum_{i=1}^{N}\rho(\vartheta_{i},\theta,x)=1$
for all $\theta\in\mathbb{M},x\in\mathbb{R}^{n}$.

  * •

$\mu:\Omega\times\mathbb{R}^{d}\to[0,1]$ is a probability measure which is
assumed to be invertible.

  * •

$p_{P}:\Omega\times\mathbb{T}\times\mathbb{R}^{d+1}\to\\{0,1\\}$ is a
homogeneous Poisson random measure on the stochastic basis, independent of
$(\theta_{0},X_{0})$. The intensity measure of $p_{P}(dt,dz)$ equals
$dt\cdot\mu_{L}(dz_{1})\cdot\mu(d\underline{z})$, where
$z=\mbox{Col}\\{z_{1},\underline{z}\\}$ and $\mu_{L}$ is the Lebesgue measure.

  * •

$W:\Omega\times\mathbb{T}\to\mathbb{R}^{h}$ such that $\\{W_{t}\\}$ is an
$h$-dimensional standard Wiener process on the stochastic basis, and
independent of $(\theta_{0},X_{0})$ and $p_{P}$.

Using these elements, the HSDE process $\\{\theta_{t}^{*},X_{t}^{*}\\}$ is
defined as follows:

$\theta_{t}^{*}=\theta_{t}^{k}\mbox{ for all
}t\in[\tau_{k}^{b},\tau_{k+1}^{b}),k=0,1,2,\ldots$ (1)
$X_{t}^{*}=X_{t}^{k}\mbox{ for all
}t\in[\tau_{k}^{b},\tau_{k+1}^{b}),k=0,1,2,\ldots$ (2)

Hence $\\{\theta_{t}^{*},X_{t}^{*}\\}$ consists of a concatenation of
processes $\\{\theta_{t}^{k},X_{t}^{k}\\}$ which are defined by (3)-(8) below.
If the system (1)-(8) has a solution in probabilistic sense, then the process
$\\{\theta_{t}^{*},X_{t}^{*}\\}$ is referred to as _HSDE solution process_ or
_GSHP_.

$d\theta_{t}^{k}=\sum_{i=1}^{N}(\vartheta_{i}-\theta_{t-}^{k})p_{P}(dt,(\Sigma_{i-1}(\theta_{t-}^{k},X_{t-}^{k}),\Sigma_{i}(\theta_{t-}^{k},X_{t-}^{k})]\times\mathbb{R}^{d})$
(3)
$dX_{t}^{k}=f(\theta_{t}^{k},X_{t}^{k})dt+g(\theta_{t}^{k},X_{t}^{k})dW_{t}+\int_{\mathbb{R}^{d}}\psi(\theta_{t}^{k},\theta_{t-}^{k},X_{t-}^{k},\underline{z})p_{P}(dt,(0,\Lambda(\theta_{t-}^{k},X_{t-}^{k})]\times
d\underline{z})$ (4)

with $\theta_{0}^{0}=\theta_{0}$, $X_{0}^{0}=X_{0}$ and with $\Sigma_{0}$
through $\Sigma_{N}$ measurable mappings satisfying, for
$\theta\in\mathbb{M}$, $\vartheta_{j}\in\mathbb{M}$, $x\in\mathbb{R}^{n}$:

$\Sigma_{i}(\theta,x)=\left\\{\begin{array}[]{ll}\Lambda(\theta,x)\sum_{j=1}^{i}\rho(\vartheta_{j},\theta,x)&\mbox{if
}i>0\\\ 0&\mbox{if }i=0\end{array}\right.$ (5)

In addition, for $k=0,1,2,\ldots$, with $\tau_{0}^{b}=0$:

$\tau_{k+1}^{b}\triangleq\inf\\{t>\tau_{k}^{b}\mid(\theta_{t}^{k},X_{t}^{k})\in\partial
E\\}$ (6)
$\mathbb{P}\\{\theta_{\tau_{k+1}^{b}}^{k+1}=\vartheta,X_{\tau_{k+1}^{b}}^{k+1}\in{A}\mid\theta_{\tau_{k+1}^{b}-}^{k}=\theta,X_{\tau_{k+1}^{b}-}^{k}=x\\}=Q(\\{\vartheta\\}\times
A;\theta,x)$ (7)

for $A\in{\cal B}(\mathbb{R}^{n})$, where $Q$ is given by

$Q(\\{\vartheta\\}\times
A;\theta,x)=\rho(\vartheta,\theta,x)\int_{\mathbb{R}^{d}}{\bf
1}_{A}(x+\psi(\vartheta,\theta,x,\underline{z}))\mu(d\underline{z})$ (8)

Next, the following proposition can be shown to hold true [17]:

#### Proposition 5.1

_Let conditions H1-H8 below hold true. Let
$(\theta_{0}^{*}(\omega),X_{0}^{*}(\omega))=(\theta_{0},X_{0})\in E$ for all
$\omega$. Then for every initial condition $(\theta_{0},X_{0})$, (1)-(8) has a
pathwise unique solution $\\{\theta_{t}^{*},X_{t}^{*}\\}$ which is càdlàg and
adapted and is a semi-martingale assuming values in the hybrid state space
$E$._

H1

For all $\theta\in\mathbb{M}$ there exists a constant $K(\theta)$ such that
for all $x\in\mathbb{R}^{n}$, $|f(\theta,x)|^{2}+\|g(\theta,x))\|^{2}\leq
K(\theta)(1+|x|^{2})$, where $|a|^{2}=\sum_{i}(a_{i})^{2}$ and
$||b||^{2}=\sum_{i,j}(b_{ij})^{2}$.

H2

For all $r\in\mathbb{N}$ and for all $\theta\in\mathbb{M}$ there exists a
constant $L_{r}(\theta)$ such that for all $x$ and $y$ in the ball
$B_{r}=\\{z\in\mathbb{R}^{n}\mid|z|\leq r+1\\}$,
$|f(\theta,x)-f(\theta,y)|^{2}+\|g(\theta,x)-g(\theta,y)\|^{2}\leq
L_{r}(\theta)|x-y|^{2}$.

H3

For each $\theta\in\mathbb{M}$, the mapping
$\Lambda(\theta,\cdot):\mathbb{R}^{n}\to[0,\infty)$ is continuous and bounded,
with upper bound a constant $C_{\Lambda}$.

H4

For all $(\theta,\vartheta)\in\mathbb{M}^{2}$, the mapping
$\rho(\vartheta,\theta,\cdot):\mathbb{R}^{n}\to[0,\infty)$ is continuous.

H5

For all $r\in\mathbb{N}$ there exists a constant $M_{r}(\theta)$ such that

$\sup_{|x|\leq
r}\int_{\mathbb{R}^{d}}|\psi(\vartheta,\theta,x,\underline{z})|\mu(d\underline{z})\leq
M_{r}(\theta),\mbox{ for all }\vartheta,\theta\in\mathbb{M}$

H6

$|\psi(\theta,\theta,x,\underline{z})|=0$ or $>1$ for all
$\theta\in\mathbb{M}$, $x\in\mathbb{R}^{n}$, $\underline{z}\in\mathbb{R}^{d}$

H7

$\\{(\theta_{t}^{*},X_{t}^{*})\\}$ hits the boundary $\partial E$ a finite
number of times on any finite time interval

H8

$|\vartheta_{i}-\vartheta_{j}|>1$ for $i\neq j$, with $|\cdot|$ a suitable
metric well defined on $\mathbb{M}$.

###  5.2 A bisimilar HSDE for the example SDCPN

Next we transform the SDCPN example model of Section 3 into a bisimilar HSDE.
This mapping follows much the same procedure as for SDCPN to GSHS, except that
the discrete state space is now referred to as $\mathbb{M}$ (rather than ${\bf
K}$) and the Markov jump rate is now referred to as $\Lambda$ (rather than
$\lambda$). The main additional difference is that the HSDE elements do not
include a transition measure $Q$ to define the size of jump, but include
functions $\psi$, $\rho$ and $\mu$ instead. The mapping of SDCPN to HSDE uses
the construction of transition measure $Q$ as an intermediate step, however.
For the particular example SDCPN in this paper, these functions can be
determined from $Q$ as follows: Since the continuous valued process jumps to
the same value with probability 1, we find that
$\psi(V^{i},V^{j},x,\underline{z})=0$ for all $V^{i}$, $V^{j}$, $x$,
$\underline{z}$. Moreover, $\rho(V^{i},V^{j},x)=P_{Q}(V^{i},x,V^{j},x)$ and
$\mu$ may be any given invertible probability measure.

With this, the SDCPN of the aircraft evolution example is uniquely mapped to
an HSDE. If in addition, we want to make use of the HSDE properties of
Proposition 5.1, i.e. the resulting HSDE solution process being adapted and a
semi-martingale, we need to make sure that HSDE conditions H1-H8 are
satisfied. It is shown below that they are, under the following sufficient
condition D1 for the example SDCPN.

D1

For $P\in\\{P_{1},P_{2}\\}$, there exist $K_{P}^{v}$, $L_{P}^{v}$, $K_{P}^{w}$
and $L_{P}^{w}$ such that for all $c,a\in{\cal C}(P)$,
$|{\cal V}_{P}(c)|^{2}\leq K_{P}^{v}(1+|c|^{2})$ and $|{\cal V}_{P}(c)-{\cal
V}_{P}(a)|^{2}\leq L_{P}^{v}|c-a|^{2}$ and
$\|{\cal W}_{P}(c)\|^{2}\leq K_{P}^{w}(1+|c|^{2})$ and $\|{\cal
W}_{P}(c)-{\cal W}_{P}(a)\|^{2}\leq L_{P}^{w}|c-a|^{2}$.

We verify that under condition D1, HSDE conditions H1-H8 hold true in this
example.

H1:

From the construction of $f$ and $g$ above we have for $\theta=V_{1}$:
$|f(\theta,x)|^{2}+\|g(\theta,x)\|^{2}=|{\cal V}_{1}(x)|^{2}+\|{\cal
W}_{1}(x)\|^{2}\leq
K_{P_{1}}^{v}(1+|x|^{2})+K_{P_{1}}^{w}(1+|x|^{2})=K(\theta)(1+|x|^{2})$, with
$K(\theta)=(K_{P_{1}}^{v}+K_{P_{1}}^{w})$. For $\theta=V_{2},V_{3},V_{4}$ the
verification is with replacing ${\cal V}_{1}$, ${\cal W}_{1}$ by ${\cal
V}_{2}$, ${\cal W}_{2}$.

H2:

From the construction of $f$ and $g$ above we have for $\theta=V_{1}$:
$|f(\theta,x)-f(\theta,y)|^{2}+\|g(\theta,x)-g(\theta,y)\|^{2}=|{\cal
V}_{1}(x)-{\cal V}_{1}(y)|^{2}+\|{\cal W}_{1}(x)-{\cal W}_{1}(y)\|^{2}\leq
L_{P_{1}}^{v}|x-y|^{2}+L_{P_{1}}^{w}|x-y|^{2}=L_{r}(\theta)|x-y|^{2}$ with
$L_{r}(\theta)=L_{P_{1}}^{v}+L_{P_{1}}^{w}$. For $\theta=V_{2},V_{3},V_{4}$
replace ${\cal V}_{1}$, ${\cal W}_{1}$ by ${\cal V}_{2}$, ${\cal W}_{2}$.

H3:

Since $\delta_{3}$–$\delta_{6}$ are constant, for all $\theta$,
$\Lambda(\theta,\cdot)$ is bounded and continuous, with upper bound
$C_{\Lambda}=\max\\{\delta_{4}+\delta_{6},\delta_{3}+\delta_{6},\delta_{3}+\delta_{5},\delta_{4}+\delta_{5}\\}$.

H4:

Since for all $\theta,\vartheta$, $P_{Q}(\vartheta,\cdot;\theta,x)$ is
constant, we find $\rho(\vartheta,\theta,x)=P_{Q}(\vartheta,x,\theta,x)$ is
continuous.

H5 and H6:

These are satisfied due to $\psi(V^{i},V^{j},x,\underline{z})=0$ for all
$V^{i}$, $V^{j}$, $x$, $\underline{z}$.

H7:

This condition holds due to $\delta_{3}$–$\delta_{6}$ being finite and the
fact that in this SDCPN example, there is no firing sequence of more than one
guard transition.

H8:

This condition holds for all $V_{1},\ldots,V_{8}$, with metric
$|a|^{2}=\sum_{i}(a_{i})^{2}$.

Thanks to this bisimilarity mapping we can now use HSDE tools to analyse the
GSHP that is defined by the execution of the SDCPN model for the example.

##  6 Conclusions

The aim of this paper was to explain bisimilarity relations between SDCPN
(stochastically and dynamically coloured Petri net), GSHS (general stochastic
hybrid system) and HSDE (hybrid stochastic differential equation), which means
that the strengths of one stochastic model formalism can be used by both of
the other stochastic model formalisms. More specifically, these bisimilarity
relations make it possible to combine the formal verification power of
automata with the analysis power of stochastic differential equations and the
compositional specification power of Petri nets.

We started in Section 2 by defining SDCPN and the resulting SDCPN stochastic
process, which is referred to as a GSHP (general stochastic hybrid process).
In Section 3 we presented a simple but illustrative SDCPN example model. In
Section 4 we studied GSHP as an execution of a GSHS and illustrated by using
the example of Section 3 that SDCPN and GSHS are bisimilar. Next, in Section 5
we studied GSHP as a stochastic process solution of HSDE and showed with an
illustrative example that SDCPN and HSDE are bisimilar.

The bisimilarities between SDCPN, GSHS and HSDE models for the example
considered mean that the resulting example model inherits the strengths of all
three formal stochastic modelling formalisms. This has been depicted in Fig. 1
in the introduction. Examples of GSHP properties are convergence in
discretisation, existence of limits, existence of event probabilities, strong
Markov properties, reachability analysis. Examples of GSHS features are their
connection to formal methods in automata theory and optimal control theory.
Examples of HSDE features are stochastic analysis tools for semi-martingales.
Examples of SDCPN features are natural expression of causal dependencies,
concurrency and synchronisation mechanism, hierarchical and modular
construction, and graphical representation. These complementary advantages of
SDCPN, GSHS, HSDE and GSHP perspectives tend to increase with the complexity
of the system considered.

An illustrative large scale application of bisimularity relations between
SDCPN, HSDE and stochastic hybrid automata has been developed in air traffic
management. Currently pilots depend of air traffic controllers in solving
potential conflicts between their flight trajectories. This places a huge
requirement on the tasks of an air traffic controller. Imagine a similar kind
of approach for road traffic; then each car driver would be blind and depends
of instructions that some road traffic controller is communicating with the
car drivers. How many cars do you think can be managed by one road traffic
controller? The number of aircraft that one air traffic controller can handle
ranges between 10 and 20, depending of the complexity of the traffic pattern.
Over a decade ago, it had been suggested by [23] that this limitation of the
air traffic controller can be solved by moving the responsibility of conflict
resolution from the air traffic controller to the pilots. Since then this
airborne self separation idea has received a lot of research attention.
Nevertheless, it still is unknown how much more air traffic can safely be
accommodated under a well designed airborne self separation way of working. In
order to add to the solution of this debate, a series of large European
studies towards solving this question have been started under the name
HYBRIDGE [19] and iFly [20] respectively. The way of working is to first
develop a well defined SDCPN model of the airborne self separation concept of
operation to be analysed, e.g. [18]. Subsequently this SDCPN model is further
analysed using a bisimilar HSDE and hybrid automation formal model
representation [6, 4], in which powerful stochastic analysis theory is
exploited for the speeding up of Monte Carlo simulations. Using this approach,
[5] have shown that the first generation of airborne self separation concept
designs falls short in safely accommodating higher air traffic demand than
conventional ATM can. The feedback of this finding to advanced air traffic
concept designers triggered the development of more advanced airborne self
separation concept of operation, e.g. see [20].

## References

  * [1]
  * [2] H.A.P. Blom (2003): _Stochastic hybrid processes with hybrid jumps_.  In: Proceedings IFAC conference on analysis and design of hybrid system (ADHS), Saint-Malo, Brittany, France. pp. 361–365.
  * [3] H.A.P. Blom, G.J. Bakker, M.H.C. Everdij & M.N.J. Van der Park (2003): _Stochastic analysis background of accident risk assessment for air traffic management_.  HYBRIDGE Report, D2.2.  http://hosted.nlr.nl/public/hosted-sites/hybridge/.
  * [4] H.A.P. Blom, G.J. Bakker & J. Krystul (2009): _Rare event estimation for a large scale stochastic hybrid system with air traffic application_.  In: G. Rubino & B. Tuffin, editors: Rare event simulation using Monte Carlo methods, forthcoming. J.Wiley.
  * [5] H.A.P. Blom, B. Klein Obbink & G.J. Bakker (2009): _Simulated collision risk of an uncoordinated airborne self separation concept of operation_.  ATC Quarterly 17, pp. 63–93.
  * [6] H.A.P. Blom, J. Krystul, G.J. Bakker, M.B. Klompstra & B. Klein Obbink (2007): _Free flight collision risk estimation by sequential Monte Carlo simulation_.  In: C.G. Cassandras & J. Lygeros, editors: Stochastic hybrid systems: recent developments and research trends, Control engineering 24, chapter 10. Taylor & Francis Group / CRC Press, pp. 247–281.
  * [7] M.L. Bujorianu & J. Lygeros (2004): _General stochastic hybrid systems: modelling and optimal control_.  In: Proceedings 43rd conference on decision and control (CDC), Nassau, Bahamas.
  * [8] M.L. Bujorianu & J. Lygeros (2006): _Toward a general theory of stochastic hybrid systems_.  In: H.A.P. Blom & J. Lygeros, editors: Stochastic hybrid systems: theory and safety critical applications, Lectures notes in control and information sciences (LNCIS) 337\. Springer, pp. 3–30.
  * [9] M.L. Bujorianu, J. Lygeros & M.C. Bujorianu (2005): _Different approaches on bisimulation for stochastic hybrid systems_.  In: M. Morari & L. Thiele, editors: Proceedings 8th international workshop on hybrid systems: computation and control (HSCC), Zürich, Switzerland, Lecture notes in computer science (LNCS) 3414\. pp. 198–214.
  * [10] M.H.A. Davis (1993): _Markov models and optimization_ , Monographs on statistics and applied probability 49\.  Chapman and Hall, London.
  * [11] R.J. Elliott (1982): _Stochastic calculus and applications_ , Applications of mathematics: Stochastic modelling and applied probability 18\.  Springer-Verlag.
  * [12] R.J. Elliott, L. Aggoun & J.B. Moore (1995): _Hidden Markov models: estimation and control_ , Applications of mathematics: stochastic modelling and applied probability 29\.  Springer-Verlag.
  * [13] S.N. Ethier & T.G. Kurtz (1986): _Markov processes, characterization and convergence_.  Wiley series in probability and mathematical statistics. John Wiley & Sons, New York.
  * [14] M.H.C. Everdij & H.A.P. Blom (2003): _Petri nets and hybrid state Markov processes in a power-hierarchy of dependability models_.  In: Proceedings IFAC conference on analysis and design of hybrid system (ADHS), Saint-Malo, Brittany, France. pp. 355–360.
  * [15] M.H.C. Everdij & H.A.P. Blom (2005): _Piecewise deterministic Markov processes represented by dynamically coloured Petri nets_.  In: S. Jacka, editor: Stochastics: an international journal of probability and stochastic processes, 77, number 1. Taylor & Francis, pp. 1–29.
  * [16] M.H.C. Everdij & H.A.P. Blom (2006): _Hybrid Petri nets with diffusion that have into-mappings with generalised stochastic hybrid processes_.  In: H.A.P. Blom & J. Lygeros, editors: Stochastic hybrid systems: theory and safety critical applications, Lectures notes in control and information sciences (LNCIS) 337\. Springer, pp. 31–63.
  * [17] M.H.C. Everdij & H.A.P. Blom (2010): _Hybrid state Petri nets which have the analysis power of stochastic hybrid systems and the formal verification power of automata_.  In: P. Pawlewski, editor: Petri nets: Applications, chapter 12. InTech, pp. 227–252.
  * [18] M.H.C. Everdij, M.B. Klompstra, H.A.P. Blom & B. Klein Obbink (2006): _Compositional specification of a multi-agent system by stochastically and dynamically coloured Petri nets_.  In: H.A.P. Blom & J. Lygeros, editors: Stochastic hybrid systems: theory and safety critical applications, Lectures notes in control and information sciences (LNCIS) 337\. Springer, pp. 325–350.
  * [19] HYBRIDGE (2002).  _EC project description_.  http://hosted.nlr.nl/public/hosted-sites/hybridge.
  * [20] iFly (2007).  _EC project description_.  http://iFly.nlr.nl.
  * [21] J. Krystul (2006): _Modelling of stochastic hybrid systems with applications to accident risk assessment_.  Ph.D. thesis, University of Twente, The Netherlands.
  * [22] J. Krystul, H.A.P. Blom & A. Bagchi (2007): _Stochastic hybrid systems_ , chapter 2: Stochastic differential equations on hybrid state spaces, pp. 15–45.  Number 24 in Control engineering series. Taylor and Francis / CRC Press.
  * [23] RTCA (1995): _Final report of RTCA Task Force 3; Free Flight implementation_.  Final, RTCA Inc., Washington DC.
  * [24] A.J. Van der Schaft (2004): _Equivalence of dynamical systems by bisimulation_.  IEEE transactions on automatic control 49(12), pp. 2160–2172.


# Parametric Amplification of a Quantum Pulse

Offek Tziperman<EMAIL_ADDRESS>Equal contributors. Technion –
Israel Institute of Technology, 32000 Haifa, Israel  Victor Rueskov
Christiansen<EMAIL_ADDRESS>Equal contributors. Department of Physics
and Astronomy, Aarhus University, Ny Munkegade 120, DK-8000 Aarhus C, Denmark
Ido Kaminer  Technion – Israel Institute of Technology, 32000 Haifa, Israel
Klaus Mølmer<EMAIL_ADDRESS>Niels Bohr Institute, University of
Copenhagen, Blegdamsvej 17, 2100 Copenhagen, Denmark

###### Abstract

Creating and manipulating quantum states of light requires nonlinear
interactions, but while nonlinear optics is inherently multi-mode, quantum
optical analyses are often done with single-mode approximations. We present a
multi-mode theory for the transformation of a quantum pulse by Hamiltonians
that are quadratic in the field creation and annihilation operators. Our
theory describes nonlinear processes, such as parametric amplification and
squeezing, as well as all linear processes, such as dispersion and beam
splitting. We show that a single input pulse feeds only two distinct output
modes and, for certain quantum states, just one. Our theory provides the
quantum states in the output modes, which are crucial for the application of
pulses in quantum optics and quantum information.

Introduction— Quantum states of light are the key ingredients for some of the
most promising quantum technologies [1]. Light pulses can propagate between
stationary components in quantum networks [2], and their applications include
quantum key distribution [3], sensing beyond the standard quantum limit [4, 5,
6], bosonic error correcting codes [7, 8], and measurement based quantum
computation [9]. However, the preparation and manipulation of quantum states
of light are more complicated for multi-mode traveling fields than for their
stationary counterparts, such as the field in a single-mode cavity. This is
because a continuum of frequency modes is available for the radiation [10].

The multi-mode nature of quantum light is of critical importance in any system
involving nonlinearities, as in this case, the transformation of a quantum
state will generally lead to a population of output pulse shapes correlated
with the photon number state content [9, 10]. The output field is thus
inherently multi-mode in nature, and the quantum properties of the state may
not be accessible to the desired quantum information processing task.

In this article, we deal with operations on traveling light pulses, governed
by Hamiltonians that are of second order in the creation and annihilation
operators. This class of Hamiltonians can describe beam splitters and
interferometers, dispersion, diffraction and polarization rotation, but also
parametric amplification, parametric down conversion and frequency conversion.

We are particularly interested in parametric amplification which is in the
single (cavity) mode case described by the Hamiltonian ($\hbar=1$),
$H=\frac{i\xi}{2}((a^{\dagger})^{2}-a^{2})$, and the time evolution:

$\displaystyle a(t)=\cosh(\xi t)a(0)+\sinh(\xi t)a^{\dagger}(0).$ (1)

This transformation is also referred to as squeezing, as it reduces the value
(and uncertainty) of one of the field quadratures while the other quadrature
is correspondingly amplified [11]. Both the squeezing and amplification
properties are widely used in quantum optics for quantum sensing [6], quantum
state tomography [12, 13] and for signal amplification and read-out of
superconducting qubits [14, 15] or nanophotonics [16]. In its single mode
version, parametric amplification has been proposed as part of the gate set
for continuous variable quantum computing [17] and for creation and
manipulation of Gottesman-Kitaev-Preskill (GKP) states for error correction
[18, 19]. As there are now proposals and attempts to employ these schemes with
traveling wave packet modes, it is pertinent to assess the multi-mode
performance of the squeezing devices and ask to what extent Eq. (1) applies to
the transformation of the quantum state of a traveling pulse of radiation.

Figure 1: Parametric amplification of quantum pulses. (a) A non-linear medium
in a cavity is driven by an arbitrarily shaped classical pump (blue) such that
a single quantum pulse (red), transmitted through the cavity, is subject to
parametric amplification. The cavity output field populates several spatio-
temporal modes with components of the amplified pulse and with squeezed vacuum
contributions. (b) When the input pulse is prepared in a Schrödinger cat
state, a single photon state, and a vacuum state (shown in the left panels),
the mostly occupied output mode is in an approximate, amplified (squeezed)
version of these input states (shown in the right panels).

Here, we answer this question by presenting a general theory for the time
evolution of quantum pulses of radiation subject to optical components that
actuate quadratic Hamiltonians, such as dispersive and parametrically
amplifying elements. Due to the multi-mode character of the problem, the
transformation of quantum pulses by quadratic Hamiltonians is fundamentally
different from situations where a cavity restricts the dynamics to a single
mode. Still, even when subject to general, time-dependent Hamiltonians that
are quadratic in creation and annihilation operators, a single-mode input
quantum state transforms into an output retained in at most two spatiotemporal
modes. We can thus calculate the state occupying these, or any other desired
output wave packet mode, in which we recover a squeezed version of part of the
initial state mixed with various amounts of squeezed vacuum.

Figure 1 shows examples of parametric amplification of incident wave packet
states. The lower panels in Figure 1 show Wigner functions of the input states
and the states of the most populated output mode. While our theory takes the
multi-mode character of traveling fields explicitly into account, it shows
that it is still possible to obtain an approximate squeezing transformation
from a single input to a single output mode.

Quantum state transformation by a parametric amplifier—Our main goal is to
determine the output of the amplifier given any arbitrary quantum state
occupying an input wave packet $u(\omega)$. We treat here the case of guided
propagation of a single transversal field mode, while further spatial and
polarization degrees of freedom are readily included, see Sec. III (SMIII)) in
[20]. We define:
$a_{u,\textup{in}}^{\dagger}=\int_{-\infty}^{\infty}u(\omega)a^{\dagger}(\omega)d\omega,$
where $a(\omega)$ are bosonic operators with the commutation relation
$\left[a(\omega),a^{\dagger}(\omega^{\prime})\right]=\delta(\omega-\omega^{\prime})$.

We consider a general quadratic Hamiltonian in the creation and annihilation
operators [21]:

$\displaystyle\begin{split}H=&\int\int d\omega
d\omega^{\prime}K(\omega,\omega^{\prime})a^{\dagger}(\omega)a(\omega^{\prime})\\\
+&\int\int d\omega
d\omega^{\prime}J(\omega,\omega^{\prime})a^{\dagger}(\omega)a^{\dagger}(\omega^{\prime})+\text{h.c.},\end{split}$
(2)

where $K(\omega,\omega^{\prime})$ alone leads to a transformation of the input
wave packet, while $J(\omega,\omega^{\prime})$ is the amplitude for
simultaneously creating a pair of photons in the output field at frequency
$\omega$ and $\omega^{\prime}$, and h.c. denotes the Hermitian conjugate.

The equations of motion can be solved to give the transformation of the
operators in the Heisenberg picture [22]:

$\displaystyle\begin{split}a_{\textup{out}}(\omega)=&\int
d\omega^{\prime}F(\omega,\omega^{\prime})a_{\textup{in}}(\omega^{\prime})\\\
+&\int
d\omega^{\prime}G^{\ast}(\omega,\omega^{\prime})a_{\textup{in}}^{\dagger}(\omega^{\prime}),\end{split}$
(3)

where $F(\omega,\omega^{\prime})$ and $G(\omega,\omega^{\prime})$ are uniquely
determined by the functions $J(\omega,\omega^{\prime})$ and
$K(\omega,\omega^{\prime})$ [21], and we shall give examples in later sections
of the article.

The ideal transformation for parametric amplification of a single mode quantum
pulse would amount to a single output mode squeezed as in Eq. (1). However,
since the amplification applies to a continuum of field modes and is
frequency-dependent, the output field may occupy mode functions that will
generally differ from the input mode. We therefore move our focus to the
calculation of the quantum state content of any given output mode with
annihilation operator
$a_{v,\textup{out}}=\int_{-\infty}^{\infty}v^{\ast}(\omega^{\prime})a_{\textup{out}}(\omega^{\prime})d\omega^{\prime}$.
First, we shall analyse the mode content of the output to find the optimal
candidate $v$ function containing the transformed input quantum state, and in
the following section, we shall calculate its actual quantum state content.

Modes at the output of the amplifier— To characterize the output field, we
evaluate and decompose the first order coherence function
$g_{1}(\omega_{1},\omega_{2})=\braket{a^{\dagger}(\omega_{1})a(\omega_{2})}$
[10]. Consider the situation where the input occupies a single wave packet,
$u(\omega)$, with a quantum state given by a density matrix $\rho_{u}$, and
vacuum in all modes orthogonal to $u(\omega)$. We can then readily replace
$\braket{a^{\dagger}_{\textup{in}}(\omega^{\prime})a_{\textup{in}}(\omega^{\prime\prime})}$
by
$\braket{a^{\dagger}_{u}a_{u}}u^{*}(\omega^{\prime})u(\omega^{\prime\prime})$
and apply similar expressions for the other field-field correlation functions.
With Eq. (3) we find

$\displaystyle\begin{split}g_{1}(\omega_{1},\omega_{2})&=\braket{a_{u}^{\dagger}a_{u}}\int
d\omega^{\prime}F^{*}(\omega_{1},\omega^{\prime})u^{*}(\omega^{\prime})\int
d\omega^{\prime\prime}F(\omega_{2},\omega^{\prime\prime})u(\omega^{\prime\prime})+\braket{a_{u}^{\dagger}a_{u}^{\dagger}}\int
d\omega^{\prime}F^{*}(\omega_{1},\omega^{\prime})u^{*}(\omega^{\prime})\int
d\omega^{\prime\prime}G^{*}(\omega_{2},\omega^{\prime\prime})u^{*}(\omega^{\prime\prime})\\\
&+\braket{a_{u}a_{u}}\int
d\omega^{\prime}G(\omega_{1},\omega^{\prime})u(\omega^{\prime})\int
F(\omega_{2},\omega^{\prime\prime})u(\omega^{\prime\prime})d\omega^{\prime\prime}+\braket{a_{u}^{\dagger}a_{u}}\int
d\omega^{\prime}G(\omega_{1},\omega^{\prime})u(\omega^{\prime})\int
d\omega^{\prime\prime}G^{*}(\omega_{2},\omega^{\prime\prime})u^{*}(\omega^{\prime\prime})\\\
&+\int
d\omega^{\prime}G(\omega_{1},\omega^{\prime})G^{*}(\omega_{2},\omega^{\prime}).\end{split}$
(4)

We observe that the coherence function is governed by both the input pulse’s
shape and quantum state. The first four terms contain correlation functions of
the input field (for example $\braket{a_{u}a_{u}}$) while the last term is the
squeezed vacuum output of the amplifier in the absence of any input pulse.

By inspection of Eq. (4), it turns out that the input quantum pulse gives rise
to population of only two output modes, while a potentially infinite number of
modes are populated by squeezed vacuum,

$\displaystyle\begin{split}g_{1}(\omega_{1},\omega_{2})=&\underbrace{\sum^{2}_{i=1}n_{i}v_{i}^{\ast}(\omega_{1})v_{i}(\omega_{2})}_{\textup{dependent
on input quantum state}}+\\\
&\underbrace{\sum_{i=1}^{\infty}m_{i}w_{i}^{\ast}(\omega_{1})w_{i}(\omega_{2})}_{\textup{independent
of input quantum state}}.\end{split}$ (5)

This is demonstrated in detail in Sec. SMI [20].

In the special case of a coherent state input
($\rho_{u}=\ket{\alpha}\bra{\alpha}$) the operator expectation values factor,
and the above expansion leads to only a single output mode seeded by the input
field. For a Fock state input ($\rho_{u}=\ket{n}\bra{n}$) the output field
occupies two modes, seeded by the input pulse. In fact, the condition for
which the output occupies only a single mode is [20]

$\displaystyle\braket{a_{u}^{\dagger}a_{u}}=|\braket{a_{u}a_{u}}|.$ (6)

We note that these findings have intriguing and practical consequences for
quantum state tomography using parametric amplification [12, 16], where both
input fed output modes must be considered, as well as the noise input from
squeezed vacuum modes.

The quantum state in an output mode— While the operation of the amplifier is
entirely described by the multi-mode Bogoliubov transformation (3), this does
not provide an immediate description of the transformation of a given quantum
state input. We are faced, in fact, with an instance of the boson sampling
problem [23], and the number state representation of the output state is an
unwieldy expression in terms of matrix permanents and the density matrix of
the input state. The problem simplifies, however, when we restrict our
interest to the quantum state of the output of the amplifier in any single
mode $v(\omega)$. Our theory applies to any output mode of interest, but a
good choice for $v(\omega)$ is the most occupied single mode function at the
output, say, $v_{1}(\omega)$ in Eq.(5). Assuming any output mode function
$v(\omega)$, and by applying (3), we can formally write

$\displaystyle a_{v,\textup{out}}=\int
v^{\ast}(\omega)a_{\textup{out}}(\omega)d\omega=\zeta a_{f,\textup{in}}+\xi
a_{g,\textup{in}}^{\dagger}$ (7)

where we have defined two new input mode functions through $f^{*}(\omega)=\int
v^{*}(\omega^{\prime})F(\omega^{\prime},\omega)d\omega^{\prime}/\zeta$ and
$g(\omega)=\int
v^{*}(\omega^{\prime})G^{*}(\omega^{\prime},\omega)d\omega^{\prime}/\xi$, and
$\zeta,\xi$ ensure their normalization (we omit the subscript
$\cdot_{\textup{in}}$ in the following).

Notice a few key properties of this transformation: First, in the purely
dispersive case, $\xi=0,\;\zeta=1$ and as a result, the annihilation operator
of the output mode function $v(\omega)$ represents the exact same quantum
state content as the one that occupied the input mode function $f(\omega)$
(which is not necessarily $u(\omega)$). This reflects that the quantum states
are unchanged while traveling through a linearly dispersive element, but their
mode functions change shape.

Second, in the amplifier case, where $\xi\neq 0$, the functions $f,g$ are
normalized but generally not orthogonal, and their relationship with the
populated input mode $u(\omega)$ is not yet specified. To find the output
quantum state in the mode $v(\omega)$, we must decompose the transformation
into one that refers specifically to the input mode $u(\omega)$ and vacuum
modes orthogonal to $u(\omega)$. Defining $\braket{f,g}=\int
f^{\ast}(\omega)g(\omega)d\omega$, we decompose the modes into parallel and
orthogonal components of the input mode $u$,

$\displaystyle\begin{split}a_{v,\textup{out}}&=\zeta\braket{f,u}a_{u}+\xi\braket{u,g}a_{u}^{\dagger}\\\
&+\zeta\sqrt{1-|\braket{f,u}|^{2}}\braket{h,k}a_{k}+\xi\sqrt{1-|\braket{u,g}|^{2}}a_{k}^{\dagger}\\\
&+\zeta\sqrt{1-|\braket{f,u}|^{2}}\sqrt{1-|\braket{k,h}|^{2}}a_{s}.\end{split}$
(8)

Explicit expressions for the mode functions for $k$, $h$, and $s$ are given in
Sec. SMII in [20].

Eq. (8) is a main result of our analysis. It shows how a single output mode
captures a squeezed version of the potentially interesting state occupying the
input pulse (first line). However, it is mixed with a squeezed vacuum
component (second line) and a vacuum component (third line). In the absence of
the terms in the second and third line, we recover unitary single-mode
squeezing as in Eq. (1), while the presence of these terms in Eq. (8) will in
general contribute added quantum noise.

The transformation is effectively described by the Bloch-Messiah reduction
[24, 25], which states that a general multi-mode Bogoliubov transformation can
be separated in 3 steps: first, a linear beam splitter transformation among
the modes, then a sequence of single mode squeezing operations, and finally
another beam splitter transformation. We thus find
$\rho_{v\textup{,out}}=\textup{Tr}_{k,s}(U\rho_{u}U^{\dagger})$ with

$\displaystyle
U=U_{u,s}(\theta_{3},\phi_{3})U_{u,k}(\theta_{2},\phi_{2})S_{u}(r_{1})S_{k}(r_{2})U_{u,k}(\theta_{1},\phi_{1}),$
(9)

where the beam-splitter and squeezing transformations and parameters are
specified in Sec. SMII in [20]. Note that these transformations act in a
particularly simple manner on the states in the Wigner function
representation, where they amount to linear transformations on the field
quadrature variables [26], given directly by the coefficients in Eq. (8). By a
straightforward extension of this procedure, one can also find the joint
quantum state of the pair of output modes $v_{1},v_{2}$ fed by the input mode,
and, e.g., study their mutual entanglement properties, see Sec. SMII in [20]
and the code repository111Link to the repository with Python code to find the
output quantum state for a given transformation
https://github.com/offektziperman/Amplifying-a-quantum-pulse. In the
following, we shall restrict the analysis to the most occupied single mode
$v_{1}$ for $v_{\textup{out}}$, and leave the study of other candidate output
modes and the two-mode output for later investigation.

Amplification by an Optical Parametric Oscillator— To demonstrate applications
of our theory, we consider as a first example a degenerate Optical Parametric
Oscillator (OPO) cavity with the cavity Hamiltonian

$\displaystyle H_{\textup{sys}}=\Delta
a^{\dagger}a+\frac{i\xi(t)}{2}\left[a^{\dagger 2}-a^{2}\right],$ (10)

where $\Delta$ is the detuning of the cavity and $\xi(t)$ is the time-
dependent parametric gain due to a pulsed drive. We assume an input quantum
state occupying a single mode of temporal shape $u(t)$, coupled with the
coherent amplitude $\sqrt{\gamma}$ to the cavity mode (the cavity loss rate is
$\gamma$). Under the assumption of a Markovian coupling, we can relate the
output field to the cavity field at time $t$ via the input-output relation
222This input-output relation assumes unidirectional propagation of the field
which may be ensured by a ring cavity, a circulator, or the phase matching in
the OPA gain material in the absence of a cavity.
$a_{\textup{out}}(t)=a_{\textup{in}}(t)+\sqrt{\gamma}a(t)$ [29], and solving
the Heisenberg equation of motion for the intra-cavity field, we obtain the
explicit Bogoliubov transformation of the field operators in time domain in
the SMV [20]. This yields directly the time domain correlation function
$g_{1}(t_{1},t_{2})$, and the minimal basis of temporal modes, that
efficiently describe the parametrically amplified input quantum state.

Figure 2: Temporal modes emitted by a parametric amplifier. (a) Vacuum input.
Occupation of most populated modes as a function of pump duration for Gaussian
pumps with constant pulse area $(\int\xi(t)dt,\;\xi_{0}=1)$. A short pump
pulse excites the cavity mode into a squeezed state that leaks into a single-
mode travelling pulse. A longer pump pulse overlaps the emission process, and
more modes become occupied. The mode shapes for the limiting cases of very
wide and narrow pumps are plotted as insets. (b) Occupation of the modes fed
by a Gaussian pump $(\xi_{0}=1.5)$ and different initial quantum states in the
input pulse. The occupation of the modes seeded by an input pulse of duration
$\tau=1/\gamma$ (two modes for Fock states, one mode for coherent and
Schrödinger cat states) are plotted as a function of the pump width, assuming
a constant pump pulse area. We assume that the pump and quantum pulses peak at
the same time.

We analyse first the output of the OPO with only vacuum seed and a time-
dependent Gaussian pump
$\xi(t)=\frac{\xi_{0}}{\sqrt{2\pi\sigma_{\xi}^{2}}}\exp({-\frac{(t-t_{0})^{2}}{2\sigma_{\xi}^{2}}})$.
Figure 2(a) shows the occupation of the modes of the output field. For a short
pulsed pump, ($\gamma\sigma_{\xi}\ll 1$) only a single mode of the output
field is excited (squeezed), as can be seen by the large occupation of the
most occupied mode in the left-hand side of Figure 2(a). The mode occupied
represents the exponential decay (shown in the inset) of the abruptly excited
cavity mode. For wider pumps, the output field populates many modes (shown in
the right-hand inset of Figure 2(a)). In the limiting case of an infinite
duration constant pump, the photons possess strong frequency correlations,
$\omega_{1}+\omega_{2}=\omega_{\textup{pump}}$, to conserve energy. In this
limit, the output field explores a continuum of modes and
$g^{(1)}(\omega_{1},\omega_{2})\propto
S(\omega_{1})\delta(\omega_{\textup{pump}}-\omega_{1}-\omega_{2})$ where
$S(\omega_{1})$ is the output spectrum [30].

Next, we consider the parametric amplification of a quantum input pulse in the
mode $u(t)\propto\exp\left(-\frac{(t-t_{0})^{2}}{2\tau^{2}}\right)$ in Figure
2(b). Notice that the coherent state and the Schrödinger cat input feed only
into one output mode as they obey the condition in Eq. (6). In contrast, Fock
states lead to a population of two modes, as shown by the red and blue pairs
of curves in Figure 2(b). The efficiency of the gain is largest for the case
of a short pump. In this case, the emission by the cavity is negligible during
the gain process, and the accumulated cavity field stimulates the strongest
gain.

Figure 3: Comparison of three parametric amplifiers acting on a single photon
input state. Upper panels: Mean photon number in the dominant output mode
$n_{1}$. Lower panels: occupation ratio $n_{1}/(n_{1}+n_{2})$ as defined in Eq
(4). (a-b) For the OPO, these quantities are plotted as a function of the
delay time between the pump and probe and the pump temporal width in units of
the cavity decay rate $\gamma$. (c-d) For the OPA, with a spatially extended
non-linear medium and no cavity, they are shown as a function of the pump-seed
detuning and the pump spectral width in units of the spectral width of the
input quantum pulse $\sigma_{u}$. (e-f) For the TWPA, they are shown as a
function of the pump detuning in units of the cavity decay rate $\gamma$ and
the total amount of gain for 3000 OPOs.

Comparison of three parametric amplifiers— We recall that the output of the
OPO contains both the transformed input quantum pulse(s) and components of the
multi-mode squeezed vacuum. This results from our general analysis and applies
to any setup with quadratic Hamiltonians. We shall consider exemplary systems
and analyse the output field modes that are seeded by a single photon input
pulse $u(t)\propto\exp\left(-\frac{t^{2}}{2\sigma_{u}^{2}}\right)$. The three
setups are illustrated above the upper panels in Figure 3: 1. the Optical
Parametric Oscillator (OPO) as explained above in Eq. (10), 2. The Optical
Parametric Amplifier (OPA), with a spatially extended non-linear medium and no
cavity, and 3. the Traveling Wave Parametric Amplifier (TWPA), represented
here as a concatenated sequence of OPOs. For all three systems, we find
regimes of close to ideal single-mode amplification, as shown by the regimes
of large occupation in the dominant mode $v_{1}$ (upper panels in Figure 3)
and the large ratio between the occupation in the dominant mode $n_{1}$ and
the sum of the two most populated modes $n_{1}+n_{2}$ (lower panels Figure 3).

For the amplification of a single photon Gaussian pulse by the OPO, a short
delay between the short pump and quantum state seed is optimal, for the OPA, a
large overlap in frequency domain between the pump and seed will provide high
gain and single mode operation, while for the TWPA, we find that resonant
cavities provide the highest gain and highest single mode operation ratio.
More information on the calculations presented in Figure 3 is given in the
SMV-VII.

Figure 4: Transformation of quantum states by an OPO. Wigner functions are
shown for the most occupied output temporal modes for different input states,
(a) the vacuum state $\ket{0}$, (b) the single photon Fock state $\ket{n=1}$
and (c) the Schrödinger cat state $\ket{\alpha}+\ket{-\alpha}$ with
$\alpha=2.5$. Assuming the same pump pulse area, short pulses lead to the
strongest squeezing. (d) The purity of the state in the mode $v_{1}$ is
plotted on the left. For each pump temporal width, we calculate the quantum
state in the mode $v_{1}$, then we find the single mode parametrically
amplified state to which it has the highest fidelity, we plot this fidelity
and the corresponding squeezing on (d) (middle,right). The squeezing is
quantified by the amplitude gain in the $p$ quadrature.

Next, we implement the theory to calculate the quantum state content of the
most occupied mode. Results of these calculations are shown in Figure 4 for
squeezing by an OPO. The output Wigner functions in the most occupied single
mode are plotted for different input states: vacuum (a), a single photon state
(b), and a cat state (c) (squeezed cat states can be mixed to create GKP
states [7, 18]). Results are shown for different pump widths, assuming the
optimal delay identified in Figure 3. We observe the effect of parametric
amplification for short pulses while, in agreement with Figure 2., for longer
pump pulses, the quantum states experience little or no amplification. When
the pump area is increased by a factor of 5 (rightmost column of Figure 4),
the gain is stronger also for the longer pulses. However, in this regime the
quantum states are polluted by the multi-mode squeezed vacuum generated by the
OPO, decohering the quantum states.

To make a quantitative comparison of the transformation of quantum pulses with
the ideal single-mode squeezing transformation, we calculate the quantum state
occupying the mode $v_{1}$, and in the leftmost panel, we plot the purity of
this state, which can be lower than unity, both due to entanglement with the
second output mode seeded by the input state $v_{2}$ and due to contamination
by the squeezed vacuum modes $w_{i}$. In the rightmost panel, we compare the
state in $v_{1}$ to the input state squeezed by the single mode transformation
of Eq. (1), optimizing over the amount of squeezing, parameterized as gain in
the p quadrature, which we plot in the central panel of Figure 4(d). Notably,
we find pure state fidelities of over 85% with more than a factor 3 of
parametric amplification on the $p$-quadrature, while close to unity fidelity
parametric amplification of a single mode pulse is only possible for much
smaller gains, except for the special case of vacuum input. We note that we
considered an input seed and pump in a Gaussian temporal mode and that
optimizing over these may increase the fidelity, purity and squeezing.

Conclusion— An optical component described by a multi-mode field Hamiltonian
that is quadratic in creation and annihilation operators causes linear
dispersion of wave packets combined with correlated creation and annihilation
of photon pairs over all modes. A multi-mode Bogoliubov transformation
accounts for the effect of the Hamiltonian on the input field operators. We
have shown that if all but a single mode of the input field are in the vacuum
state, despite the multi-mode character of the problem, only two output modes
will contain quantum states that depend on the non-trivial input state. These
output modes will generally also contain components of squeezed vacuum from
the parametric amplifier. While these results are at variance with the ideal
parametric amplification as a unitary operation on a single-mode quantum
field, our analysis can identify optimal parameter settings for reducing
additional noise and spreading over more than a single mode.

Our theory applies also to multiple input and output modes, and the joint
quantum state of two or more modes can be readily found. Optimization of the
fidelity of the parametric amplification process and the more general
properties of the output quantum state of the multi-mode Bogoliubov
transformation constitute a promising topic for further study.

V. Rueskov Christiansen acknowledges support from the Danish National Research
Foundation through the Center of Excellence for Complex Quantum Systems (Grant
agreement No. DNRF156) and discussions with Daniel Holleufer. K. Mølmer
acknowledges support from the Carlsberg Foundation through the “Semper Ardens”
Research Project QCooL.

## References

  * O’Brien  _et al._ [2009] J. L. O’Brien, A. Furusawa, and J. Vučković, Photonic quantum technologies, Nature Photonics 3, 687 (2009).
  * Kimble [2008] H. J. Kimble, The quantum internet, Nature 453, 1023 (2008).
  * Korzh  _et al._ [2015] B. Korzh, C. C. W. Lim, R. Houlmann, N. Gisin, M. J. Li, D. Nolan, B. Sanguinetti, R. Thew, and H. Zbinden, Provably secure and practical quantum key distribution over 307 km of optical fibre, Nature Photonics 9, 163 (2015).
  * Caves [1982] C. M. Caves, Quantum limits on noise in linear amplifiers, Phys. Rev. D 26, 1817 (1982).
  * Caves [2020] C. M. Caves, Reframing su(1,1) interferometry, Advanced Quantum Technologies 3, 1900138 (2020).
  * Lawrie  _et al._ [2019] B. J. Lawrie, P. D. Lett, A. M. Marino, and R. C. Pooser, Quantum sensing with squeezed light, ACS Photonics 10.1021/acsphotonics.9b00250 (2019).
  * Gottesman  _et al._ [2001] D. Gottesman, A. Kitaev, and J. Preskill, Encoding a qubit in an oscillator, Phys. Rev. A 64, 012310 (2001).
  * Sivak  _et al._ [2023] V. V. Sivak, A. Eickbusch, B. Royer, S. Singh, I. Tsioutsios, S. Ganjam, A. Miano, B. L. Brock, A. Z. Ding, L. Frunzio, S. M. Girvin, R. J. Schoelkopf, and M. H. Devoret, Real-time quantum error correction beyond break-even, Nature 616, 50 (2023).
  * Briegel  _et al._ [2009] H. J. Briegel, D. E. Browne, W. Dür, R. Raussendorf, and M. Van den Nest, Measurement-based quantum computation, Nature Physics 5, 19 (2009).
  * Fabre and Treps [2020] C. Fabre and N. Treps, Modes and states in quantum optics, Rev. Mod. Phys. 92, 035005 (2020).
  * Caves [1981] C. M. Caves, Quantum-mechanical noise in an interferometer, Phys. Rev. D 23, 1693 (1981).
  * Kalash and Chekhova [2023] M. Kalash and M. V. Chekhova, Wigner function tomography via optical parametric amplification, Optica 10, 1142 (2023).
  * Li  _et al._ [2017] T. Li, B. E. Anderson, T. Horrom, B. L. Schmittberger, K. M. Jones, and P. D. Lett, Improved measurement of two-mode quantum correlations using a phase-sensitive amplifier, Opt. Express 25, 21301 (2017).
  * Macklin  _et al._ [2015] C. Macklin, K. O’Brien, D. Hover, M. E. Schwartz, V. Bolkhovsky, X. Zhang, W. D. Oliver, and I. Siddiqi, A near-quantum-limited josephson traveling-wave parametric amplifier, Science 350, 307 (2015).
  * Qiu  _et al._ [2023] J. Y. Qiu, A. Grimsmo, K. Peng, B. Kannan, B. Lienhard, Y. Sung, P. Krantz, V. Bolkhovsky, G. Calusine, D. Kim, A. Melville, B. M. Niedzielski, J. Yoder, M. E. Schwartz, T. P. Orlando, I. Siddiqi, S. Gustavsson, K. P. O’Brien, and W. D. Oliver, Broadband squeezed microwaves and amplification with a josephson travelling-wave parametric amplifier, Nature Physics 19, 706 (2023).
  * Nehra  _et al._ [2022] R. Nehra, R. Sekine, L. Ledezma, Q. Guo, R. M. Gray, A. Roy, and A. Marandi, Few-cycle vacuum squeezing in nanophotonics, Science 377, 1333 (2022).
  * Lloyd and Braunstein [1999] S. Lloyd and S. L. Braunstein, Quantum computation over continuous variables, Phys. Rev. Lett. 82, 1784 (1999).
  * Weigand and Terhal [2018] D. J. Weigand and B. M. Terhal, Generating grid states from schrödinger-cat states without postselection, Phys. Rev. A 97, 022341 (2018).
  * Yanagimoto  _et al._ [2023] R. Yanagimoto, R. Nehra, R. Hamerly, E. Ng, A. Marandi, and H. Mabuchi, Quantum nondemolition measurements with optical parametric amplifiers for ultrafast universal quantum information processing, PRX Quantum 4, 010333 (2023).
  * [20] See supplemental material at (url to be inserted by publisher) for derivations of central results .
  * Christ  _et al._ [2013] A. Christ, B. Brecht, W. Mauerer, and C. Silberhorn, Theory of quantum frequency conversion and type-ii parametric down-conversion in the high-gain regime, New Journal of Physics 15, 053038 (2013).
  * Bogolyubov [1958] N. N. Bogolyubov, A New method in the theory of superconductivity. I, Sov. Phys. JETP 7, 41 (1958).
  * Aaronson and Arkhipov [2011] S. Aaronson and A. Arkhipov, The computational complexity of linear optics, in  _Proceedings of the Forty-Third Annual ACM Symposium on Theory of Computing_, STOC ’11 (Association for Computing Machinery, New York, NY, USA, 2011) p. 333–342.
  * Braunstein [2005] S. L. Braunstein, Squeezing as an irreducible resource, Phys. Rev. A 71, 055801 (2005).
  * Cariolaro and Pierobon [2016] G. Cariolaro and G. Pierobon, Reexamination of bloch-messiah reduction, Phys. Rev. A 93, 062115 (2016).
  * Ekert and Knight [1989] A. K. Ekert and P. L. Knight, The wigner function of two-mode squeezed states; free and dissipative evolution, in  _Coherence and Quantum Optics VI_ , edited by J. H. Eberly, L. Mandel, and E. Wolf (Springer US, Boston, MA, 1989) pp. 255–259.
  * Note [1] Link to the repository with Python code to find the output quantum state for a given transformation https://github.com/offektziperman/Amplifying-a-quantum-pulse.
  * Note [2] This input-output relation assumes unidirectional propagation of the field which may be ensured by a ring cavity, a circulator, or the phase matching in the OPA gain material in the absence of a cavity.
  * Gardiner and Collett [1985] C. W. Gardiner and M. J. Collett, Input and output in damped quantum systems: Quantum stochastic differential equations and the master equation, Phys. Rev. A 31, 3761 (1985).
  * Collett and Gardiner [1984] M. J. Collett and C. W. Gardiner, Squeezing of intracavity and traveling-wave light fields produced in parametric amplification, Phys. Rev. A 30, 1386 (1984).