# A Markov jump process associated with the matrix-exponential distribution

Oscar Peralta111The University of Adelaide, School of Mathematical Sciences,
SA 5005, Australia<EMAIL_ADDRESS>
###### Abstract

Let $f$ be the density function associated to a matrix-exponential
distribution of parameters $(\bm{\alpha},T,\bm{s})$. By exponentially tilting
$f$, we find a probabilistic interpretation which generalises the one
associated to phase-type distributions. More specifically, we show that for
any sufficiently large $\lambda\geq 0$, the function
$x\mapsto\left(\int_{0}^{\infty}e^{-\lambda
s}f(s)\mathrm{d}s\right)^{-1}e^{-\lambda x}f(x)$ can be described in terms of
a Markov jump process whose generator is tied to $T$. Finally, we show how to
revert the exponential tilting in order to assign a probabilistic
interpretation to $f$ itself.

##  1 Introduction

A phase-type distribution corresponds to the law of $Y:=\inf\\{t\geq
0:J_{t}=\star\\}$ where $\\{J_{t}\\}_{t\geq 0}$ is a Markov jump process with
state space $\\{1,\dots,p\\}\cup\\{\star\\}$, with $\\{1,\dots,p\\}$ assumed
to be transient states and $\\{\star\\}$ absorbing. If $\\{J_{t}\\}_{t\geq 0}$
has a block-partitioned initial distribution $(\bm{\pi},0)$ and intensity
matrix given by

$\left[\begin{matrix}A&\bm{b}\\\
\bm{0}&0\end{matrix}\right]\quad\mbox{with}\quad\bm{b}=-A\bm{1},$ (1.1)

where $\bm{0}$ represents a $p$-dimensional row vector of $0$’s and $\bm{1}$ a
$p$-dimensional column vector of $1$’s, then we say that the phase-type
distribution is of parameters $(\bm{\pi},A)$. Via simple probabilistic
arguments, it can be shown that the density function of a phase-type
distribution of parameters $(\bm{\pi},A)$ is of the form

$g(x)=\bm{\pi}e^{Ax}\bm{b},\quad x\geq 0.$ (1.2)

Indeed, the vector $\bm{\pi}e^{Ax}$ yields the probabilities of
$\\{J_{t}\\}_{t\geq 0}$ being in some state $\\{1,\dots,p\\}$ at time $x$, and
$\bm{b}$ corresponds to the intensity vector of an absorption happening
immediately after. Phase-type distributions were first introduced in [12] with
the aim of constructing a robust and tractable class of distributions on
$\mathds{R}_{+}$ to be used in econometric problems. A more comprehensive
study of phase-type distributions was carried on by Neuts [13, 14], whose work
popularized their use in more general stochastic models.

On the other hand, a _matrix-exponential distribution_ of dimension $p\geq 1$
is an absolutely continuous distribution on $(0,\infty)$ whose density
function can be written as

$f(x)=\bm{\alpha}e^{Tx}\bm{s},\quad x\geq 0,$ (1.3)

where $\bm{\alpha}=(\alpha_{1},\dots,\alpha_{p})$ is a $p$-dimensional row
vector, $T=\\{t_{ij}\\}_{i,j\in\\{1,\dots,p\\}}$ is a $(p\times
p)$-dimensional square matrix, and $\bm{s}=(s_{1},\dots,s_{p})^{\intercal}$ is
a $p$-dimensional column vector, all with complex entries. If the dimension
need not be specified, we refer to such a distribution simply as _matrix-
exponential_. It follows from (1.2) and (1.3) that the class of phase-type
distributions is a subset of those that are matrix-exponential, with the
inclusion being strict (see [15] for details on the latter).

Matrix-exponential distributions were first studied in [7, 8] through the
concept of complex-valued transition probabilities. More precisely, it was
shown that certain systems with complex-valued elements can be formally
studied by analytical means without assigning a specific physical
interpretation to their components. While their method provided mathematical
rigour to systems “driven” by complex-valued intensity matrices, it failed to
provide a physical meaning to each individual component, as opposed to the
case of Markov jump processes with genuine intensity matrices. Later on, it
was proved in [15, 5] that matrix-exponential distributions have an
interpretation in terms of a Markov process with continuous state space, as
opposed to the finite-state-space one that phase-type distributions enjoy.
Even after the previous physical interpretations of matrix-exponential
distributions were discovered, properties of this class of distributions are
still not as well understood as they are for its phase-type counterpart. One
of the main reasons for this is that continuous state space processes are more
difficult to handle, so that studying matrix-exponential distributions by
physical means requires a more sophisticated framework. For example, this is
the case in [2, 4, 3] where the theory of piecewise-deterministic Markov
processes is used to study models with matrix-exponential components. Thus,
having a finite-state system interpretation for matrix-exponential
distributions available may potentially lead to the discovery of new
properties, as it has traditionally been the case for phase-type
distributions.

In this paper we give a physical interpretation to _each element_ of the
parameters $(\bm{\alpha},T,\bm{s})$ satisfying:

  1. A1.

The elements of $\bm{\alpha}$, $T$ and $\bm{s}$ are real,

  2. A2.

The dominant eigenvalue of $T$, denoted $\sigma_{0}$, is real and strictly
negative.

Since it can be shown that for a given matrix-exponential density of the form
(1.3) the parameters $(\bm{\alpha},T,\bm{s})$ can be chosen is such a way that
A1 and A2 hold (see [1]), the interpretation that we develop essentially
completes the picture laid out in [7, 8]. Our method, inspired by the recent
work in [16], provides a transparent interpretation of
$(\bm{\alpha},T,\bm{s})$ in terms of a finite-state Markov jump process. To do
so, we employ the technique known as _exponential tilting_ , which means that
we focus on the density proportional to $e^{-\lambda\cdot}f(\cdot)$ for large
enough $\lambda>0$. After we perform this transformation, we construct a
Markov jump process on a finite state space formed by two groups: the
_original states_ and the _anti-states_ , the latter being a copy of the
former. Heuristically, jumps within the set of original states or within the
set of anti-states occur according to the off-diagonal nonnegative “jump
intensities” of $T$, while jumps between the original and the anti-states
occur according to the negative “jump intensities” of $T$. Eventual absorption
or termination happens, and each realization “carries” a positive or negative
sign depending only on its initial and final state. Our main contribution is
showing that this mechanism yields the exponentially-tilted matrix-exponential
distribution, and by reverting the exponential tilting, providing some
probabilistic insight of the original matrix-exponential distribution as well.

The structure of the paper is as follows. In Section 2 we provide a brief
exposition on exponential tilting and how it affects the representation of a
matrix-exponential distribution. In Section 3 we present our main results,
Theorem 3.2 and Corollary 3.4, where we give a precise interpretation of an
exponentially-tilted matrix-exponential density in terms of a Markov jump
process. Finally, in Section 4 we provide methods to recover formulae and
probabilistic interpretations for any matrix-exponential distribution based on
the results on their exponentially-tilted version.

##  2 Preliminaries

Exponential tilting, also known as the _Escher transform_ , is a technique
which transforms any probability density function $f$ with support on
$[0,\infty)$ into a new probability density function $f_{\lambda}$ defined by

$f_{\lambda}(x)=\frac{e^{-\lambda x}f(x)}{\int_{0}^{\infty}e^{-\lambda
r}f(r)\mathrm{d}r},\quad x\geq 0,$

where $\lambda\geq 0$ is the _tilting rate_. The use of exponential tilting
goes back at least to [10], where it was used to build upon Crámer’s classical
actuarial models [9]. Later on, the exponential tilting method played a
prominent role in the theory of option pricing [11].

The exponentially-tilted version of a matrix-exponential distribution has a
simple form which happens to be matrix-exponential itself. To see this, notice
that if $f$ is of the form (1.3), then

$\displaystyle\int_{0}^{\infty}e^{-\lambda r}f(r)\mathrm{d}r$
$\displaystyle=\int_{0}^{\infty}e^{-\lambda
r}(\bm{\alpha}e^{Tr}\bm{s})\mathrm{d}r=\bm{\alpha}(\lambda I-T)^{-1}\bm{s},$

where we used the fact that $T-\lambda I$ has eigenvalues with strictly
negative real parts and thus $e^{(T-\lambda I)r}$ vanishes as
$r\rightarrow\infty$. Thus,

$\displaystyle f_{\lambda}(x)=\frac{e^{-\lambda
x}(\bm{\alpha}e^{Tx}\bm{s})}{\bm{\alpha}(\lambda
I-T)^{-1}\bm{s}}=\left(\frac{\bm{\alpha}}{\bm{\alpha}(\lambda
I-T)^{-1}\bm{s}}\right)e^{(T-\lambda I)x}\bm{s},\quad x\geq 0,$ (2.1)

implying that $f_{\lambda}$ corresponds to the density function of a matrix-
exponential distribution of parameters
$\left(\tfrac{\bm{\alpha}}{\bm{\alpha}(\lambda I-T)^{-1}\bm{s}},T-\lambda
I,\bm{s}\right)$.

Recall that the parameters $(\bm{\alpha},T,\bm{s})$ need not have a
probabilistic meaning in terms of a finite state space Markov chain, as
opposed to the parameters associated to phase-type distributions. For
instance, the parameters

$\bm{\alpha}=\left(1,0,0\right),\quad T=\left[\begin{array}[]{rrr}-1&-1&2/3\\\
1&-1&-2/3\\\
0&0&-1\end{array}\right],\quad\bm{s}=\left[\begin{array}[]{r}4/3\\\ 2/3\\\
1\end{array}\right]$ (2.2)

yield a valid matrix-exponential distribution whose density function is given
by $f(x)=\tfrac{2}{3}e^{-x}(1+\cos(x))$, and where the dominant eigenvalue of
$T$ is $-1$ (see [6, Example 4.5.21] for details). In the following section we
show how to assign a probabilistic meaning to the exponentially-tilted version
of (2.2), and more generally to those attaining properties A1 and A2, in terms
of a finite-state Markov jump process.

##  3 Main results

Let $(\bm{\alpha},T,\bm{s})$ be parameters associated to a $p$-dimensional
matrix-exponential distribution which attain properties A1 and A2. For $1\leq
i,j\leq p$ denote by $t_{ij}$ the $(i,j)$-entry of $T$, and denote by $s_{i}$
the $i$-th entry of $\bm{s}$. For $\ell\in\\{+,-\\}$, define the $(p\times
p)$-dimensional matrix $T^{\ell}=\\{t_{ij}^{\ell}\\}_{1\leq i,j\leq p}$ and
the $p$-dimensional column vector
$\bm{s}^{\ell}=(s^{\ell}_{1},\dots,s^{\ell}_{p})^{\intercal}$ where

$\displaystyle t_{ij}^{\pm}$ $\displaystyle=\max\\{0,\pm
t_{ij}\\}\quad\forall\,i\neq j,$ $\displaystyle t_{ii}^{\pm}$
$\displaystyle=\pm\min\\{0,\pm t_{ii}\\}\quad\forall\,i,\mbox{ and}$
$\displaystyle s_{i}^{\pm}$ $\displaystyle=\max\\{0,\pm
s_{i}\\}\quad\forall\,i.$

It follows that $T^{+}$ has nonnegative off-diagonal elements and nonpositive
diagonal elements, $T^{-}$ is a nonnegative matrix, $\bm{s}^{+}$ and
$\bm{s}^{-}$ are nonnegative column vectors, $T=T^{+}-T^{-}$ and
$\bm{s}=\bm{s}^{+}-\bm{s}^{-}$. Now, let

$\lambda_{0}=\min\left\\{r\geq
0:\sum_{j=1}^{p}s^{+}_{ij}+s^{-}_{ij}+s^{+}_{i}+s^{-}_{i}-\lambda\leq 0\mbox{
for all }1\leq i\leq p\right\\}.$

For some fixed $\lambda\geq\lambda_{0}$, consider a (possibly) terminating
Markov jump process $\\{\varphi_{t}^{\lambda}\\}_{t\geq 0}$ driven by the
block-partitioned subintensity matrix

