0$
\end_inset
.
\end_layout
\begin_layout Itemize
Each event occurring between consecutive
\begin_inset Formula $\bar{Z}$
\end_inset
and
\begin_inset Formula $Z$
\end_inset
events occurs at
\begin_inset Formula $x<0$
\end_inset
.
\end_layout
\begin_layout Standard
If the solved system of equations satisfies the set of inequalities then
the sequence is legal, for a given
\begin_inset Formula $\left(b,\tau\right)$
\end_inset
.
These techniques allow us to use symbolic sequences to study solutions
systematically.
Solving for the times and positions associated with events in a sequence,
we can plot the solution represented by the sequence.
Additionally, by changing the inequalities to equalities, we obtain the
bifurcation curves shown in Fig.
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:bifurcations"
\end_inset
.
\end_layout
\begin_layout Subsection
Calculation of the bifurcation curve of
\begin_inset Formula $\mathbb{P}$
\end_inset
\begin_inset CommandInset label
LatexCommand label
name "subsec:Calculation-of-the"
\end_inset
\end_layout
\begin_layout Standard
For
\begin_inset Formula $\tau>\frac{1}{2}$
\end_inset
,
\begin_inset Formula $\mathbb{P}$
\end_inset
can written as
\begin_inset Formula
\begin{equation}
\mathbb{P}\left(S_{z}\right)=AS_{z}+B
\end{equation}
\end_inset
Due to the sparsity of
\begin_inset Formula $A$
\end_inset
, the characteristic equation can be readily calculated as
\begin_inset Formula
\begin{equation}
(-1-\lambda)(-\lambda)^{n-1}+(-1)^{n-1}\frac{2(-1)^{n-1}}{b+(-1)^{n-1}}=0
\end{equation}
\end_inset
This simplifies to
\begin_inset Formula
\begin{equation}
\lambda^{n}+\lambda^{n-1}-\frac{2(-1)^{n-1}}{b+(-1)^{n-1}}=0
\end{equation}
\end_inset
By substituting
\begin_inset Formula $\lambda=e^{i\rho}$
\end_inset
, we can solve for
\begin_inset Formula $b$
\end_inset
to determine
\begin_inset Formula $b=b_{\text{bif }}$
\end_inset
for which the fixed point of
\begin_inset Formula $\mathbb{P}$
\end_inset
is bifurcating as
\begin_inset Formula
\begin{equation}
e^{i\rho n}+e^{i\rho\left(n-1\right)}-\frac{2(-1)^{n-1}}{b_{\text{bif}}+(-1)^{n-1}}=0\label{eq:a}
\end{equation}
\end_inset
Note that as
\begin_inset Formula $\frac{2(-1)^{n-1}}{b_{\text{bif}}+(-1)^{n-1}}\in\mathbb{R}$
\end_inset
,
\begin_inset Formula
\begin{equation}
e^{i\rho n}=\overline{e^{i\rho\left(n-1\right)}}\label{eq:b}
\end{equation}
\end_inset
So,
\begin_inset Formula $e^{i\rho\left(2n-1\right)}=1=e^{i2\pi k},k\in\mathbb{Z}$
\end_inset
, and so
\begin_inset Formula $\rho=\frac{2\pi k}{2n-1}$
\end_inset
.
In order to satisfy Eq.
(
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:b"
\end_inset
),
\begin_inset Formula $k=n-1$
\end_inset
.
Substituting
\begin_inset Formula $\rho$
\end_inset
back into Eq.
(
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:a"
\end_inset
),
\begin_inset Formula
\begin{equation}
\cos\left(\frac{2\pi n(n-1)}{2n-1}\right)+\cos\left(\frac{2\pi\left(n-1\right)^{2}}{2n-1}\right)-\frac{2(-1)^{n-1}}{b_{\text{bif}}+(-1)^{n-1}}=0
\end{equation}
\end_inset
which simplifies under a sum-to-product cosine identity to
\begin_inset Formula
\begin{equation}
\cos\left(\frac{\pi\left(n-1\right)}{2n-1}\right)-\frac{1}{b_{\text{bif}}+(-1)^{n-1}}=0
\end{equation}
\end_inset
Therefore,
\begin_inset Formula
\begin{equation}
b_{\text{bif}}(n)=\frac{1}{\cos\left(\frac{\pi\left(n-1\right)}{2n-1}\right)}-(-1)^{n-1}
\end{equation}
\end_inset
where
\begin_inset Formula $n=\left\lceil 2\tau\right\rceil $
\end_inset
.
\end_layout
\begin_layout Section
Dynamics of a band-pass filter system with switched time-delayed feedback
\end_layout
\begin_layout Subsection
Derivation of the bandpass-filter system
\begin_inset CommandInset label
LatexCommand label
name "subsec:Derivation-of-the"
\end_inset
\end_layout
\begin_layout Standard
A low-pass filter is a system or device that passes frequencies lower than
the corner frequency
\begin_inset Formula $\omega_{L}$
\end_inset
and attenuates signals with frequency higher than
\begin_inset Formula $\omega_{L}$
\end_inset
\begin_inset CommandInset citation
LatexCommand cite
key "udaltsov2002bandpass,blakely2004high,illing2005hopf"
literal "false"
\end_inset
.
A high-pass filter is a system or device that passes frequencies higher
than the corner frequency
\begin_inset Formula $\omega_{H}$
\end_inset
and attenuates signals with frequency lower than
\begin_inset Formula $\omega_{H}$
\end_inset
.
The dynamics of low-pass and high-pass filter systems are described in
Eq.
(
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:lowpass"
plural "false"
caps "false"
noprefix "false"
\end_inset
) and Eq.(
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:lowpass"
plural "false"
caps "false"
noprefix "false"
\end_inset
) respectively
\begin_inset CommandInset citation
LatexCommand cite
key "udaltsov2002bandpass,illing2005hopf,blakely2004high"
literal "false"
\end_inset
\end_layout
\begin_layout Standard
\begin_inset Formula
\begin{align}
\tau_{L}\dot{x}_{L}+x_{L} & =f\left(t\right)\label{eq:lowpass}\\
\dot{x}_{H}+\frac{x_{H}}{\tau_{H}} & =\frac{d}{dt}\left[g\left(t\right)\right]\label{eq:highpass}
\end{align}
\end_inset
where
\begin_inset Formula $\tau_{L}=\omega_{L}^{-1}$
\end_inset
and
\begin_inset Formula $\tau_{H}=\omega_{H}^{-1}$
\end_inset
, and
\begin_inset Formula $f$
\end_inset
and
\begin_inset Formula $g$
\end_inset
are the input signals which drive the filters.
When undriven, the filters relax to
\begin_inset Formula $x_{L}=x_{H}=0$
\end_inset
.
When we drive Eq.
(
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:lowpass"
plural "false"
caps "false"
noprefix "false"
\end_inset
) and Eq.
(
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:highpass"
plural "false"
caps "false"
noprefix "false"
\end_inset
) with input signal
\begin_inset Formula $f\left(t\right)=g\left(t\right)=\cos\left(\omega t\right)$
\end_inset
, the outputs
\begin_inset Formula $x_{L}\left(t\right)$
\end_inset
and
\begin_inset Formula $x_{H}\left(t\right)$
\end_inset
are signals signal with frequency
\begin_inset Formula $\omega$
\end_inset
whose amplitude depends on
\begin_inset Formula $\omega$
\end_inset
relative to the corner frequency, as shown in Fig.
\begin_inset CommandInset ref
LatexCommand ref
reference "fig:High-pass-and-low-pass"
plural "false"
caps "false"
noprefix "false"
\end_inset
\begin_inset Formula $\left(a\right)$
\end_inset
and
\begin_inset Formula $\left(b\right)$
\end_inset
respectively.
\begin_inset Float figure
wide false
sideways false
status open
\begin_layout Plain Layout
\align center
\begin_inset Graphics
filename ../lucas_system/fig_lowhighpass_filter.pdf
width 75text%
\end_inset
\end_layout
\begin_layout Plain Layout
\begin_inset Caption Standard
\begin_layout Plain Layout
Low-pass and high-pass filter output amplitude with input
\begin_inset Formula $f$
\end_inset
\begin_inset Formula $\left(t\right)=\cos\left(\omega t\right)$
\end_inset
for varying input frequency
\begin_inset Formula $\omega$
\end_inset
, where the corner frequency
\begin_inset Formula $\omega_{L}=\omega_{H}=1$
\end_inset
.
\begin_inset CommandInset label
LatexCommand label
name "fig:High-pass-and-low-pass"
\end_inset
\end_layout
\end_inset
\end_layout
\end_inset
A band-pass filter can be constructed by coupling a low-pass filter with
a high-pass filter
\begin_inset CommandInset citation
LatexCommand cite
key "udaltsov2002bandpass,blakely2004high,illing2005hopf"
literal "false"
\end_inset
.
This is achieved by letting
\begin_inset Formula $g\left(t\right)=x_{L}\left(t\right)$
\end_inset
, yielding the equations
\begin_inset Formula
\begin{align}
\tau_{L}\dot{x}_{L}+x_{L} & =f\left(t\right)\label{eq:lowpass-2}\\
\dot{x}_{H}+\frac{x_{H}}{\tau_{H}} & =\dot{x}_{L}\label{eq:highpass-2}
\end{align}
\end_inset
The system can be simplified by first integrating Eq.
(
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:highpass-2"
plural "false"
caps "false"
noprefix "false"
\end_inset
) to yield
\begin_inset Formula
\begin{equation}
x_{H}+\frac{1}{\tau_{H}}\intop_{0}^{t}x_{H}\left(t'\right)dt'=x_{L}\label{eq:integrated}
\end{equation}
\end_inset
with the integration constant being satisfied by the lower bound of the
integral.
Substituting Eq.
(
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:highpass-2"
plural "false"
caps "false"
noprefix "false"
\end_inset
) and Eq.
(
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:integrated"
plural "false"
caps "false"
noprefix "false"
\end_inset
) into Eq.
(
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:lowpass-2"
plural "false"
caps "false"
noprefix "false"
\end_inset
) yields
\begin_inset Formula
\[
\tau_{L}\left[\dot{x}_{H}+\frac{x_{H}}{\tau_{H}}\right]+\left[x_{H}+\frac{1}{\tau_{H}}\intop_{0}^{t}x_{H}\left(s\right)ds\right]=f\left(t\right)
\]
\end_inset
which we rewrite as
\begin_inset Formula
\begin{align*}
x_{H}+\frac{\tau_{H}\tau_{L}}{\left(\tau_{H}+\tau_{L}\right)}\dot{x}_{H}+\frac{1}{\left(\tau_{H}+\tau_{L}\right)}\intop_{0}^{t}x_{H}\left(s\right)ds & =\frac{\tau_{H}}{\left(\tau_{H}+\tau_{L}\right)}f\left(t\right)
\end{align*}
\end_inset
where
\begin_inset Formula $\frac{\tau_{H}}{\left(\tau_{H}+\tau_{L}\right)}$
\end_inset
is the signal gain due to passing through the band-pass filter.
We substitute centre frequency
\begin_inset Formula $\Omega=\frac{1}{\sqrt{\tau_{H}\tau_{L}}}$
\end_inset
and quality factor
\begin_inset Formula $Q=\frac{\sqrt{\tau_{H}\tau_{L}}}{\left(\tau_{H}+\tau_{L}\right)}$
\end_inset
to yield
\begin_inset Formula
\begin{align}
x_{H}+\frac{Q}{\Omega}\frac{dx_{H}}{d\bar{t}}+Q\Omega\intop_{0}^{t}x_{H}\left(s\right)ds & =\frac{\tau_{H}}{\left(\tau_{H}+\tau_{L}\right)}f\left(t\right)\label{eq:bandpass}
\end{align}
\end_inset
We rescale
\begin_inset Formula $x_{H}$
\end_inset
to get dimensionless output
\begin_inset Formula $x=x_{H}\left(\frac{\tau_{H}}{\left(\tau_{H}+\tau_{L}\right)}\right)^{-1}$
\end_inset
, yielding the equation
\begin_inset Formula
\begin{equation}
x+\frac{Q}{\Omega}\frac{dx}{dt}+Q\Omega\intop^{t}x\left(s\right)ds=f\left(t\right)
\end{equation}
\end_inset
Finally, we introduce a variable
\begin_inset Formula $y=Q\Omega\intop^{t}x\left(s\right)ds$
\end_inset
such that
\begin_inset Formula $\frac{dy}{dt}=Q\Omega x$
\end_inset
, yielding the non-smooth delay differential equation
\begin_inset Formula
\begin{equation}
\begin{aligned}Q\Omega^{-1}\dot{x} & =-x-y+f\left(t\right)\\
\dot{y} & =Q\Omega x
\end{aligned}
\label{eq:dde-1}
\end{equation}
\end_inset
previously studied in
\begin_inset CommandInset citation
LatexCommand cite
key "udaltsov2002bandpass,illing2005hopf,blakely2004high"
literal "false"
\end_inset
with different input signals
\begin_inset Formula $f$
\end_inset
\begin_inset Formula $\left(t\right)$
\end_inset
.
\end_layout
\begin_layout Subsection
Useful properties of Pauli matrices
\end_layout
\begin_layout Standard
The Pauli matrices, with the identify matrix
\begin_inset Formula $I$
\end_inset
, form a basis for the vector space of
\begin_inset Formula $2\times2$
\end_inset
real matrices
\begin_inset CommandInset citation
LatexCommand cite
key "condon1929quantum"
literal "false"
\end_inset
.
They are
\begin_inset Formula
\begin{align}
I & =\left(\begin{array}{cc}
1 & 0\\
0 & 1
\end{array}\right)\\
\sigma_{x} & =\left(\begin{array}{cc}
0 & 1\\
1 & 0
\end{array}\right)\\
\sigma_{y} & =\left(\begin{array}{cc}
0 & -i\\
i & 0
\end{array}\right)\\
\sigma_{z} & =\left(\begin{array}{cc}
1 & 0\\
0 & -1
\end{array}\right)
\end{align}
\end_inset
Pauli matrices have some very useful properties.
Firstly, the square of any Pauli matrix
\begin_inset Formula $\sigma_{j}$
\end_inset
yields the identity matrix
\begin_inset Formula
\begin{equation}
\left(\sigma_{i}\right)^{2}=I
\end{equation}
\end_inset
Secondly, Pauli matrices anti-commute; for a pair of different Pauli matrices
\begin_inset Formula $\sigma_{j}$
\end_inset
and
\begin_inset Formula $\sigma_{k}$
\end_inset
,
\begin_inset Formula $j\neq k$
\end_inset
,
\begin_inset Formula
\begin{equation}
\sigma_{j}\sigma_{k}=-\sigma_{k}\sigma_{j}
\end{equation}
\end_inset
These two properties can be summarised in Einstein notation
\begin_inset CommandInset citation
LatexCommand cite
key "einstein1916die"
literal "false"
\end_inset
as
\begin_inset Formula
\begin{equation}
\sigma_{j}\sigma_{k}=\delta_{jk}I+i\epsilon_{jkl}\sigma_{l}
\end{equation}
\end_inset
As the Pauli matrices and the identify matrix form a basis for the vector
space of
\begin_inset Formula $2\times2$
\end_inset
real matrices, a
\begin_inset Formula $2\times2$
\end_inset
real matrix
\begin_inset Formula $A$
\end_inset
can be written as
\begin_inset Formula
\begin{align}
A & =a_{0}I+\bar{a}\cdot\bar{\sigma}\\
& =a_{0}I+a_{x}\sigma_{x}+a_{y}\sigma_{y}+a_{z}\sigma_{z}
\end{align}
\end_inset
where
\begin_inset Formula
\begin{align}
\bar{a} & =\left(a_{x},a_{y},a_{z}\right)\\
\bar{\sigma} & =\left(\sigma_{x},\sigma_{y},\sigma_{z}\right)
\end{align}
\end_inset
A further useful property can be obtained from the Einstein notation
\begin_inset Formula
\begin{align}
\left(\bar{a}\cdot\bar{\sigma}\right)^{2} & =a_{j}\sigma_{j}a_{k}\sigma_{k}\\
& =a_{j}a_{k}\sigma_{j}\sigma_{k}\\
& =a_{j}a_{k}\left(\delta_{jk}I+i\epsilon_{jkl}\sigma_{l}\right)\\
& =a_{j}a_{k}\delta_{jk}I+ia_{j}a_{k}\epsilon_{jkl}\sigma_{l}\\
& =a_{j}a_{j}I+ia_{j}a_{k}\epsilon_{jkl}\sigma_{l}
\end{align}
\end_inset
The second term term collapses to
\begin_inset Formula $0$
\end_inset
as
\begin_inset Formula
\begin{align}
a_{j}a_{k}\epsilon_{jkl}\sigma_{l} & =\epsilon_{jkl}a_{j}a_{k}\sigma_{l}\\
& =\epsilon_{kjl}a_{k}a_{j}\sigma_{l}\text{ by swapping order of summation}\\
& =-\epsilon_{jkl}a_{k}a_{j}\sigma_{l}\text{ by permutation of \ensuremath{\epsilon_{jkl}}}\\
& =-\epsilon_{jkl}a_{j}a_{k}\sigma_{l}\text{ by commutativity of scalar multiplication}\\
& \rightarrow\epsilon_{jkl}a_{j}a_{k}\sigma_{l}=-\epsilon_{jkl}a_{j}a_{k}\sigma_{l}=0
\end{align}
\end_inset
Thus,
\begin_inset Formula
\begin{align}
\left(\bar{a}\cdot\bar{\sigma}\right)^{2} & =a_{j}a_{j}I\label{eq:magnitude pauli vector}
\end{align}
\end_inset
\end_layout
\begin_layout Subsection
Calculating the matrix exponential
\begin_inset Formula $e^{At}$
\end_inset
\begin_inset CommandInset label
LatexCommand label
name "subsec:Calculating-the-matrix"
\end_inset
\end_layout
\begin_layout Standard
We can decompose
\begin_inset Formula $A$
\end_inset
from Eq.
(
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:matrixA"
plural "false"
caps "false"
noprefix "false"
\end_inset
) in terms of Pauli Matrices as
\begin_inset Formula
\begin{equation}
A=-\frac{\Omega}{2Q}I+\frac{\Omega}{2Q}\left(Q^{2}-1\right)\sigma_{x}-\frac{i\Omega}{2Q}\left(Q^{2}+1\right)\sigma_{y}-\frac{\Omega}{2Q}\sigma_{z}
\end{equation}
\end_inset
or
\end_layout
\begin_layout Standard
\begin_inset Formula
\begin{equation}
A=a_{0}I+\bar{a}\cdot\bar{\sigma}
\end{equation}
\end_inset
where
\begin_inset Formula
\begin{equation}
\begin{aligned}\bar{a} & =\left(\frac{\Omega}{2Q}\left(Q^{2}-1\right),-\frac{i\Omega}{2Q}\left(Q^{2}+1\right),-\frac{\Omega}{2Q}\right)\\
a_{0} & =-\frac{\Omega}{2Q}
\end{aligned}
\label{eq:A_vectors}
\end{equation}
\end_inset
For our convenience later on, we rewrite
\begin_inset Formula $\bar{a}$
\end_inset
to get
\begin_inset Formula
\begin{equation}
A=a_{0}I+\bar{a}\cdot\bar{\sigma}=a_{0}I+c_{a}\hat{a}\cdot\bar{\sigma}\label{eq:rewrite abar}
\end{equation}
\end_inset
where we require
\begin_inset Formula $\hat{a}$
\end_inset
to be a unit vector such that
\begin_inset Formula
\begin{equation}
\left(\hat{a}\cdot\bar{\sigma}\right)^{2}=I
\end{equation}
\end_inset
Using Eq.
(
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:magnitude pauli vector"
plural "false"
caps "false"
noprefix "false"
\end_inset
) and Eq.
(
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:rewrite abar"
plural "false"
caps "false"
noprefix "false"
\end_inset
), we can calculate
\begin_inset Formula $c_{a}$
\end_inset
as
\begin_inset Formula
\begin{align}
\left(\bar{a}\cdot\bar{\sigma}\right)^{2} & =c_{a}^{2}\left(\hat{a}\cdot\bar{\sigma}\right)^{2}\\
a_{j}a_{j}I & =c_{a}^{2}I\\
c_{a} & =\sqrt{a_{j}a_{j}}\\
& =\sqrt{\left[\frac{\Omega}{2Q}\left(Q^{2}-1\right)\right]^{2}+\left[-\frac{i\Omega}{2Q}\left(Q^{2}+1\right)\right]^{2}+\left(-\frac{\Omega}{2Q}\right)^{2}}\\
& =\frac{\Omega}{2Q}\sqrt{1-4Q^{2}}\label{eq:ca}
\end{align}
\end_inset
Then
\begin_inset Formula
\begin{align}
\hat{a} & =\frac{1}{c_{a}}\bar{a}\\
& =\frac{1}{\frac{\Omega}{2Q}\sqrt{1-4Q^{2}}}\left(\frac{\Omega}{2Q}\left(Q^{2}-1\right),-\frac{i\Omega}{2Q}\left(Q^{2}+1\right),-\frac{\Omega}{2Q}\right)\\
& =\frac{1}{\sqrt{1-4Q^{2}}}\left(Q^{2}-1,-i\left(Q^{2}+1\right),-1\right)
\end{align}
\end_inset
so that
\begin_inset Formula
\begin{align}
\hat{a}\cdot\bar{\sigma} & =\frac{1}{\sqrt{1-4Q^{2}}}\left(Q^{2}-1,-i\left(Q^{2}+1\right),-1\right)\cdot\left(\sigma_{x},\sigma_{y},\sigma_{z}\right)\\
& =\frac{1}{\sqrt{1-4Q^{2}}}\left[\left(Q^{2}-1\right)\left(\begin{array}{cc}
0 & 1\\
1 & 0
\end{array}\right)-i\left(Q^{2}+1\right)\left(\begin{array}{cc}
0 & -i\\
i & 0
\end{array}\right)-1\left(\begin{array}{cc}
1 & 0\\
0 & -1
\end{array}\right)\right]\\
& =\frac{1}{\sqrt{1-4Q^{2}}}\left(\begin{array}{cc}
-1 & -2\\
2Q^{2} & 1
\end{array}\right)\label{eq:normalised pauli vector}
\end{align}
\end_inset
Using Eq.
(
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:rewrite abar"
plural "false"
caps "false"
noprefix "false"
\end_inset
) and Eq.
(
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:normalised pauli vector"
plural "false"
caps "false"
noprefix "false"
\end_inset
), we can now calculate the matrix exponential
\end_layout
\begin_layout Standard
\begin_inset Formula
\begin{align}
e^{At} & =\exp\left(ta_{0}I+t\bar{a}\cdot\bar{\sigma}\right)\\
& =\exp\left(ta_{0}\right)\exp\left(t\bar{a}\cdot\bar{\sigma}\right)\\
& =\exp\left(ta_{0}\right)\exp\left(tc_{a}\hat{a}\cdot\bar{\sigma}\right)\\
& =\exp\left(ta_{0}\right)\sum_{k=0}^{\infty}\frac{1}{k!}\left(tc_{a}\hat{a}\cdot\bar{\sigma}\right)^{k}\\
& =\exp\left(ta_{0}\right)\left\{ \sum_{n=0}^{\infty}\frac{1}{\left(2n\right)!}\left(tc_{a}\hat{a}\cdot\bar{\sigma}\right)^{2n}+\sum_{n=0}^{\infty}\frac{1}{\left(2n+1\right)!}\left(tc_{a}\hat{a}\cdot\bar{\sigma}\right)^{2n+1}\right\} \\
& =\exp\left(ta_{0}\right)\left\{ \sum_{n=0}^{\infty}\frac{1}{\left(2n\right)!}\left(tc_{a}\right)^{2n}\left(\hat{a}\cdot\bar{\sigma}\right)^{2n}+\sum_{n=0}^{\infty}\frac{1}{\left(2n+1\right)!}\left(tc_{a}\right)^{2n+1}\left(\hat{a}\cdot\bar{\sigma}\right)^{2n}\left(\hat{a}\cdot\bar{\sigma}\right)\right\} \\
& =\exp\left(ta_{0}\right)\left\{ \sum_{n=0}^{\infty}\frac{1}{\left(2n\right)!}\left(tc_{a}\right)^{2n}I^{n}+\sum_{n=0}^{\infty}\frac{1}{\left(2n+1\right)!}\left(tc_{a}\right)^{2n+1}I^{n}\left(\hat{a}\cdot\bar{\sigma}\right)\right\} \\
& =\exp\left(ta_{0}\right)\left\{ I\sum_{n=0}^{\infty}\frac{1}{\left(2n\right)!}\left(tc_{a}\right)^{2n}+\left(\hat{a}\cdot\bar{\sigma}\right)\sum_{n=0}^{\infty}\frac{1}{\left(2n+1\right)!}\left(tc_{a}\right)^{2n+1}\left(\hat{a}\cdot\bar{\sigma}\right)\right\} \\
& =\exp\left(ta_{0}\right)\left\{ I\cosh\left(tc_{a}\right)+\left(\hat{a}\cdot\bar{\sigma}\right)\sinh\left(tc_{a}\right)\right\}
\end{align}
\end_inset
We can substitute from
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:normalised pauli vector"
plural "false"
caps "false"
noprefix "false"
\end_inset
and
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:ca"
plural "false"
caps "false"
noprefix "false"
\end_inset
to get the matrix exponential in the overdamped regime
\begin_inset Formula
\begin{equation}
\begin{aligned}e^{At} & =e^{-\frac{\Omega}{2Q}t}\cosh\left(t\frac{\Omega}{2Q}\sqrt{1-4Q^{2}}\right)\left(\begin{array}{cc}
1 & 0\\
0 & 1
\end{array}\right)\\
& +e^{-\frac{\Omega}{2Q}t}\frac{1}{\sqrt{1-4Q^{2}}}\sinh\left(t\frac{\Omega}{2Q}\sqrt{1-4Q^{2}}\right)\left(\begin{array}{cc}
-1 & -2\\
2Q^{2} & 1
\end{array}\right)
\end{aligned}
\label{eq:exponential_overdamped}
\end{equation}
\end_inset
In the underdamped regime, the
\begin_inset Formula $\sinh$
\end_inset
and
\begin_inset Formula $\cosh$
\end_inset
terms switch over to
\begin_inset Formula $\sin$
\end_inset
and
\begin_inset Formula $\cos$
\end_inset
terms, as
\begin_inset Formula
\begin{align}
\sinh x & =-i\sin\left(ix\right)\\
\cosh x & =\cos\left(ix\right)
\end{align}
\end_inset
yielding
\begin_inset Formula
\begin{equation}
\begin{aligned}e^{At} & =e^{-\frac{\Omega}{2Q}t}\cos\left(t\frac{\Omega}{2Q}\sqrt{4Q^{2}-1}\right)\left(\begin{array}{cc}
1 & 0\\
0 & 1
\end{array}\right)\\
& +e^{-\frac{\Omega}{2Q}t}\frac{1}{\sqrt{4Q^{2}-1}}\sin\left(t\frac{\Omega}{2Q}\sqrt{4Q^{2}-1}\right)\left(\begin{array}{cc}
-1 & -2\\
2Q^{2} & 1
\end{array}\right)
\end{aligned}
\label{eq:exponential_underdamped}
\end{equation}
\end_inset
\end_layout
\begin_layout Paragraph
A note on convergence
\end_layout
\begin_layout Standard
At first glance, the second term in Eq.
(
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:exponential_overdamped"
plural "false"
caps "false"
noprefix "false"
\end_inset
) and Eq.
(
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:exponential_underdamped"
plural "false"
caps "false"
noprefix "false"
\end_inset
) respectively may have a convergence issue as
\begin_inset Formula $Q\rightarrow\frac{1}{2}$
\end_inset
and
\begin_inset Formula $\pm\left(1-4Q^{2}\right)\rightarrow0$
\end_inset
.
We can doublecheck this by calculating the limit of
\begin_inset Formula $\frac{\sinh\left(rx\right)}{x}$
\end_inset
as
\begin_inset Formula $x\rightarrow\infty$
\end_inset
where
\begin_inset Formula $r=t\frac{\Omega}{2Q}$
\end_inset
.
\begin_inset Formula
\begin{align}
\underset{x\rightarrow0}{\lim}\frac{\sinh\left(rx\right)}{x} & =\underset{x\rightarrow0}{\lim}\frac{1}{x}\sum_{n=0}^{\infty}\frac{1}{\left(2n+1\right)!}\left(rx\right)^{2n+1}\\
& =\underset{x\rightarrow0}{\lim}\sum_{n=0}^{\infty}\frac{r^{2n+1}x^{2n}}{\left(2n+1\right)!}\\
& =\underset{x\rightarrow0}{\lim}\left[r+\frac{r^{3}x^{2}}{6}+...\right]\\
& =r
\end{align}
\end_inset
and likewise for
\begin_inset Formula $\frac{\sin\left(rx\right)}{x}$
\end_inset
.
Thus, there is no issue with convergence, and the matrix exponential along
the line of critical damping
\begin_inset Formula $Q=\frac{1}{2}$
\end_inset
is
\end_layout
\begin_layout Standard
\begin_inset Formula
\begin{equation}
e^{At}=e^{-\Omega t}\left\{ \left(\begin{array}{cc}
1 & 0\\
0 & 1
\end{array}\right)+t\Omega\left(\begin{array}{cc}
-1 & -2\\
\frac{1}{2} & 1
\end{array}\right)\right\}
\end{equation}
\end_inset
\end_layout
\begin_layout Section
Stochastic models of Time-To-Live caching systems
\end_layout
\begin_layout Subsection
Derivation of Poisson distribution from the Poisson process
\begin_inset CommandInset label
LatexCommand label
name "subsec:Derivations-of-Poisson"
\end_inset
\end_layout
\begin_layout Standard
The Poisson process is defined by the transition probabilities
\end_layout
\begin_layout Standard
\begin_inset Formula
\begin{align}
\frac{d}{dt}P_{0}\left(t\right) & =-\lambda_{U}P_{0}\left(t\right)\label{eq:poisson_process0}\\
\frac{d}{dt}P_{k+1}\left(t\right) & =\lambda_{U}P_{k}\left(t\right)-\lambda_{U}P_{k+1}\left(t\right)\,\,k>0\label{eq:poisson_processk}\\
P_{k}\left(0\right) & =\begin{cases}
1 & k=0\\
0 & k>0
\end{cases}
\end{align}
\end_inset
Eq.
(
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:poisson_process0"
plural "false"
caps "false"
noprefix "false"
\end_inset
) has solution
\begin_inset Formula $P_{0}\left(t\right)=e^{-\lambda_{U}t}$
\end_inset
; the probability
\begin_inset Formula $P\left[N_{U}\left(t\right)=0\right]$
\end_inset
decays exponentially with constant rate
\begin_inset Formula $\lambda_{U}$
\end_inset
as time evolves.
More generally, if we assume
\begin_inset Formula
\[
P_{k}\left(t\right)=f\left(t,k\right)e^{-\lambda_{U}t}
\]
\end_inset
for some function
\begin_inset Formula $f$
\end_inset
and substitute into Eq.
(
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:poisson_processk"
plural "false"
caps "false"
noprefix "false"
\end_inset
), we get
\begin_inset Formula
\begin{align}
\frac{d}{dt}\left[f\left(t,k+1\right)e^{-\lambda_{U}t}\right] & =\lambda_{U}\left[f\left(t,k\right)e^{-\lambda_{U}t}\right]-\lambda_{U}\left[f\left(t,k+1\right)e^{-\lambda_{U}t}\right]\\
e^{-\lambda_{U}t}\left[f'\left(t,k+1\right)-\lambda_{U}f\left(t,k+1\right)\right] & =e^{-\lambda_{U}t}\left[\lambda_{U}f\left(t,k\right)-\lambda_{U}f\left(t,k+1\right)\right]\\
f'\left(t,k+1\right) & =\lambda_{U}f\left(t,k\right)
\end{align}
\end_inset
Hence,
\begin_inset Formula $f\left(t,k,\lambda_{U}\right)=C\left(k\right)\left(\lambda_{U}t\right)^{k}$
\end_inset
is polynomial in
\begin_inset Formula $\lambda_{U}t$
\end_inset
so that
\begin_inset Formula
\begin{align}
\frac{d}{dt}\left[C\left(k+1\right)\left(\lambda_{U}t\right)^{k+1}\right] & =\lambda_{U}C\left(k\right)\left(\lambda_{U}t\right)^{k}\\
\left(k+1\right)\lambda_{U}C\left(k+1\right)\left(\lambda_{U}t\right)^{k} & =\lambda_{U}C\left(k\right)\left(\lambda_{U}t\right)^{k}\\
C\left(k+1\right) & =\frac{C\left(k\right)}{\left(k+1\right)}\\
& \Rightarrow C\left(k\right)=\frac{1}{k!}
\end{align}
\end_inset
yielding
\begin_inset Formula
\begin{equation}
P_{k}\left(t\right)=\frac{\left(\lambda_{U}t\right)^{k}e^{-\lambda_{U}t}}{k!}\label{eq:poisson_distribution}
\end{equation}
\end_inset
Eq.
(
\begin_inset CommandInset ref
LatexCommand ref
reference "eq:poisson_distribution"
plural "false"
caps "false"
noprefix "false"
\end_inset
) is the
\emph on
probability mass function
\emph default
for the Poisson distribution.
For constant
\begin_inset Formula $t=T$
\end_inset
, if
\begin_inset Formula
\begin{equation}
P\left[N_{U}\left(T\right)=k\right]=\frac{\left(\lambda_{U}T\right)^{k}e^{-\lambda_{U}T}}{k!}
\end{equation}
\end_inset
then
\begin_inset Formula $N_{U}\left(T\right)$
\end_inset
follows a Poisson distribution
\begin_inset CommandInset citation
LatexCommand cite
key "gardiner2009stochastic"
literal "false"
\end_inset
with rate
\begin_inset Formula $\lambda_{U}T$
\end_inset
; that is,
\begin_inset Formula
\begin{equation}
N_{U}\left(T\right)\sim\text{Poi}\left(\lambda_{U}T\right)
\end{equation}
\end_inset
\end_layout
\begin_layout Subsection
Expected value of Poisson distribution
\begin_inset CommandInset label
LatexCommand label
name "subsec:Exp_value_poisson"
\end_inset
\end_layout
\begin_layout Standard
If
\begin_inset Formula
\begin{equation}
X\sim\text{Poi}\left(\lambda\right)
\end{equation}
\end_inset
then the expected value of
\begin_inset Formula $X$
\end_inset
is given by
\begin_inset Formula
\begin{align}
\left\langle X\right\rangle & =\sum_{k=0}^{\infty}k\frac{\lambda^{k}e^{-\lambda}}{k!}\\
& =\lambda e^{-\lambda}\sum_{k=1}^{\infty}\frac{\lambda^{k-1}}{\left(k-1\right)!}\\
& =\lambda e^{-\lambda}e^{\lambda}\\
& =\lambda
\end{align}
\end_inset
\end_layout
\begin_layout Subsection
Inter-event times of a Poisson process
\begin_inset CommandInset label
LatexCommand label
name "subsec:Inter-event-times"
\end_inset
\end_layout
\begin_layout Standard
The inter-event times of a Poisson process with rate
\begin_inset Formula $\lambda$
\end_inset
are independent
\begin_inset Formula $\text{Exp}\left(\lambda\right)$
\end_inset
random variables.
\end_layout
\begin_layout Paragraph
Proof
\end_layout
\begin_layout Standard
Consider the case where a Poisson process is in the state
\begin_inset Formula $N\left(T\right)=m$
\end_inset
at some time
\begin_inset Formula $T$
\end_inset
.
Let
\begin_inset Formula $t_{0}>0$
\end_inset
be the inter-event time, or the time until the next event occurs; that
is,
\begin_inset Formula $N\left(T+t_{0}\right)=m+1$
\end_inset
and
\begin_inset Formula $N\left(T+t\right)=m$
\end_inset
for
\begin_inset Formula $t