$\displaystyle G=\left[\begin{array}[]{cccc}T^{+}-\lambda
I&T^{-}&\bm{s}^{+}&\bm{s}^{-}\\\ T^{-}&T^{+}-\lambda
I&\bm{s}^{-}&\bm{s}^{+}\\\ \bm{0}&\bm{0}&0&0\\\
\bm{0}&\bm{0}&0&0\end{array}\right]$ (3.5)

evolving on the state space
$\mathcal{E}=\mathcal{E}^{o}\cup\mathcal{E}^{a}\cup\\{\Delta^{o}\\}\cup\\{\Delta^{a}\\}$
where $\mathcal{E}^{o}:=\\{1^{o},2^{o},\dots,p^{o}\\}$ and
$\mathcal{E}^{a}=\\{1^{a},2^{a},\dots,p^{a}\\}$. The state space $\mathcal{E}$
may be thought as the union of two sets: a collection of _original states_
$\mathcal{E}^{o}\cup\\{\Delta^{o}\\}$ and a collection of _anti-states_
$\mathcal{E}^{a}\cup\\{\Delta^{a}\\}$, where both $\Delta^{o}$ and
$\Delta^{a}$ are absorbing. In the case $\lambda>\lambda_{0}$, the process
$\\{\varphi_{t}^{\lambda}\\}_{t\geq 0}$ alternates between sojourn times in
$\mathcal{E}^{o}$ and $\mathcal{E}^{a}$ up until one of the following happens:
(a) get absorbed into $\Delta^{o}$, (b) get absorbed into $\Delta^{a}$ or (c)
undergo termination due to the defect of (3.5). If $\lambda=\lambda_{0}$, the
states $\mathcal{E}^{o}\cup\mathcal{E}^{a}$ may or may not be transient, their
status depending on the values of $T$.

In Theorem 3.1 we establish a link between the absoprtion probabilities of
$\\{\varphi_{t}^{\lambda}\\}_{t\geq 0}$ and the vector $e^{(T-\lambda
I)x}\bm{s}$ appearing in the exponentially-tilted matrix-exponential density
(2.1). More specifically, we express each element in $e^{(T-\lambda
I)x}\bm{s}$ as the sum of some positive density functions and some negative
density function, where the positive density is associated to an absorption of
$\\{\varphi_{t}^{\lambda}\\}_{t\geq 0}$ to $\Delta^{o}$, while the negative
density function corresponds to an absorption of
$\\{\varphi_{t}^{\lambda}\\}_{t\geq 0}$ to $\Delta^{a}$. To shorten notation,
from now on we denote by $\mathds{P}_{j}$ ($\mathds{E}_{j}$) ,
$j\in\mathcal{E}$, the probability measure (expectation) associated to
$\\{\varphi^{\lambda}_{t}\\}_{t\geq 0}$ conditional on the event
$\\{\varphi^{\lambda}_{0}=j\\}$.

######  Theorem 3.1.

Let $\lambda\geq\lambda_{0}$ be such that the states
$\mathcal{E}^{o}\cup\mathcal{E}^{a}$ are transient. Define

$\tau=\inf\\{x\geq
0:\varphi_{x}^{\lambda}\notin\mathcal{E}^{o}\cup\mathcal{E}^{a}\\}.$ (3.6)

Then, for $i\in\\{1,\dots,p\\}$ and $x\geq 0$,

$\displaystyle\big{(}\bm{e}_{i}^{\intercal}e^{(T-\lambda
I)x}\bm{s}\big{)}\mathrm{d}x$
$\displaystyle=\mathds{E}_{i^{o}}\left[\mathds{1}\\{\tau\in[x,x+\mathrm{d}x]\\}\beta(\varphi_{\tau}^{\lambda})\right]$
(3.7)
$\displaystyle=(\bm{e}_{i}^{\intercal},\bm{0})\exp\left(\left[\begin{matrix}T^{+}-\lambda
I&T^{-}\\\ T^{-}&T^{+}-\lambda
I\end{matrix}\right]x\right)\left[\begin{matrix}\bm{s}\\\
-\bm{s}\end{matrix}\right]\mathrm{d}x,$ (3.8)

where $\bm{e}_{i}$ denotes the column vector with $1$ on its $i$-th entry and
$0$ elsewhere, and
$\beta(j):=\mathds{1}\\{j=\Delta^{o}\\}-\mathds{1}\\{j=\Delta^{a}\\}$.
Moreover,

$\displaystyle\big{(}-\bm{e}_{i}^{\intercal}e^{(T-\lambda
I)x}\bm{s}\big{)}\mathrm{d}x$
$\displaystyle=\mathds{E}_{i^{a}}\left[\mathds{1}\\{\tau\in[x,x+\mathrm{d}x]\\}\beta(\varphi_{\tau}^{\lambda})\right]$
(3.9)
$\displaystyle=(\bm{0},\bm{e}_{i}^{\intercal})\exp\left(\left[\begin{matrix}T^{+}-\lambda
I&T^{-}\\\ T^{-}&T^{+}-\lambda
I\end{matrix}\right]x\right)\left[\begin{matrix}\bm{s}\\\
-\bm{s}\end{matrix}\right]\mathrm{d}x.$ (3.10)

###### Proof.

The block structure of (3.5) implies that

$\displaystyle\mathds{P}_{i^{o}}(\tau\in[x,x+\mathrm{d}x],\varphi_{\tau}^{\lambda}=\Delta^{o})$
$\displaystyle=(\bm{e}_{i}^{\intercal},\bm{0})\exp\left(\left[\begin{matrix}T^{+}-\lambda
I&T^{-}\\\ T^{-}&T^{+}-\lambda
I\end{matrix}\right]x\right)\left[\begin{matrix}\bm{s}^{+}\\\
\bm{s}^{-}\end{matrix}\right],$
$\displaystyle\mathds{P}_{i^{o}}(\tau\in[x,x+\mathrm{d}x],\varphi_{\tau}^{\lambda}=\Delta^{a})$
$\displaystyle=(\bm{e}_{i}^{\intercal},\bm{0})\exp\left(\left[\begin{matrix}T^{+}-\lambda
I&T^{-}\\\ T^{-}&T^{+}-\lambda
I\end{matrix}\right]x\right)\left[\begin{matrix}\bm{s}^{-}\\\
\bm{s}^{+}\end{matrix}\right],$

therefore, the r.h.s. of (3.7) is equal to (3.8). We prove that (3.7) holds
next.

Define the collection of $(p\times p)$-dimensional matrices
$\\{\Phi_{oa}(x)\\}_{x\geq 0}$, $\\{\Phi_{ao}(x)\\}_{x\geq 0}$,
$\\{\Phi_{oo}(x)\\}_{x\geq 0}$ and $\\{\Phi_{oa}(x)\\}_{x\geq 0}$ by

$\displaystyle(\Phi_{oa}(x))_{ij}$
$\displaystyle=\mathds{P}_{i^{o}}(\tau>x,\varphi_{x}^{\lambda}=j^{a}),\qquad(\Phi_{ao}(x))_{ij}=\mathds{P}_{i^{a}}(\tau>x,\varphi_{x}^{\lambda}=j^{o}),$
$\displaystyle(\Phi_{oo}(x))_{ij}$
$\displaystyle=\mathds{P}_{i^{o}}(\tau>x,\varphi_{x}^{\lambda}=j^{o}),\qquad(\Phi_{aa}(x))_{ij}=\mathds{P}_{i^{a}}(\tau>x,\varphi_{x}^{\lambda}=j^{a}),$

for all $i,j\in\\{1,\dots,p\\}$. By the symmetry of the subintensity matrix
$G$ it is clear that for all $x\geq 0$, $\Phi_{oa}(x)=\Phi_{ao}(x)$ and
$\Phi_{oo}(x)=\Phi_{aa}(x)$, even if their probabilistic interpretations
differ. For all $x\geq 0$, let $\Phi_{o}(x):=\Phi_{oo}(x)-\Phi_{oa}(x)$ and
$\Phi_{a}(x):=\Phi_{aa}(x)-\Phi_{ao}(x)$. Define

$\gamma=\inf\\{r\geq 0:\varphi_{r}^{\lambda}\notin\mathcal{E}^{o}\\}.$

Then, for $i,j\in\\{1,\dots,p\\}$

$\displaystyle\bm{e}_{i}^{\intercal}\Phi_{o}(x)\bm{e}_{j}$
$\displaystyle=\mathds{P}_{i^{o}}(\tau>x,\varphi_{x}^{\lambda}=j^{o})-\mathds{P}_{i^{o}}(\tau>x,\varphi_{x}^{\lambda}=j^{a})$
$\displaystyle=\left\\{\mathds{P}_{i^{o}}(\gamma>x,\varphi_{x}^{\lambda}=j^{o})+\int_{0}^{x}\mathds{P}_{i^{o}}(\gamma\in[r,r+\mathrm{d}r],\tau>x,\varphi_{x}^{\lambda}=j^{o})\right\\}$
$\displaystyle\quad-\left\\{\mathds{P}_{i^{o}}(\gamma>x,\varphi_{x}^{\lambda}=j^{a})+\int_{0}^{x}\mathds{P}_{i^{o}}(\gamma\in[r,r+\mathrm{d}r],\tau>x,\varphi_{x}^{\lambda}=j^{a})\right\\}$
$\displaystyle=\mathds{P}_{i^{o}}(\gamma>x,\varphi_{x}^{\lambda}=j^{o})$
$\displaystyle\quad+\int_{0}^{x}\sum_{k=1}^{p}\mathds{P}_{i^{o}}(\gamma\in[r,r+\mathrm{d}r],\varphi_{\gamma}^{\lambda}=k^{a})\mathds{P}_{k^{a}}(\tau>x-r,\varphi_{x-r}^{\lambda}=j^{o})$
$\displaystyle\quad-\int_{0}^{x}\sum_{k=1}^{p}\mathds{P}_{i^{o}}(\gamma\in[r,r+\mathrm{d}r],\varphi_{\gamma}^{\lambda}=k^{a})\mathds{P}_{k^{a}}(\tau>x-r,\varphi_{x-r}^{\lambda}=j^{a}),$
(3.11)

where in the last equality we used that
$\\{\gamma>x,\varphi_{x}^{\lambda}=j^{a}\\}=\varnothing$ and the Markov
property of $\\{\varphi_{x}^{\lambda}\\}_{x\geq 0}$. Note that all the
elements in (3.11) correspond to transition probabilities or intensities that
can be expressed in matricial form as follows:

$\displaystyle\mathds{P}_{i^{o}}(\gamma>x,\varphi_{x}^{\lambda}=j^{o})$
$\displaystyle=\bm{e}_{i}^{\intercal}e^{(T^{+}-\lambda I)x}\bm{e}_{j},$
$\displaystyle\mathds{P}_{i^{o}}(\gamma\in[r,r+\mathrm{d}r],\varphi_{\gamma}^{\lambda}=k^{a})$
$\displaystyle=\bm{e}_{i}^{\intercal}e^{(T^{+}-\lambda
I)r}T^{-}\bm{e}_{k}\mathrm{d}r,$
$\displaystyle\mathds{P}_{k^{a}}(\tau>x-r,\varphi_{x-r}^{\lambda}=j^{o})$
$\displaystyle=\bm{e}_{k}^{\intercal}\Phi_{ao}(x-r)\bm{e}_{j},$
$\displaystyle\mathds{P}_{k^{a}}(\tau>x-r,\varphi_{x-r}^{\lambda}=j^{a})$
$\displaystyle=\bm{e}_{k}^{\intercal}\Phi_{aa}(x-r)\bm{e}_{j}.$

Substituting these expressions into (3.11) gives

$\displaystyle\bm{e}_{i}^{\intercal}\Phi_{o}(x)\bm{e}_{j}$
$\displaystyle=\bm{e}_{i}^{\intercal}e^{(T^{+}-\lambda
I)x}\bm{e}_{j}+\int_{0}^{x}\sum_{k=1}^{p}\left(\bm{e}_{i}^{\intercal}e^{(T^{+}-\lambda
I)r}T^{-}\bm{e}_{k}\right)\left(\bm{e}_{k}^{\intercal}\Phi_{ao}(x-r)\bm{e}_{j}\right)\mathrm{d}r$
$\displaystyle\quad-\int_{0}^{x}\sum_{k=1}^{p}\left(\bm{e}_{i}^{\intercal}e^{(T^{+}-\lambda
I)r}T^{-}\bm{e}_{k}\right)\left(\bm{e}_{k}^{\intercal}\Phi_{aa}(x-r)\bm{e}_{j}\right)\mathrm{d}r$
$\displaystyle=\bm{e}_{i}^{\intercal}\left(e^{(T^{+}-\lambda
I)x}+\int_{0}^{x}e^{(T^{+}-\lambda
I)r}T^{-}\left[\Phi_{ao}(x-r)-\Phi_{aa}(x-r)\right]\mathrm{d}r\right)\bm{e}_{j}$
$\displaystyle=\bm{e}_{i}^{\intercal}\left(e^{(T^{+}-\lambda
I)x}+\int_{0}^{x}e^{(T^{+}-\lambda
I)r}T^{-}\left[\Phi_{oa}(x-r)-\Phi_{oo}(x-r)\right]\mathrm{d}r\right)\bm{e}_{j}$
$\displaystyle=\bm{e}_{i}^{\intercal}\left(e^{(T^{+}-\lambda
I)x}+\int_{0}^{x}e^{(T^{+}-\lambda
I)r}(-T^{-})\Phi_{o}(x-r)\mathrm{d}r\right)\bm{e}_{j},$

so that $\\{\Phi_{o}(x)\\}_{x\geq 0}$ is the bounded solution to the matrix-
integral equation

$\Phi_{o}(x)=e^{(T^{+}-\lambda I)x}+\int_{0}^{x}e^{(T^{+}-\lambda
I)r}(-T^{-})\Phi_{o}(x-r)\mathrm{d}r.$

By [3, Theorem 3.10],

$\Phi_{o}(x)=e^{[(T^{+}-\lambda I)+(-T^{-})]x}=e^{(T-\lambda I)x}.$

The Markov property implies that

$\displaystyle\mathds{P}_{i^{o}}(\tau\in[x,x+\mathrm{d}x],\varphi_{\tau}^{\lambda}=\Delta^{o})$
$\displaystyle\quad=\sum_{k=1}^{p}\mathds{P}_{i^{o}}(\tau>x,\varphi_{x}^{\lambda}=k^{o})\mathds{P}_{k^{o}}(\tau\in[x,x+\mathrm{d}x],\varphi_{\tau}^{\lambda}=\Delta^{o})$
$\displaystyle\quad\quad+\sum_{k=1}^{p}\mathds{P}_{i^{o}}(\tau>x,\varphi_{x}^{\lambda}=k^{a})\mathds{P}_{k^{a}}(\tau\in[x,x+\mathrm{d}x],\varphi_{\tau}^{\lambda}=\Delta^{o})$
$\displaystyle\quad=\sum_{k=1}^{p}(\bm{e}_{i}^{\intercal}\Phi_{oo}(x)\bm{e}_{k})(\bm{e}_{k}^{\intercal}\bm{s}^{+})\mathrm{d}x+\sum_{k=1}^{p}(\bm{e}_{i}^{\intercal}\Phi_{oa}(x)\bm{e}_{k})(\bm{e}_{k}^{\intercal}\bm{s}^{-}\mathrm{d}x)$
$\displaystyle\quad=\bm{e}_{i}^{\intercal}(\Phi_{oo}(x)\bm{s}^{+}+\Phi_{oa}(r)\bm{s}^{-})\mathrm{d}x.$

Similarly,

$\mathds{P}_{i^{o}}(\tau\in[x,x+\mathrm{d}x],\varphi_{\tau}^{\lambda}=\Delta^{a})=\bm{e}_{i}^{\intercal}(\Phi_{oa}(x)\bm{s}^{+}+\Phi_{oo}(x)\bm{s}^{-})\mathrm{d}x.$

Thus,

$\displaystyle\mathds{P}_{i^{o}}(\tau\in[x,x+\mathrm{d}x],\varphi_{\tau}^{\lambda}=\Delta^{o})-\mathds{P}_{i^{o}}(\tau\in[x,x+\mathrm{d}x],\varphi_{\tau}^{\lambda}=\Delta^{a})$
$\displaystyle\quad=\bm{e}_{i}^{\intercal}([\Phi_{oo}(x)\bm{s}^{+}+\Phi_{oa}(x)\bm{s}^{-}]-[\Phi_{oa}(x)\bm{s}^{+}+\Phi_{oo}(x)\bm{s}^{-}])\mathrm{d}x$
$\displaystyle\quad=\bm{e}_{i}^{\intercal}([\Phi_{oo}(x)-\Phi_{oa}(x)]\bm{s}^{+}-[\Phi_{oo}(x)-\Phi_{oa}(x)]\bm{s}^{-})\mathrm{d}x$
$\displaystyle\quad=\bm{e}_{i}^{\intercal}\Phi_{o}(x)\bm{s}\mathrm{d}x=\bm{e}_{i}^{\intercal}e^{(T-\lambda
I)x}\bm{s}\mathrm{d}x,$

so that (3.7) holds. Analogous arguments follow for (3.9) and (3.10), which
completes the proof. ∎

Heuristically, Equations (3.7) and (3.9) imply that initiating
$\\{\varphi^{\lambda}_{t}\\}_{t\geq 0}$ in the anti-state $i^{a}$ has the
opposite effect, in terms of sign, to initiating in the original state
$i^{o}$. In the following we exploit this fact to provide a probabilistic
interpretation not only for the elements of $e^{(T-\lambda I)x}\bm{s}$, but
for the exponentially-tilted matrix-exponential density
$\bm{\alpha}e^{(T-\lambda I)x}\bm{s}$.

Define $w^{+}$ and $w^{-}$ by

$w^{\pm}=\sum_{i=1}^{p}\max\\{0,\pm\alpha_{i}\\},$

and define $\bm{\alpha}^{+}=(\alpha_{1}^{+},\dots,\alpha_{p}^{+})$ and
$\bm{\alpha}^{-}=(\alpha_{1}^{-},\dots,\alpha_{p}^{-})$ by

$\displaystyle\alpha^{\pm}_{i}=\left\\{\begin{array}[]{ccc}\frac{1}{w^{\pm}}\max\\{0,\pm\alpha_{i}\\}&\mbox{if}&w^{\pm}>0\\\
0&\mbox{if}&w^{\pm}=0.\end{array}\right.$

If $w^{\pm}>0$ then $\bm{\alpha}^{\pm}$ is a probability vector and in
general,

$\bm{\alpha}=w^{+}\bm{\alpha}^{+}-w^{-}\bm{\alpha}^{-}.$ (3.12)

In some sense, $(w^{+}+w^{-})^{-1}\bm{\alpha}$ can be thought as a mixture of
the probability vectors $\bm{\alpha}^{+}$ and $\bm{\alpha}^{-}$, with the
latter contributing “negative mass”. Fortunately, this “negative mass” in the
context of $\bm{\alpha}e^{(T-\lambda I)x}\bm{s}$ can be given a precise
probabilistic interpretation by means of anti-states as follows.

######  Theorem 3.2.

Let $f_{\lambda}(x)=(\bm{\alpha}(\lambda
I-T)^{-1}\bm{s})^{-1}\bm{\alpha}e^{(T-\lambda I)x}\bm{s}$, $x\geq 0$, be the
density of the exponentially-tilted matrix-exponential distribution of
parameters $(\bm{\alpha},T,\bm{s})$. Define the vectors

$\widehat{\bm{\alpha}}^{+}:=\tfrac{w^{+}}{w^{+}+w^{-}}\bm{\alpha}^{+}\quad\mbox{and}\quad\widehat{\bm{\alpha}}^{-}:=\tfrac{w^{-}}{w^{+}+w^{-}}\bm{\alpha}^{-},$

and suppose
$\varphi_{0}^{\lambda}\sim(\widehat{\bm{\alpha}}^{+},\widehat{\bm{\alpha}}^{-})$.
Then,

$\displaystyle f_{\lambda}(x)\mathrm{d}x$
$\displaystyle=\frac{w^{+}+w^{-}}{\bm{\alpha}(\lambda
I-T)^{-1}\bm{s}}\mathds{E}\left[\mathds{1}\\{\tau\in[x,x+\mathrm{d}x]\\}\beta(\varphi_{\tau}^{\lambda})\right]$
(3.13) $\displaystyle=\frac{w^{+}+w^{-}}{\bm{\alpha}(\lambda
I-T)^{-1}\bm{s}}\left((\widehat{\bm{\alpha}}^{+},\widehat{\bm{\alpha}}^{-})\exp\left(\left[\begin{matrix}T^{+}-\lambda
I&T^{-}\\\ T^{-}&T^{+}-\lambda
I\end{matrix}\right]x\right)\left[\begin{matrix}\bm{s}\\\
-\bm{s}\end{matrix}\right]\right)\mathrm{d}x,$ (3.14)

where $\tau$ and $\beta(\cdot)$ are defined as in Theorem 3.1.

###### Proof.

Equation (3.12) implies that

$\displaystyle f_{\lambda}(x)$ $\displaystyle=\frac{1}{\bm{\alpha}(\lambda
I-T)^{-1}\bm{s}}\left(\sum_{i=1}^{p}w^{+}\alpha^{+}_{i}\left(\bm{e}_{i}^{\intercal}e^{(T-\lambda
I)x}\bm{s}\right)+\sum_{i=1}^{p}w^{-}\alpha^{-}_{i}\left(-\bm{e}_{i}^{\intercal}e^{(T-\lambda
I)x}\bm{s}\right)\right)$
$\displaystyle=\frac{w^{+}+w^{-}}{\bm{\alpha}(\lambda
I-T)^{-1}\bm{s}}\left(\sum_{i=1}^{p}\tfrac{w^{+}}{w^{+}+w^{-}}\alpha^{+}_{i}\left(\bm{e}_{i}^{\intercal}e^{(T-\lambda
I)x}\bm{s}\right)+\sum_{i=1}^{p}\tfrac{w^{-}}{w^{+}+w^{-}}\alpha^{-}_{i}\left(-\bm{e}_{i}^{\intercal}e^{(T-\lambda
I)x}\bm{s}\right)\right).$ (3.15)

Equality (3.13) follows from (3.15), (3.7) and (3.9). Equality (3.14) follows
from (3.15), (3.8) and (3.10). ∎

######  Example 3.3.

Let $(\bm{\alpha},T,\bm{s})$ be the matrix-exponential parameters
corresponding to (2.2). As noted previously, these parameters by themselves
lack a probabilistic interpretation, so we apply Theorem 3.1 to construct one.
For such parameters we take the tilting parameter $\lambda:=\lambda_{0}=2$,
leading to the block-partitioned matrices

$\displaystyle\left[\begin{matrix}T^{+}-\lambda I&T^{-}\\\ T^{-}&T^{+}-\lambda
I\end{matrix}\right]=\left[\begin{array}[]{rrr|rrr}-3&0&2/3&0&1&0\\\
1&-3&0&0&0&2/3\\\ 0&0&-3&0&0&0\\\ \hline\cr 0&1&0&-3&0&2/3\\\
0&0&2/3&1&-3&0\\\ 0&0&0&0&0&-3\end{array}\right],$ (3.22)
$(\widehat{\bm{\alpha}}^{+},\widehat{\bm{\alpha}}^{-})=(1,0,0,0,0,0),\qquad\left[\begin{matrix}\bm{s}\\\
-\bm{s}\end{matrix}\right]=\left[\begin{array}[]{r}4/3\\\ 2/3\\\ 1\\\
\hline\cr-4/3\\\ -2/3\\\ -1\end{array}\right],$

and $w^{+}=1$, $w^{-}=0$ and $\bm{\alpha}(\lambda I-T)^{-1}\bm{s}=4$. We can
then verify that

$\displaystyle\frac{w^{+}+w^{-}}{\bm{\alpha}(\lambda
I-T)^{-1}\bm{s}}\left((\widehat{\bm{\alpha}}^{+},\widehat{\bm{\alpha}}^{-})\exp\left(\left[\begin{matrix}T^{+}-\lambda
I&T^{-}\\\ T^{-}&T^{+}-\lambda
I\end{matrix}\right]x\right)\left[\begin{matrix}\bm{s}\\\
-\bm{s}\end{matrix}\right]\right)$
$\displaystyle\qquad=\frac{1}{6}e^{-3x}(1+\cos(x))=\frac{e^{-2x}}{4}\left(\frac{2}{3}e^{-x}(1+\cos(x))\right),$

the latter corresponding to the exponentially-tilted matrix-exponential
density function $f(x)=\tfrac{2}{3}e^{-x}(1+\cos(x))$.∎

A probabilistic interpretation of $f_{\lambda}$ alternative to that of (3.13)
is the following.

######  Corollary 3.4.

Define $\bm{d}=(d_{1},\dots,d_{p})^{\intercal}:=-(T^{+}-\lambda
I)\bm{1}-T^{-}\bm{1}$ to be the termination intensities vector from
$\mathcal{E}^{o}$ or $\mathcal{E}^{a}$, and define
$\bm{q}^{\pm}=(q_{1}^{\pm},\dots,q_{p}^{\pm})^{\intercal}$ by

$q^{\pm}_{i}=\left\\{\begin{array}[]{ccc}\frac{s^{\pm}_{i}}{d_{i}}&\mbox{if}&d_{i}>0,\\\
0&\mbox{if}&d_{i}=0.\end{array}\right.$

Let $\bar{q}:\mathcal{E}^{o}\cup\mathcal{E}^{a}\mapsto\mathds{R}$ be defined
by

$\displaystyle\bar{q}(i^{o})=q^{+}_{i}-q^{-}_{i}\quad\mbox{ and
}\quad\bar{q}(i^{a})=q^{-}_{i}-q^{+}_{i}\quad\mbox{ for }\quad
i\in\\{1,\dots,p\\}.$

Then,

$\displaystyle f_{\lambda}(x)\mathrm{d}x$
$\displaystyle=\left(\frac{w^{+}+w^{-}}{\bm{\alpha}(\lambda
I-T)^{-1}\bm{s}}\right)\mathds{E}\left[\mathds{1}\\{\tau\in[x,x+\mathrm{d}x]\\}\bar{q}\big{(}\varphi^{\lambda}_{\tau^{-}}\big{)}\right]$
(3.23)

where $\\{\varphi^{\lambda}_{t}\\}_{t\geq 0}$ and $\tau$ are defined as in
Theorem 3.2.

###### Proof.

First, notice that the jump mechanism of $\\{\varphi^{\lambda}_{t}\\}_{t\geq
0}$ described in (3.5) imply that for $i\in\\{1,\dots,p\\}$,

$\displaystyle\mathds{P}(\varphi^{\lambda}_{\tau}=\Delta^{o}\mid\tau,\varphi^{\lambda}_{\tau^{-}}=i^{o})$
$\displaystyle=q_{i}^{+},$
$\displaystyle\mathds{P}(\varphi^{\lambda}_{\tau}=\Delta^{a}\mid\tau,\varphi^{\lambda}_{\tau^{-}}=i^{o})$
$\displaystyle=q_{i}^{-},$
$\displaystyle\mathds{P}(\varphi^{\lambda}_{\tau}=\Delta^{o}\mid\tau,\varphi^{\lambda}_{\tau^{-}}=i^{a})$
$\displaystyle=q_{i}^{-},$
$\displaystyle\mathds{P}(\varphi^{\lambda}_{\tau}=\Delta^{a}\mid\tau,\varphi^{\lambda}_{\tau^{-}}=i^{a})$
$\displaystyle=q_{i}^{+},$

which in turn implies that

$\displaystyle\mathds{E}\left[\beta(\varphi_{\tau}^{\lambda})\mid\tau,\varphi^{\lambda}_{\tau^{-}}\right]=\mathds{P}\left(\varphi_{\tau}^{\lambda}=\Delta^{o}\mid\tau,\varphi^{\lambda}_{\tau^{-}}\right)-\mathds{P}\left(\varphi_{\tau}^{\lambda}=\Delta^{a}\mid\tau,\varphi^{\lambda}_{\tau^{-}}\right)=\bar{q}(\varphi^{\lambda}_{\tau^{-}}).$

Consequently,

$\displaystyle\mathds{E}\left[\mathds{1}\\{\tau\in[x,x+\mathrm{d}x]\\}\beta(\varphi_{\tau}^{\lambda})\right]$
$\displaystyle=\mathds{E}\left[\mathds{E}\left[\mathds{1}\\{\tau\in[x,x+\mathrm{d}x]\\}\beta(\varphi_{\tau}^{\lambda})\mid\tau,\varphi^{\lambda}_{\tau^{-}}\right]\right]$
$\displaystyle=\mathds{E}\left[\mathds{1}\\{\tau\in[x,x+\mathrm{d}x]\\}\mathds{E}\left[\beta(\varphi_{\tau}^{\lambda})\mid\tau,\varphi^{\lambda}_{\tau^{-}}\right]\right]$
$\displaystyle=\mathds{E}\left[\mathds{1}\\{\tau\in[x,x+\mathrm{d}x]\\}\bar{q}\big{(}\varphi^{\lambda}_{\tau^{-}}\big{)}\right]$

and the result follows from (3.13). ∎

Though closely related, the interpretation provided by Corollary 3.4 is more
suitable than that of Theorem 3.2 for Monte Carlo applications. Indeed, a
realization of $\\{\varphi^{\lambda}_{t}\\}_{t\geq 0}$ may get absorbed in
$\Delta^{o}$, $\Delta^{a}$ or terminated. If termination is the case, such
realization contributes nothing to the term in the r.h.s. of (3.13). In
contrast, by observing the process until its exit time of
$\mathcal{E}^{o}\cup\mathcal{E}^{a}$ and disregard its landing point as in
Corollary 3.4, we make sure that each realization contributes towards the mass
in the r.h.s. of (3.23).

##  4 Recovering the untilted distribution

Once the exponentially-tilted density $f_{\lambda}$ of a matrix-exponential
distribution of parameters $(\bm{\alpha},T,\bm{s})$ has a tractable known
form, say as in (3.14), in principle it is straightforward to recover the
original untilted density $f$ by taking

$\displaystyle f(x)$ $\displaystyle=(\bm{\alpha}(\lambda I-T)\bm{s})e^{\lambda
x}f_{\lambda}(x)$
$\displaystyle=(w^{+}+w^{-})(\widehat{\bm{\alpha}}^{+},\widehat{\bm{\alpha}}^{-})\exp\left(\left[\begin{matrix}T^{+}&T^{-}\\\
T^{-}&T^{+}\end{matrix}\right]x\right)\left[\begin{matrix}\bm{s}\\\
-\bm{s}\end{matrix}\right],\qquad x\geq 0.$ (4.1)

While (4.1) is a legitimate matrix-exponential representation of $f$, it has
two drawbacks:

  1. 1.

The matrix $\left[\begin{smallmatrix}T^{+}&T^{-}\\\
T^{-}&T^{+}\end{smallmatrix}\right]$ may no longer be a subintensity matrix.

  2. 2.

The dominant eigenvalue of $\left[\begin{smallmatrix}T^{+}&T^{-}\\\
T^{-}&T^{+}\end{smallmatrix}\right]$ may be nonnegative.

The first item may impact the probabilistic interpretation of $f$, while the
second one may make integration of certain functions (with respect to the
density $f$) more difficult to handle. For instance, in the context of Example
3.3, the matrix

$\left[\begin{matrix}T^{+}&T^{-}\\\
T^{-}&T^{+}\end{matrix}\right]=\left[\begin{array}[]{rrr|rrr}-1&0&2/3&0&1&0\\\
1&-1&0&0&0&2/3\\\ 0&0&-1&0&0&0\\\ \hline\cr 0&1&0&-1&0&2/3\\\
0&0&2/3&1&-1&0\\\ 0&0&0&0&0&-1\end{array}\right]$

is not a subintensity matrix since some row sums are strictly positive, and it
has $0$ as its dominant eigenvalue. Having $0$ as an eigenvalue implies that
some entries of $\exp\left(\left(\begin{smallmatrix}T^{+}&T^{-}\\\
T^{-}&T^{+}\end{smallmatrix}\right)x\right)$ can potentially be of order
$e^{0\cdot x}=1$, meaning that the matrix-integral

$\int_{0}^{\infty}h(x)\exp\left(\left[\begin{smallmatrix}T^{+}&T^{-}\\\
T^{-}&T^{+}\end{smallmatrix}\right]x\right)\mathrm{d}x$ (4.2)

may only be well-defined for functions $h:\mathds{R}_{+}\mapsto\mathds{R}_{+}$
that decrease to $0$ fast enough. In comparison, $\exp(Tx)$ with $T$ as in
(2.2) has entries of at most $e^{\sigma_{0}x}=e^{-x}$, so that

$\int_{0}^{\infty}h(x)\exp\left(Tx\right)\mathrm{d}x$ (4.3)

is well-defined and finite for every function
$h:\mathds{R}_{+}\mapsto\mathds{R}_{+}$ of the order $e^{\sigma x}$ for any
$\sigma<1$. This apparent disagreement between the applicability of (4.2) and
(4.3) vanishes when we multiply
$\exp\left(\left[\begin{smallmatrix}T^{+}&T^{-}\\\
T^{-}&T^{+}\end{smallmatrix}\right]x\right)$ with
$\left[\begin{smallmatrix}\bm{s}\\\ -\bm{s}\end{smallmatrix}\right]$. Indeed,
in the context of Example 3.3 it can be verified that the elements of the
vector $\exp\left(\left[\begin{smallmatrix}T^{+}&T^{-}\\\
T^{-}&T^{+}\end{smallmatrix}\right]x\right)\left[\begin{smallmatrix}\bm{s}\\\
-\bm{s}\end{smallmatrix}\right]$ are at most of order $e^{-x}$, with the
higher order terms of $\exp\left(\left[\begin{smallmatrix}T^{+}&T^{-}\\\
T^{-}&T^{+}\end{smallmatrix}\right]x\right)$ cancelling each other when we
multiply the matrix-function by $\left[\begin{smallmatrix}\bm{s}\\\
-\bm{s}\end{smallmatrix}\right]$. In the general case, this “cancellation” of
higher order terms than those of $e^{\sigma_{0}x}$ occuring _after_ the matrix
multiplication can be directly deduced from Theorem 3.1.

In terms of expectations, (3.13) and (3.23) provide alternative ways to
recover properties of any matrix-exponential density $f$ of parameters
$(\bm{\alpha},T,\bm{s})$ in terms of the exponentially-tilted density
$f_{\lambda}$. Indeed, for any function
$h:\mathds{R}_{+}\mapsto\mathds{R}_{+}$ with
$\int_{0}^{\infty}h(x)f(x)\mathrm{d}x<\infty$, we have that

$\displaystyle\int_{0}^{\infty}h(x)f(x)\mathrm{d}x$
$\displaystyle=(\bm{\alpha}(\lambda I-T)\bm{s})\int_{0}^{\infty}h(x)e^{\lambda
x}f(x)\mathrm{d}x$ (4.4)
$\displaystyle=(w^{+}+w^{-})\mathds{E}\left[h(\tau)e^{\lambda\tau}\beta(\varphi^{\lambda}_{\tau})\right]$
(4.5)
$\displaystyle=(w^{+}+w^{-})\mathds{E}\left[h(\tau)e^{\lambda\tau}\bar{q}(\varphi^{\lambda}_{\tau^{-}})\right],$
(4.6)

where $\\{\varphi^{\lambda}_{t}\\}_{t\geq 0}$ and $\tau$ are as in Theorem
3.2. Notice that, as oppossed to the formula in (4.1), representations (4.5)
and (4.6) still have probabilistic interpretations in terms of the Markov jump
process $\\{\varphi^{\lambda}_{t}\\}_{t\geq 0}$.

## Acknowledgement

The author acknowledges the funding of the Australian Research Council
Discovery Project DP180103106.

## References

  * [1] S. Asmussen and M. Bladt.  Renewal theory and queueing algorithms for matrix-exponential distributions, pages 313–341.  Lecture notes in pure and applied mathematics. Marcel Dekker, 1996.
  * [2] S. Asmussen and M. Bladt.  Point processes with finite-dimensional conditional probabilities.  Stochastic Processes and their Applications, 82(1):127–142, 1999\.
  * [3] N. G. Bean, G. T. Nguyen, B. F. Nielsen, and O. Peralta.  RAP-modulated fluid processes: First passages and the stationary distribution.  arXiv preprint arXiv:2101.03242, 2021.
  * [4] N. G. Bean and B. F. Nielsen.  Quasi-birth-and-death processes with rational arrival process components.  Stochastic Models, 26(3):309–334, 2010.
  * [5] M. Bladt and M. F. Neuts.  Matrix-exponential distributions: calculus and interpretations via flows.  Stochastic Models, 19(1):113–124, 2003.
  * [6] M. Bladt and B. F. Nielsen.  Matrix-exponential distributions in applied probability, volume 81.  Springer, 2017.
  * [7] D. R. Cox.  A use of complex probabilities in the theory of stochastic processes.  Mathematical Proceedings of the Cambridge Philosophical Society, 51(02):313–319, 1955.
  * [8] D. R. Cox.  The analysis of non-Markovian stochastic processes by the inclusion of supplementary variables.  Mathematical Proceedings of the Cambridge Philosophical Society, 51(03):433, 1955.
  * [9] H. Cramér.  Collective risk theory: A survey of the theory from the point of view of the theory of stochastic processes.  In Collected Works II, pages 1028–1115. Springer Berlin Heidelberg, Berlin, Heidelberg, 2013.
  * [10] F. Escher.  On the probability function in the collective theory of risk.  Skand. Aktuarie Tidskr., 15:175–195, 1932.
  * [11] H. Gerber and E. Shiu.  Option pricing by esscher transforms.  Transactions of the Society of Actuaries, 46, 1994.
  * [12] A. Jensen.  A Distribution Model, Applicable to Economics.  Copenhagen, 1954.
  * [13] M. Neuts.  Probability distributions of phase type.  In Liber Amicorum Prof. Emeritus H. Florin, pages 173–206. Department of Mathematics, University of Louvain, Belgium, 1975.
  * [14] M. Neuts.  Matrix-geometric Solutions in Stochastic Models: An Algorithmic Approach.  Dover Publications, 1981.
  * [15] C. A. O’Cinneide.  Characterization of phase-type distributions.  Stochastic Models, 6(1):1–57, 1990.
  * [16] F. Völlering.  Markov process representation of semigroups whose generators include negative rates.  Electronic Communications in Probability, 25, 2020.


# Cavity-enhanced polarization rotation measurements for low-disturbance
probing of atoms

Chiara Mazzinghi1  Daniel Benedicto Orenes1  Pau Gomez1  Vito G. Lucivero1
Enes Aybar1  Stuti Gugnani1 and Morgan W. Mitchell1,2  1ICFO - Institut de
Ciencies Fotoniques, The Barcelona Institute of Science and Technology,
Castelldefels (Barcelona) 08860, Spain  2ICREA - Institució Catalana de
Recerca i Estudis Avançats, 08010 Barcelona, Spain<EMAIL_ADDRESS>
###### Abstract

We propose and demonstrate cavity-enhanced polarization-rotation measurement
as a means to detect magnetic effects in transparent media with greater
sensitivity at equal optical disturbance to the medium. Using the Jones
calculus, we compute the effective polarization rotation effect in a Fabry-
Perot cavity containing a magnetic medium, including losses due to enclosure
windows or other sources. The results show that when measuring polarization
rotation, collecting the transmitted light has advantages in simplicity and
linearity relative to collecting the reflected light. We demonstrate the
technique by measuring Faraday rotation in a 87Rb atomic ensemble in the
single-pass and cavity-enhanced geometries, and observe enhancement in good
agreement with the theoretical predictions. We also demonstrate shot-noise-
limited operation of the enhanced rotation scheme in the small-angle regime.

††journal: oe††articletype: Research Article

##  1 Introduction

Dispersive measurement of atomic variables using near-resonant light is both a
practical method of non-destructive probing [1, 2] and a versatile technique
for preparing non-classical states of matter [3, 4, 5]. When the atomic system
is contained within an optical resonator, both the measurement uncertainty and
the disturbance to the measured variable can be far below the intrinsic
quantum noise, leading to strong squeezing and entanglement [6, 7, 8]. To
date, most cavity-enhanced probing techniques have been applied to non-
magnetic transitions in atomic systems, of interest for atomic clocks [9, 10].
There has been less work with cavity-enhanced measurement of atomic
magnetization [11, 12, 13], of interest to magnetometers [14], gyroscopes
[15], and instruments to search for physics beyond the standard model [16, 17,
18].

Polarization rotation (PR) measurements, in which a probe beam’s linear
polarization rotates upon propagating through a medium, are widely applied to
detect magnetic effects in atomic media. These include Faraday rotation, in
which the magnetic field produces circular birefringence via the Zeeman effect
[19], and the so-called paramagnetic Faraday rotation, in which the spin
polarization or magnetization of the medium is responsible for the rotation
[20]. While it is to be expected that PR measurements can be enhanced by the
Purcell effect created by cavity resonance, the methods used to date,
including transmission [6, 7] and Pound-Drever-Hall (PDH) [8, 10] measurements
of the line shift of a single cavity mode, are not directly applicable to PR,
which involves two modes of different polarization. This intrinsically two-
mode character of PR probing presents novel challenges in the design of
cavity-enhancement (CE) methods.

In this work we propose and demonstrate a novel cavity-enhancement scheme to
boost the sensitivity of PR measurements. We first describe the atom-light
system using the Jones matrix formalism [21, 22] to obtain analytic
expressions for the cavity’s output field in transmission and in reflection.
We note various inconveniences of a PDH-like PR measurement in reflection, and
focus on the transmission geometry. We then report an experimental
implementation, in which a resonant optical cavity is built around a 87Rb
vapor cell. We compare the polarization rotation angle with and without the
cavity, and observe an enhancement factor in agreement with the theoretical
model. We measure the noise arising in this implementation, and conclude that
the probing system is shot-noise limited in a relevant regime for pulsed
probing of cold atomic ensembles [23, 24]. The new probing scheme will enable
cavity-enhanced quantum-non-demolition measurement [25] of magnetic degrees of
freedom in hot and cold atomic systems, creation of non-classical states of
matter [24, 26], and magnetic sensitivity beyond the standard quantum limit
[27].

Figure 1: Model for cavity-enhanced polarization-rotation measurements. A
quasi-monochromatic field $E_{\rm in}$ impinges upon an optical cavity formed
by mirrors CM1 and CM2, with intensity reflection (transmission) coefficients
$R_{1}(T_{1})$ and $R_{2}(T_{2})$, respectively. An atomic medium in the
cavity introduces a phase shift $\mp\phi_{\mathrm{at}}$ for the $\sigma_{\pm}$
polarizations on each pass (blue oval) between cavity mirrors. Optical windows
with intensity transmission coefficients $T_{G}$ scatter light out of the
cavity mode (red arrows), but do not produce back-reflection. Reflected or
transmitted output field $E_{\cal R}$ or $E_{\cal T}$, respectively, is
analyzed by measurement of the Stokes parameter $S_{2}$.

##  2 Model of cavity-enhanced polarization rotation

The components of the model are illustrated in Figure 1: an atomic medium is
placed within a Fabry-Perot optical resonator that, apart from possible
birefringent effects arising from the atomic medium itself, is polarization-
independent. Horizontally polarized input light, with electric field $E_{\rm
in}$, impinges upon the cavity and output fields $E_{\cal R}$ (reflected) and
$E_{\cal T}$ (transmitted), which are in general elliptically polarized,
containing information about the magnetic conditions within the atomic medium.
A balanced polarimeter (not shown in Figure 1) detects a Stokes component of
$E_{\cal R}$ or $E_{\cal T}$.

###  2.1 Jones Matrix description

We treat the probe light as quasi-monochromatic, and assume it is spatially
matched to the cavity’s fundamental TEM00 mode. At any given point in space
and in a given direction (forward or backward propagating), the field is
completely specified by the amplitudes of circular left (L) and circular right
(R) polarizations. Because the cavity per se is polarization-insensitive and
the atoms are not, it is more convenient to describe the field in terms of its
amplitudes for $\sigma_{+}$ and $\sigma_{-}$ polarization, i.e. the circular
polarizations with respect to the atomic quantization axis, chosen to be the
forward propagation direction.

We thus define Jones vectors
$E_{\mathrm{\alpha}}\equiv(E_{\mathrm{\alpha}}^{(+)},E_{\mathrm{\alpha}}^{(-)})^{T}$,
where $\alpha$ indicates a location and direction, and
$E_{\mathrm{\alpha}}^{(\pm)}$ indicates the electric field amplitude for the
$\sigma_{\pm}$ polarization at $\alpha$. For example, the input field,
immediately before reaching the input mirror and propagating in the forward
direction, is $E_{\mathrm{in}}=E_{0}\ (1/\sqrt{2},1/\sqrt{2})^{T}$, where
$E_{0}$ is the amplitude of the input field. Each optical transformation, e.g.
transmission through the medium, propagation by a distance, or reflection from
a mirror, is represented by a matrix that acts on the Jones vector. These
matrices are given in Table 1.

optical transformation | symbol | Jones matrix
---|---|---
free space propagation (round trip) | $P_{\mathrm{L}}$ | $e^{i\phi_{0}}\mathbbm{1}$
transmission through mirror $m$ | $T_{\mathrm{C}_{m}}$ | $\sqrt{T_{m}}\mathbbm{1}$
reflection from mirror $m$ | $R_{\mathrm{C}_{m}}$ | $\sqrt{R_{m}}\mathbbm{1}$
transmission through window | $T_{\mathrm{W}}$ | $\sqrt{T_{G}}\mathbbm{1}$
reflection from window | N/A | $0\cdot\mathbbm{1}$
transmission through medium | $F_{A}$ | $\left(\begin{array}[]{cc}e^{-i\phi_{\mathrm{at}}/2}&0\\\ 0&e^{i\phi_{\mathrm{at}}/2}\end{array}\right)$
polarization component | Stokes component | detection matrix
total power | $S_{0}$ | $\frac{1}{2}\mathbbm{1}$
linear | $S_{1}$ | $\frac{1}{2}\left(\begin{array}[]{cc}0&1\\\ 1&0\\\ \end{array}\right)$
diagonal | $S_{2}$ | $\frac{1}{2}\left(\begin{array}[]{cc}0&-i\\\ i&0\\\ \end{array}\right)$
Table 1: Jones matrices used to describe various steps of field propagation in
the model. $\mathbbm{1}$ indicates the $2\times 2$ identity matrix. $\phi_{0}$
is the phase shift acquired on a round trip of the cavity, including
polarization-independent phases acquired in propagation through the atomic
medium, cell windows, and upon reflection from mirrors. $R_{m}$ and $T_{m}$
are the intensity reflection and transmission coefficients of the $m$th
mirror. $T_{G}$ is the intensity transmission coefficient of the windows,
assumed equal. Reflection from a window is assumed to fall outside the TEM00
mode, and thus does not propagate further within the model.
$\phi_{\mathrm{at}}$ is the relative phase between $\sigma_{+}$ and
$\sigma_{-}$ polarizations upon passing through the atomic medium, equal to
half the geometrical polarization rotation angle. Because the PR is of
magnetic origin, the same matrix $F_{A}$ applies for either forward or
backward propagation through the medium.

It is convenient to use Stokes components [28, 22] to describe the
polarization at the detection stage. Considering a pulse of duration $\tau$
and effective beam cross-section $A_{\rm eff}$, the Stokes parameter
$S_{\beta}$, in photon-number units, is

${S}_{\beta}=\frac{c\tau A_{\rm eff}\epsilon_{0}}{4}E_{\rm
out}^{*}\hat{S}_{\beta}E_{\rm out},$ (1)

where $\hat{S}_{\beta}$ is the corresponding matrix, given in Table 1.
Considering that the input light is fully polarized in the horizontal
direction, i.e., $S_{1}^{\mathrm{in}}=S_{0}^{\mathrm{in}}$, the angle of
polarization rotation is
$\psi=\arcsin({S}_{2}^{\mathrm{out}}/S_{0}^{\mathrm{out}})$. The superscript
“in” refers to the forward-propagating field that enters the system, and “out”
refers to the output field, which could be in the forward or backward
direction, depending on where the detection is placed.

###  2.2 Single pass and cavity schemes

We consider three cases: single pass (i.e. without cavity), cavity in
reflection and in transmission. We will be interested in two figures of merit,
the rotation gain $\kappa\equiv\psi/\phi_{\mathrm{at}}$ and the system
efficiency $\eta_{\rm sys}\equiv S_{0}^{\mathrm{out}}/S_{0}^{\mathrm{in}}$.

###  2.3 Single pass PR measurement

The output field is

$E_{\mathrm{SP}}=T_{\rm W}F_{\rm A}T_{\rm W}E_{\rm
in}=\frac{E_{0}T_{G}}{\sqrt{2}}\begin{pmatrix}e^{-i\phi_{\mathrm{at}}/2}\\\
e^{i\phi_{\mathrm{at}}/2}\end{pmatrix}$ (2)

from which $\psi^{\mathrm{SP}}=\phi_{\mathrm{at}}$ and thus
$\kappa^{\mathrm{SP}}=1$. The system efficiency is $\eta_{\rm
sys}^{\mathrm{SP}}=T_{G}^{2}$.

###  2.4 Cavity enhanced PR measurement

We repeat the calculation of the output field as in previous section, but
include now the cavity mirrors. The input field is considered as in the
previous section, and the output field $E_{\cal T}$ ($E_{\cal R}$) from cavity
in transmission (in reflection) is:

$\displaystyle E_{\cal T}$ $\displaystyle=$ $\displaystyle T_{\rm C1}T_{\rm
C2}T_{\rm W}^{2}P_{\rm L}F_{\rm A}\sum_{p=0}^{\infty}\Bigl{(}R_{\rm C1}R_{\rm
C2}T_{\rm W}^{4}P_{\rm L}^{2}F_{\rm A}^{2}\Bigr{)}^{p}E_{\rm in}$ (3)
$\displaystyle E_{\cal R}$ $\displaystyle=$ $\displaystyle\biggl{[}-R_{\rm
C1}+T_{\rm C1}^{2}R_{\rm C2}T_{\rm W}^{4}P_{\rm L}^{2}F_{\rm
A}^{2}\sum_{p=0}^{\infty}\Bigl{(}R_{\rm C1}R_{\rm C2}T_{\rm W}^{4}P_{\rm
L}^{2}F_{\rm A}^{2}\Bigr{)}^{p}\biggr{]}E_{\rm in}.$ (4)

The geometric series can be evaluated analytically, to find

$\displaystyle E_{\mathcal{T}}^{(\pm)}$ $\displaystyle=$
$\displaystyle\frac{E_{0}}{\sqrt{2}}\biggl{(}\frac{\sqrt{T_{1}T_{2}}T_{G}\exp{[i\phi_{0}\pm
i\phi_{\mathrm{at}}/2]}}{1-T_{G}^{2}\sqrt{R_{2}R_{1}}\exp{[i\phi_{0}\pm
i\phi_{\mathrm{at}}]}}\biggr{)}$ (5) $\displaystyle E_{\mathcal{R}}^{(\pm)}$
$\displaystyle=$
$\displaystyle\frac{E_{0}}{\sqrt{2}}\biggl{(}-\sqrt{R_{1}}+\frac{\sqrt{R_{2}}T_{1}T_{G}^{2}\exp{[i\phi_{0}\pm
i\phi_{\mathrm{at}}]}}{1-T_{G}^{2}\sqrt{R_{2}R_{1}}\exp{[i\phi_{0}\pm
i\phi_{\mathrm{at}}]}}\biggr{)}$ (6)

For simplicity and since it is the highest-sensitivity scenario, in what
follows we consider only the on-resonance case, that is $\phi_{0}=0$. To
quantify the performance of the cavity in boosting the rotation signal, we
define

$\psi^{\mathrm{CE}}\equiv\arcsin({S_{2}^{\mathcal{T},\mathcal{R}}}/{S_{0}^{\mathcal{T},\mathcal{R}}})=\kappa_{\mathcal{R,T}}\psi^{\mathrm{SP}}$
(7)

where $\kappa_{\mathcal{R,T}}$ is the cavity enhancement factor, i.e. rotation
gain factor. The system efficiency is $\eta_{\rm sys}^{{\cal R},{\cal
T}}\equiv S_{0}^{\mathcal{T},\mathcal{R}}/S_{0}^{\mathrm{in}}$, in analogy
with the single-pass case.

###  2.5 Sensitivity and sensitivity enhancement at equal disturbance

The shot-noise limited sensitivity to rotation angle can be quantified using
the propagation of error formula. Writing $\hat{\phi}_{\mathrm{at}}$ for the
estimate of $\phi_{\mathrm{at}}$, the mean squared error (MSE) of
$\hat{\phi}_{\mathrm{at}}$ is

$\displaystyle\mathrm{MSE}(\hat{\phi}_{\mathrm{at}})$ $\displaystyle=$
$\displaystyle\left(\frac{\partial
S_{2}^{\mathrm{out}}}{\partial\phi_{\mathrm{at}}}\right)^{-2}\mathrm{var}(S_{2}^{\mathrm{out}})$
(8) $\displaystyle\approx$
$\displaystyle\left(\frac{S_{0}^{\mathrm{out}}\partial\psi}{\partial\phi_{\mathrm{at}}}\right)^{-2}\frac{S_{0}^{\mathrm{out}}}{2},$

where the approximation assumes $\psi\ll 1$ and that shot noise dominates over
electronic and technical noise.

Meanwhile, for linear-optical disturbance effects such as optical pumping, the
disturbance to the atoms is proportional to the spatially-averaged probe
intensity. Defining as $2S_{0}^{\rm int}$ the photon number seen by the atoms,
$S_{0}^{\rm int}$ can be expressed in terms of $S_{0}^{\mathrm{in}}$, which
through $\eta_{\rm sys}$ is related also to $S_{0}^{\mathrm{out}}$. In the SP
case $(S_{0}^{\rm int})_{\rm SP}=T_{G}(S_{0}^{\mathrm{in}})_{\rm SP}$. In the
CE case, the forward-propagating intracavity field $E_{\rm int}$ is

$\displaystyle E_{\rm int}$ $\displaystyle=T_{\rm C1}T_{\rm
W}\sum_{p=0}^{\infty}\Bigl{(}R_{\rm C1}R_{\rm C2}T_{\rm W}^{4}P_{\rm
L}^{2}F_{\rm A}^{2}\Bigr{)}^{p}E_{\rm in},$ (9)

Assuming $\phi_{\mathrm{at}}\ll 1$ and $\phi_{0}=0$, i.e., cavity resonance,
evaluating again analytically the geometric series, using Equation 1 and
multiplying by two to account for the fact that the intra-cavity light passes
the atoms both in the forward and backward directions, we obtain

$\displaystyle(S_{0}^{\rm int})_{\rm CE}$
$\displaystyle=\frac{2T_{1}T_{G}}{(1-\sqrt{R_{1}R_{2}}T_{G}^{2})^{2}}(S_{0}^{\mathrm{in}})_{\rm
CE}$ (10)

For a comparison at equal disturbance, we take $(S_{0}^{\rm int})_{\rm
CE}=(S_{0}^{\rm int})_{\rm SP}$, which implies the relation
$(S_{0}^{\mathrm{in}})_{\rm CE}=(1-\sqrt{R_{1}R_{2}}T_{G}^{2})^{2}\
(S_{0}^{\mathrm{in}})_{\rm SP}/(2T_{1})$ between the input powers. Using
Equation 8, the CE versus SP MSE ratio at equal disturbance is then

$\frac{\mathrm{MSE}(\hat{\phi}_{\mathrm{at}})_{\rm
CE}}{\mathrm{MSE}(\hat{\phi}_{\mathrm{at}})_{\rm SP}}=\frac{2\eta_{\rm
sys}^{\rm SP}T_{1}}{\eta_{\rm sys}^{{\cal R},{\cal
T}}(1-\sqrt{R_{1}R_{2}}T_{G}^{2})^{2}\kappa_{\mathcal{R,T}}^{2}}\equiv\zeta_{\cal
R,T}.$ (11)

Figure 2:  Rotation gain, linearity and system efficiency, for cavity-enhanced
PR schemes in transmission and reflection. Graphs show rotation gain $\kappa$,
system efficiency $\eta_{\rm sys}$ and sensitivity enhancement $\zeta_{\cal
R,T}$ for a range of atomic-medium rotation angles $\phi_{\mathrm{at}}$, as
calculated with Equation 1, Equation 5 and Equation 6 in the resonant
condition $\phi_{0}=0$. We take $T_{G}=0.975$ for the window transmission, a
typical number for anti-reflection coated glass windows. Left column shows
transmitted signal with $R_{1}=0.99$, which ensures a large forward escape
probability except at the highest values of $R_{2}$. Right column shows
reflected signal computed with $R_{2}=1$, to maximize reflected power.

###  2.6 Results and comparison between cavity schemes

Figure 2 shows the rotation gain $\kappa$, the system efficiency $\eta_{\rm
sys}$ and sensitivity enhancement $\zeta_{\cal R,T}$ for the CE PR
measurement, as a function of mirror reflectivities and for different values
of $\phi_{\mathrm{at}}$. The $\phi_{\mathrm{at}}$ values studied are between
$10^{-3}$ and $10^{-2}$ radians. As one might expect, and as predicted by
Equation 7, the cavity produces a strong enhancement when the $\sigma_{+}$ and
$\sigma_{-}$ polarizations are simultaneously resonant with the cavity. This
is possible when $\mathcal{F}\phi_{\mathrm{at}}\lesssim 1$. We note that in
applications to atomic media, $\phi_{\mathrm{at}}$ can usually be made small
by choice of detuning from resonance, with no reduction in sensitivity to the
rotation-producing effect, e.g., spin polarization or magnetic field. This is
because increased detuning reduces both $\phi_{\mathrm{at}}$ and the
disturbance due to inelastic scattering, in such a way that the constant-
disturbance sensitivity [29] is maintained [30]. We consider one-sided
cavities ($R_{1}=0.99$ for cavity in transmission or $R_{2}=1$ in reflection)
to measure the PR from one cavity output. The greater complexity of the
reflected signal is evident. Moreover, the reflected one shows a non-linear
relationship of $\psi$ to $\phi_{\mathrm{at}}$, which manifests as a
$\phi_{\mathrm{at}}$-dependent gain $\kappa_{\cal R}$, whereas $\kappa_{\cal
T}$ is nearly $\phi_{\mathrm{at}}$-independent.

We can see that the two configurations have quite different behaviors. In
reflection, $\kappa$ shows a large enhancement around a sharp, multi-lobed
resonance, which moreover changes both shape and “position,” i.e., the $R_{1}$
value at which the resonance occurs, with $\phi_{\mathrm{at}}$. In other
words, $\psi$ has a notably nonlinear relationship to $\phi_{\mathrm{at}}$ and
this is reflected also in the sensitivity enhancement $\zeta_{\cal R}$. The
resonance is due to destructive interference between the immediate reflection
at the in-coupling mirror, which does not carry any information about the
atoms, and the field exiting the cavity, which depends on
$\phi_{\mathrm{at}}$.

In the transmitted signal, the rotation gain grows monotonically with $R_{2}$
and is only weakly dependent on $\phi_{\mathrm{at}}$. This means $\psi$ is a
nearly linear function of $\phi_{\mathrm{at}}$. The system efficiency
$\eta_{\rm sys}^{{\cal T}}$ shows a broad and nearly
$\phi_{\mathrm{at}}$-independent resonance, and coincides with fairly large
values of rotation gain $\kappa_{\cal T}$ and large sensitivity gain
$\zeta_{\cal T}$.

The maximum rotation gain in reflection is larger than that obtained in
transmission. This apparent advantage is perhaps illusory: the strongest
rotation gains occur together with the smallest efficiencies. In this regime,
the rotation angle $\psi=\arcsin(S_{2}^{\mathrm{out}}/S_{0}^{\mathrm{out}})$
is large because $S_{0}$ is small, due to the destructive interference upon
reflection from a cavity near critical coupling. Such “dark fringe”
interferometric techniques can be advantageous, especially in high-power
interferometry, when the detectors’ optimal power levels are below the input
power [31]. In probing atomic systems, however, the input power is limited by
saturation of the atomic medium [29, 32, 33], leading to power levels that are
easily detected with typical photodetectors. In this scenario, there is no
advantage to reducing $S_{0}$.

##  3 Experimental implementation of cavity-enhanced polarization rotation

###  3.1 Experimental system

Figure 3: Experimental diagram. Red lines show path of
$780\text{\,}\mathrm{nm}$ probe light, blue lines show path of
$1560\text{\,}\mathrm{nm}$ locking light. Thick lines indicate propagation in
fiber. PBS - polarizing beam splitter; BS - non-polarizing beam splitter; DPD
- differential photodetector; RPD - reference photodetector; FPD - fast
photodiode; DM - dichroic mirror; CM - cavity mirror; FM - flip mirror. Dashed
components indicates elements present only in the single-pass measurement. See
text for a description of the experimental strategy.

The experimental system is shown in Figure 3. Isotopically-enriched 87Rb vapor
with no buffer gas is housed in a glass vapor cell with a single-pass length
(interior dimension of the cell) of $15\text{\,}\mathrm{cm}$ and with wedged
windows, which are anti-reflection coated for $780\text{\,}\mathrm{nm}$. The
cell is housed in a solenoid to produce an adjustable magnetic field along the
cell axis (not shown in the diagram).

Cavity mirrors CM1 (input) and CM2 (output), each have a
$0.35\text{\,}\mathrm{m}$ radius of curvature, and are spaced to create a
near-concentric cavity with an optical length of
$0.661(3)\text{\,}\mathrm{m}$, measured from free-spectral range (FSR)
measurements. CM2 is mounted on a piezo-electric transducer (PZT) for precise
control of the cavity length. The nominal TEM00 mode waists at the center of
the resonator are $141(3)\text{\,}\mathrm{\SIUnitSymbolMicro m}$ for
$780\text{\,}\mathrm{nm}$ and $200(4)\text{\,}\mathrm{\SIUnitSymbolMicro m}$
for $1560\text{\,}\mathrm{nm}$. Due to the quasi-concentric geometry, higher
transverse modes are separated in frequency with respect to the TEM00 mode
[34]. In what follows, we use superscripts $p$ and $l$ to indicate probe and
locking light, at wavelengths $\lambda_{p}=$780\text{\,}\mathrm{nm}$$ and
$\lambda_{l}=$1560\text{\,}\mathrm{nm}$$, respectively. The mirrors have
$R_{1}^{(p)}=0.990(3)$ and $R_{2}^{(p)}=0.860(1.5)$ and $R_{1}^{(l)}=0.80(2)$
and $R_{2}^{(l)}>0.999$ . The measured finesses, inferred from the width and
FSR observed are $\mathcal{F}^{(p)}=$23(2)\text{\,}$$ and
$\mathcal{F}^{(l)}=$5(1)\text{\,}$$, implying window transmission
$T_{G}^{(p)}=$0.972(6)\text{\,}$$ and $T_{G}^{(l)}=$0.76(6)\text{\,}$$.
$\mathcal{F}^{(p)}$ is measured with a photodiode placed in transmission,
instead $\mathcal{F}^{(l)}$ with a photodiode placed in reflection (not shown
in the schematic). A removable “flip mirror” inside the cavity is used to
bypass CM2 and thereby switch from cavity-based to single-pass measurement.

The cavity length is stabilized using an auxiliary Pound-Drever-Hall (PDH)
lock. $1560\text{\,}\mathrm{nm}$ light from a single-frequency erbium-doped
fiber laser (EDFL) is phase modulated at $50\text{\,}\mathrm{MHz}$ with a
waveguide electro-optic modulator (EOM), circularly polarized, mode-matched to
the cavity and injected into CM1. The reflected power is collected on a fast
photodiode and demodulated to obtain the PDH error signal. A proportional-
integral-derivative (PID) controller feeds back to the PZT to maintain cavity
resonance. The EDFL is frequency-doubled and its second harmonic is stabilized
by modulation-transfer spectroscopy [35] to the 85Rb D2 line ($F=3\to
F^{{}^{\prime}}=4$).

The probe light, from a external-cavity diode laser (ECDL), is stabilized with
respect to the EDFL second harmonic by an offset lock, tuneable over
$\approx$2\text{\,}\mathrm{GHz}$$ around the 87Rb D2 line ($F=1\to
F^{{}^{\prime}}=2$). An acousto-optic modulator (AOM) is used to chop the
probe light into pulses. The polarization is purified with a pair of
polarizing beamsplitter cubes, and half of the probe power is split off and
detected on a reference photodiode (RPD) to infer $S_{0}^{\mathrm{in}}$. The
remaining beam is spatially matched to the TEM00 cavity mode. The locking and
probe beams are shaped by independent telescopes and combined at a dichroic
mirror before the cavity (not shown in the diagram).

Figure 4: Measured rotation gain of cavity-enhanced probing. Points show the
observed cavity-enhanced rotation angle $\psi^{\mathrm{CE}}$ and single-pass
rotation angle $\psi^{\mathrm{SP}}$ for a variety of magnetic field strengths.
Green dashed line and green shaded region show best least-squares fit with a
line passing through the origin, and plus/minus one standard error of that
fit. The fit gives the experimental cavity enhancement factor $\kappa_{\cal
T}^{\mathrm{exp}}=$16(2)\text{\,}$$. Purple line and shaded region show the
model prediction $\psi^{\mathrm{CE}}=\kappa_{\cal T}\psi^{\mathrm{SP}}$
$\kappa_{\cal T}^{\mathrm{th}}$ is calculated using Equation 1, Equation 5 and
Equation 6, with experimental values for $R_{2}$ and $T_{G}$, and in the
resonant condition $\phi_{0}=0$. This finds $\kappa_{\cal
T}^{\mathrm{th}}=$14.7(13)\text{\,}$$, in agreement with the experimental
result. Inset: PR angle versus magnetic field in single-pass (red squares) and
in the cavity-enhanced scheme (blue squares). All error bars show plus/minus
one mean standard error on 120 measurements.

###  3.2 Rotation gain measurements

To measure polarization rotation, horizontally-polarized probe pulses of
$\simeq$20\text{\,}\mathrm{\SIUnitSymbolMicro s}$$ duration, tuned to
resonance with the cavity and detuned $\simeq$1.7\text{\,}\mathrm{GHz}$$ to
the red of the $F=1\to F^{\prime}=0$ transition in 87Rb, are sent through the
system. $S_{2}^{\mathrm{out}}$ is detected with the differential photodiode
(DPD) [36] and $S_{0}^{\mathrm{in}}$ is detected with the reference
photodiode. The PR angle is then estimated as
$\hat{\psi}=\arcsin({S_{2}}/{S_{0}})$.

Probe pulses of
$1.4\text{\times}{10}^{7}\text{\,}\mathrm{p}\mathrm{h}\mathrm{o}\mathrm{t}\mathrm{o}\mathrm{n}\mathrm{s}$
were used for the single-pass probing, and
$1.1\text{\times}{10}^{7}\text{\,}\mathrm{p}\mathrm{h}\mathrm{o}\mathrm{t}\mathrm{o}\mathrm{n}\mathrm{s}$
for the cavity-enhanced probing. With these numbers, the spatially-averaged
probe intensity seen by the atoms is the same in the two scenarios. We note
that, in the low-power, off-resonance probing regime employed here, the PR
angles are negligibly affected by probe-induced nonlinearities such as
saturation of the optical transition or optical pumping, so the photon number
does not affect the rotation gain. Keeping the mean intensity equal is
important, however, for a fair comparison of signal to noise ratios under
conditions of equal illumination and thus equal perturbation to the atomic
medium.

To induce polarization rotation via the Faraday effect, a magnetic field B is
applied to the atoms in a range of 0 to 10 Gauss. No optical pumping is
applied to the vapor. Figure 4 shows a comparison of the measured rotation
angles for single-pass and cavity-enhanced probing. The observed rotation gain
is linear in the applied field and in agreement with predictions from the
model of subsection 2.2.

With the setup coefficients given in subsection 3.1, the parameters $\eta_{\rm
sys}^{\rm SP}$ and $\eta_{\rm sys}^{\cal T}$ can be estimated. We used the
definition given in subsection 2.3 and subsection 2.4, where we consider
$\phi_{0}=0$ (resonant cavity condition) and $\psi^{\rm SP}\ll 1$, in
agreement to the experimental data (inset Figure 4). We obtain $\eta_{\rm
sys}^{\rm SP}=$0.945(12)\text{\,}$$ and $\eta_{\rm sys}^{\cal
T}=$0.08(3)\text{\,}$$. The parameter $\zeta_{\cal T}$ can be also calculated,
with the values obtain for $\eta_{\rm sys}^{\rm SP}$ and $\eta_{\rm sys}^{\cal
T}$ and with the values of $\kappa_{\cal T}^{\mathrm{exp}}$ or $\kappa_{\cal
T}^{\mathrm{th}}$ written in the caption of Figure 4. Using Equation 11, we
have $\zeta_{\cal T}^{\rm th}=$-11.8(10)\text{\,}\mathrm{d}\mathrm{B}$$ and
$\zeta_{\cal T}^{\rm exp}=$-12.5(14)\text{\,}\mathrm{d}\mathrm{B}$$, in good
agreement with our prediction.

Figure 5: PR noise in the single-pass and cavity-enhanced scenarios. Red
(blue) points show $\mathrm{var}(S_{2}^{\mathrm{out}})$ as a function of
$\braket{S_{0}^{\mathrm{out}}}$ in the SP (CE) scheme, as described in the
text. Curves show least square fits to the data using the model in Equation
12. Solid black lines show electronic and shot noise contributions common to
both SP and CE scenarios. Light (dark) green dot-dashed (dashed) line shows
technical noise in SP (CE) case. Horizontal and vertical error bars show
plus/minus standard error of $2\braket{S_{0}^{\mathrm{out}}}$ and variance,
respectively. The red (blue) shaded area marks the SN limited region of the
detection system in SP (CE), which is found to be
$$6\text{\times}{10}^{5}\text{\,}$\leq
2\braket{S_{0}^{\mathrm{out}}}\leq$1.3\text{\times}{10}^{8}\text{\,}$$ (
$$6\text{\times}{10}^{5}\text{\,}$\leq
2\braket{S_{0}^{\mathrm{out}}}\leq$9\text{\times}{10}^{5}\text{\,}$$ ).

###  3.3 Noise measurements

The resonant cavity, in addition to enhancing the polarization rotation
signal, can be expected to introduce additional technical noise, for example
fluctuation of the cavity transmission due to fluctuations in the cavity
length. To assess these noise effects, we measure the detection noise in the
CE and SP cases as a function of input photon number. We follow a similar
procedure to the one detailed in [36]: we send light pulses of
$100\text{\,}\mathrm{\SIUnitSymbolMicro s}$ duration with different beam
powers. There is no magnetic field applied on the atoms, so the expected FR
angle is zero. We partition each pulse in segments of different lengths from
$1\text{\,}\mathrm{\SIUnitSymbolMicro s}$ to
$100\text{\,}\mathrm{\SIUnitSymbolMicro s}$, and extract the variance of the
detected $S_{2}^{\mathrm{out}}$ and the mean of the detected power
$S_{0}^{\mathrm{out}}$ (inferred from the RPD). Then we fit with the model

$\mathrm{var}(S_{2}^{\mathrm{out}})=A+B\ \braket{S_{0}^{\mathrm{out}}}+C\
\braket{S_{0}^{\mathrm{out}}}^{2},$ (12)

where $A,$ $B,$ and $C$ describe the strength of the electronic noise, shot
noise, and technical noise contributions, respectively. Figure 5 shows the
results in both SP and CE scenarios. We can appreciate a clear increase of the
technical noise in the CE scheme, which narrows the region of shot-noise
limited detection. The CE scheme is nonetheless shot-noise limited in a range
of photon numbers suitable for probing of cold atomic ensembles [37].

##  4 Conclusions and outlook

In this work we have analyzed and demonstrated cavity-enhanced polarization
rotation measurements. We found general expressions, computed using the Jones
calculus, for the rotation enhancement for both the forward and backward
collection scenarios. Using forward collection, we demonstrate experimentally
an enhancement of $\kappa_{\cal T}^{\mathrm{exp}}=$16(2)\text{\,}$$, in good
agreement with theory predictions. Noise characterization measurements show
the system to be shot noise limited in a relevant range of pulse energies. The
technique is readily applicable to a variety of material systems, including
hot, cold and ultra-cold atomic and molecular ensembles, and will enable CE
quantum non-demolition measurements of magnetic degrees of freedom, with
potential applications in atomic magnetometers [14], gyroscopes [15], and
instruments to search for physics beyond the standard model [16, 17, 18].

Funding This work was supported by the H2020 Future and Emerging Technologies
Quantum Technologies Flagship projects MACQSIMAL (Grant Agreement No. 820393)
and QRANGE (Grant Agreement No. 820405); H2020 Marie Skłodowska-Curie Actions
project ITN ZULF-NMR (Grant Agreement No. 766402); PROBIST (Grant Agreement
No. 754510); Spanish Ministry of Science projects OCARINA (Grant No.
PGC2018-097056-B-I00 project funded by MCIN/ AEI /10.13039/501100011033/ FEDER
“A way to make Europe”), Q-CLOCKS (Grant No. PCI2018-092973 project funded by
MCIN/ AEI /10.13039/501100011033/ FEDER “A way to make Europe”), and “Severo
Ochoa” Center of Excellence CEX2019-000910-S Generalitat de Catalunya through
the CERCA program; Agència de Gestió d’Ajuts Universitaris i de Recerca Grant
No. 2017-SGR-1354; Secretaria d’Universitats i Recerca del Departament
d’Empresa i Coneixement de la Generalitat de Catalunya, co-funded by the
European Union Regional Development Fund within the ERDF Operational Program
of Catalunya (project QuantumCat, ref. 001-P-001644); Fundació Privada Cellex;
Fundació Mir-Puig; 17FUN03 USOQS, which has received funding from the EMPIR
programme co-financed by the Participating States and from the European
Union’s Horizon 2020 research and innovation programme. This project has
received funding from the European Union’s Horizon 2020 research and
innovation programme under the Marie Skłodowska-Curie grant agreement No
665884.

Disclosures The authors declare no conflicts of interest.

Data availability Data underlying the results presented in this paper are not
publicly available at this time but may be obtained from the authors upon
reasonable request.

## References

  * [1] I. Kominis, T. Kornack, J. Allred, and M. Romalis, “A subfemtotesla multichannel atomic magnetometer,” Nature 422, 596–599 (2003).
  * [2] W. C. Griffith, S. Knappe, and J. Kitching, “Femtotesla atomic magnetometry in a microfabricated vapor cell,” Optics Express 18, 27167–27172 (2010).
  * [3] B. Julsgaard, A. Kozhekin, and E. S. Polzik, “Experimental long-lived entanglement of two macroscopic objects,” Nature 413, 400–403 (2001).
  * [4] T. Takano, M. Fuyama, R. Namiki, and Y. Takahashi, “Spin squeezing of a cold atomic ensemble with the nuclear spin of one-half,” Physical Review Letters 102, 033601 (2009).
  * [5] R. J. Sewell, M. Koschorreck, M. Napolitano, B. Dubost, N. Behbood, and M. W. Mitchell, “Magnetic sensitivity beyond the projection noise limit by spin squeezing,” Physical Review Letters 109, 253605 (2012).
  * [6] M. H. Schleier-Smith, I. D. Leroux, and V. Vuletic, “States of an ensemble of two-level atoms with reduced quantum uncertainty,” Physical Review Letters 104, 073604 (2010).
  * [7] J. G. Bohnet, K. C. Cox, M. A. Norcia, J. M. Weiner, Z. Chen, and J. K. Thompson, “Reduced spin measurement back-action for a phase sensitivity ten times beyond the standard quantum limit,” Nat Photon 8, 731–736 (2014).
  * [8] O. Hosten, N. J. Engelsen, R. Krishnakumar, and M. A. Kasevich, “Measurement noise 100 times lower than the quantum-projection limit using entangled atoms,” Nature 529, 505–508 (2016).
  * [9] I. D. Leroux, M. H. Schleier-Smith, and V. Vuletić, “Orientation-dependent entanglement lifetime in a squeezed atomic clock,” Phys. Rev. Lett. 104, 250801 (2010).
  * [10] G. Vallet, E. Bookjans, U. Eismann, S. Bilicki, R. L. Targat, and J. Lodewyck, “A noise-immune cavity-assisted non-destructive detection for an optical lattice clock in the quantum regime,” New Journal of Physics 19, 083002 (2017).
  * [11] H. Crepaz, L. Y. Ley, and R. Dumke, “Cavity enhanced atomic magnetometry,” Scientific Reports 5, 15448 (2015).
  * [12] C.-Y. Chang and J.-T. Shy, “Cavity-enhanced faraday rotation measurement with auto-balanced photodetection,” Applied optics 54, 8526–8530 (2015).
  * [13] V. Lucivero, W. Lee, N. Dural, and M. Romalis, “Femtotesla direct magnetic gradiometer using a single multipass cell,” Physical Review Applied 15, 014004 (2021).
  * [14] D. Budker and M. Romalis, “Optical magnetometry,” Nature Phys. 3, 227–234 (2007).
  * [15] T. W. Kornack, R. K. Ghosh, and M. V. Romalis, “Nuclear spin gyroscope based on an atomic comagnetometer,” Physical Review Letters 95, 230801 (2005).
  * [16] J. Lee, A. Almasi, and M. Romalis, “Improved limits on spin-mass interactions,” Physical Review Letters 120, 161801 (2018).
  * [17] P. Gomez, F. Martin, C. Mazzinghi, D. B. Orenes, S. Palacios, and M. W. Mitchell, “Bose-einstein condensate comagnetometer,” Physical review letters 124, 170401 (2020).
  * [18] M. Safronova, D. Budker, D. DeMille, D. F. J. Kimball, A. Derevianko, and C. W. Clark, “Search for new physics with atoms and molecules,” Reviews of Modern Physics 90, 025008 (2018).
  * [19] D. Budker, W. Gawlik, D. F. Kimball, S. M. Rochester, V. V. Yashchuk, and A. Weis, “Resonant nonlinear magneto-optical effects in atoms,” Rev. Mod. Phys. 74, 1153–1201 (2002).
  * [20] Y. Takahashi, K. Honda, N. Tanaka, K. Toyoda, K. Ishikawa, and T. Yabuzaki, “Quantum nondemolition measurement of spin via the paramagnetic faraday rotation,” Phys. Rev. A 60, 4974–4979 (1999).
  * [21] R. C. Jones, “A new calculus for the treatment of optical systems. i. description and discussion of the calculus,” Josa 31, 488–493 (1941).
  * [22] E. Hecht and A. Zajac, _Optics_ , vol. 4 (Addison Wesley San Francisco, 2002).
  * [23] R. J. Sewell, M. Koschorreck, M. Napolitano, B. Dubost, N. Behbood, and M. W. Mitchell, “Magnetic Sensitivity Beyond the Projection Noise Limit by Spin Squeezing,” Physical Review Letters 109, 253605 (2012).
  * [24] G. Colangelo, C. F. Martin, B. L. C., S. R. J., and M. M. W., “Simultaneous tracking of spin angle and amplitude beyond classical limits,” Nature 543, 525–528 (2017).
  * [25] R. J. Sewell, M. Napolitano, N. Behbood, G. Colangelo, and M. W. Mitchell, “Certified quantum non-demolition measurement of a macroscopic material system,” Nat Photon 7, 517–520 (2013).
  * [26] J. Kong, R. Jiménez-Martínez, C. Troullinou, V. G. Lucivero, G. Tóth, and M. W. Mitchell, “Measurement-induced, spatially-extended entanglement in a hot, strongly-interacting atomic system,” Nature Communications 11, 2415 (2020).
  * [27] F. Martin Ciurana, G. Colangelo, L. Slodička, R. J. Sewell, and M. W. Mitchell, “Entanglement-enhanced radio-frequency field detection and waveform sensing,” Physical Review Letters 119, 043603 (2017).
  * [28] G. G. Stokes, “On the composition and resolution of streams of polarized light from different sources,” Transactions of the Cambridge Philosophical Society 9, 399 (1851).
  * [29] F. Wolfgramm, C. Vitelli, F. A. Beduini, N. Godbout, and M. W. Mitchell, “Entanglement-enhanced probing of a delicate material system,” Nat Photon 7, 28–32 (2013). 10.1038/nphoton.2012.300.
  * [30] M. Koschorreck, M. Napolitano, B. Dubost, and M. W. Mitchell, “Sub-Projection-Noise Sensitivity in Broadband Atomic Magnetometry,” Physical Review Letters 104, 093602 (2010).
  * [31] T. L. S. Collaboration, “Advanced ligo,” Classical and Quantum Gravity 32, 074001 (2015).
  * [32] M. W. Mitchell, “Number-unconstrained quantum sensing,” Quantum Science and Technology 2, 044005 (2017).
  * [33] V. G. Lucivero, A. Dimic, J. Kong, R. Jiménez-Martínez, and M. W. Mitchell, “Sensitivity, quantum limits, and quantum enhancement of noise spectroscopies,” Phys. Rev. A 95, 041803(R) (2017).
  * [34] B. E. Saleh and M. C. Teich, _Fundamentals of photonics_ (john Wiley & sons, 2019).
  * [35] Y. N. M. de Escobar, S. P. Álvarez, S. Coop, T. Vanderbruggen, K. T. Kaczmarek, and M. W. Mitchell, “Absolute frequency references at 1529 and 1560 nm using modulation transfer spectroscopy,” Optics letters 40, 4731–4734 (2015).
  * [36] F. M. Ciurana, G. Colangelo, R. J. Sewell, and M. W. Mitchell, “Real-time shot-noise-limited differential photodetection for atomic quantum control,” Opt. Lett. 41, 2946–2949 (2016).
  * [37] M. Kubasik, M. Koschorreck, M. Napolitano, S. R. de Echaniz, H. Crepaz, J. Eschner, E. S. Polzik, and M. W. Mitchell, “Polarization-based light-atom quantum interface with an all-optical trap,” Phys. Rev. A 79, 043815 (2009).