\section{Introduction}
Rare-earth ions are appealing qubits to store and process quantum information. When embedded in crystalline matrices, they combine well protected spin and optically active transitions. The interplay between the spin and atomic like optical excitations places the rare-earth species at the interface between solid-state and atomic physics. Within the lanthanide series, the
Kramers ions (odd number of electrons) exhibit a strong paramagnetic sensitivity because of the crystal field levels admixture leading to a large anisotropic effective electron spin \cite{abragam2012electron}. They also hold a lot of promises in quantum information because a moderate bias magnetic field induces electron paramagnetic resonances (EPR) in the few GHz range. This typical splitting falls in the range of superconducting qubits to form a hybrid interface between the actively investigated quantum electronics circuits and solid-state spins as a memory buffer \cite{grezes, probst}. Additionally, the noticeably narrow inhomogeneous line ($\sim$ 1GHz) is well resolved at low field, thus allowing the spectral optical addressing of spin-selective transitions. The convergence of the appealing spin and optical properties has recently attracted a lot of attention to promote Kramers paramagnetic ions as a quantum spin-photon interface with highly diluted samples \cite{PhysRevLett.113.203601, o2014interfacing, fernandez2015coherent} or more recently with single ions \cite{dibos_isolating_2017, zhong2018optically}. Among them, erbium is emblematic because its optical transition falls in the telecom range and it exhibits large $g$-factor values ($g\sim 15$) in some
crystals such as \eryso or Er$^{3+}$:LiNbO$_3$\xspace \cite{bottger2006spectroscopy, sun2008magnetic,milori_optical_1995}. The tribute to pay to these attractive properties is a rapid loss of decoherence due to the enhanced sensitivity to the magnetic fluctuations. Whatever the viewing angle, optics or EPR, the spin dynamical properties are absolutely critical. They obviously directly impact the outcome of the EPR spectroscopy but they also influence the optical coherence time by the differential magnetic fluctuations between the ground and optically excited states \cite{bottger_decoherence}.

The spin dynamics is essentially governed by two mechanisms \cite{bottger_decoherence}. On the one side, spin lattice relaxations (SLR) describe the interaction of the spin with the phonon bath \cite{Orbach458}. They can be detailed in different mechanisms, namely direct, Raman and Orbach, depending on the process order (one or two-phonon) and the level diagram (on or off-resonance phonon interaction).
However, at liquid helium temperatures, the two-phonon processes are usually negligible. Thus, the SLR is dominated by the direct process. On the other side, the cross-relaxation between adjacent impurities with opposite spin orientations, also known as the flip-flop mechanism (FF), appears as an important process even in low doping samples \cite{van1948dipolar, portis1956spectral}. Both SLR an FF mechanisms primarily describing the population decay, or in other words the longitudinal relaxation, indirectly induce decoherence. Indeed, the fluctuations of the spins generate a background noise leading to decoherence, or transverse relaxation. SLR and FF have very different dependencies as a function of the applied magnetic field \cite{bottger_decoherence, saglamyurek2015efficient, Cruzeiro}.
Practically, SLR dominates at large field and FF at low field. As the magnetic field energy splitting increases and becomes larger than $k_B T$ (where $k_B$ is the Boltzmann constant and $T$ the temperature), the spins become fully polarized, thus preventing the FF. At the same time, the phonon density increases, making the SLR more efficient. At the end, SLR and FF are well separated and can be considered and fought independently to prevent relaxation.
At large magnetic fields, SLR can be controlled by adjusting the level structure with respect to the photon density. This can be done for example by a clever orientation of magnetic field relying on the highly anisotropic Zeeman interaction \cite{bottger2009effects, ranvcic2018coherence} or by a direct nanostructuration of the material to control the intrinsic phonon density\cite{Lutz, yang1999one, yang1999electron}. FF can be also avoided by cooling the spins to few tens of mK thus benefiting from a fully polarized sample at thermal equilibrium. By going deep into the SLR regime, the phonon-bottleneck relaxation has been observed very recently at 20 mK \cite{Budoyo}.

The low field limit at larger temperature, typically 2-4~K range, also deserves consideration because it offers relaxed experimental conditions. As previously mentioned, with a few GHz Zeeman splitting and well resolved optical transitions, this is a region of interest for quantum information applications. The FF dominates the spin dynamics in that case. A partial optical spin polarization can be obtained by frequency selective optical pumping, also known as spectral holeburning \cite{Hastings, saglamyurek2015efficient, Cruzeiro}. The imprinted spectral pattern, usually burned during the preparation stage of quantum memory protocols \cite[and references therein]{heshami2016quantum}, has a lifetime limited by the FF mechanism. Motivated by the potential applications and more generally by the possibility to optically control the out-of-equilibrium ensemble polarization, we optically study the FF cross-relaxation process between erbium spins.

The paper is organized as follows. We first remind the theoretical background to evaluate the FF rate between two identical spins. Even in the case of anisotropic $g$-tensors, we derive an analytical formula for the cross-relaxation lifetime. We then use optical techniques (transient excitation and accumulated spectral hole burning) to measure the anisotropy of the FF rate of Er$^{3+}$\xspace in \yso as a function of the magnetic field orientation. We verify that the lifetime varies by several orders of magnitude as expected from the $g$-tensor theoretical calculation. We finally extend our analysis to another Kramers dopant Nd$^{3+}$\xspace and to higher symmetry crystals (as \cawo and \linbo), which are also investigated for quantum information storage and processing.

\section{Theory}
\label{sec:theory}
The identical spins flip-flop dynamics has been investigated in the context of spin diffusion in a broad sense (spatial and spectral) \cite[chap.4, Spin diffusion in solids]{asakura1998solid}. Starting from the seminal work of Bloembergen \cite{BLOEMBERGEN1949386}, this is still an active subject of research widely stimulated by the perspectives in quantum information for which spin impurities in solids appear as a useful resource \cite[and references]{PhysRevApplied.6.044001}. We will follow the historical perturbative framework to describe the exponential relaxation of the out-of-equilibrium Er$^{3+}$\xspace spin polarization.

\subsection{Flip-flop return to equilibrium}
\label{sec:theory_general}
We assume the $N$ spin-$1/2$ identical rare-earth impurity ions randomly substituted in host crystalline matrix. The energy level degeneracy is lifted by an external magnetic field $\vec{B}$. The resulting Zeeman splitting is set much smaller than the thermal energy $k_BT$ and, although weakly coupled to the heat bath, the spins are evenly distributed over the $\ket{+}$ and $\ket{-}$ eigenstates of the Zeeman Hamiltonian. The spins are coupled to each other by magnetic dipole-dipole interaction. The interaction strength is assumed to be much smaller than the Zeeman splitting, which only allows flip-flop transitions, where two spins simultaneously flip in opposite ways, with no total energy change.

In a spectral hole burning experiment, a subset $W_+$, containing $N_0<<N$ impurity ions, can be prepared in state $\ket{+}$ at time $t=0$. That population imbalance is destroyed by FF interactions between the optically probed ions and the rest of the ensemble. The corresponding exponential-decay rate $R$ is given by the Fermi's Golden Rule \cite[W.B. Mims, Chap.4, p.294]{geschwind1972electron}:
\begin{equation}
R =  \frac{2 \pi }{ \hbar}\left\{\sum_{j\in\left\{\ket{-}\right\}}\left|\bra{+-}\bf{H}_{ij}\ket{-+}\right|^2\right\}_{i\in W_+} \frac{1}{\hbar\Gamma^{\rm spin}_{\rm inh}} ,
\label{eq:FGR_rate}
\end{equation}
where $\Gamma^{\rm spin}_{\rm inh}$ represents the inhomogeneous broadening of the spin transition, in $s^{-1}$, $\left\{\right\}_{i\in W_+}$ expresses averaging over the subset $W_+$, and $\bf{H}_{ij}$ stands for the magnetic dipole-dipole Hamiltonian of the $i,j$ ion pair:
\begin{equation}
\bf{H}_{ij} = - \frac{\mu_0}{4\pi} \frac{1}{r_{ij}^3}  \left[ 3 \left( \vec{\mu}_{i} \cdot \vec{u}_{ij} \right) \left( \vec{\mu}_{j} \cdot \vec{u}_{ij} \right)-\vec{\mu}_{i} \cdot \vec{\mu}_{j} \right]  \;,
\label{eq:Hint}
\end{equation}
where $\vec{r_{ij}}$, $\vec{u}_{ij}$, and $\vec{\mu}_{i,j}$ respectively represent the inter-ion vector, the unit vector along $\vec{r_{ij}}$, and the magnetic dipole moments.

The $\bf{H}_{ij}$ structure suggests that the sum over $j$ in Eq.~\ref{eq:FGR_rate} can be approximately evaluated as a combination of angular averaging over $\vec{u}_{ij}$ direction and sum over the inter-ion distance. Therefore:
\begin{multline}
\sum_{j\in\left\{\ket{-}\right\}}\left|\bra{+-}\bf{H}_{ij}\ket{-+}\right|^2\approx\\
\left(\frac{\mu_0}{4\pi}\mu_{\rm B}^2\right)^2\times\Xi\left(\bar{\bar{g}},\vec{B}\right)\times\sum_{j\in\left\{\ket{-}\right\}} \frac{1}{r_{ij}^6},
\label{eq:sum_over_ions}
\end{multline}
where $\mu_{\rm B}$ denotes the Bohr magneton, and the dimension-less coupling factor $\Xi\left(\bar{\bar{g}},\vec{B}\right)$, defined as:
\begin{multline}
\Xi\left(\bar{\bar{g}},\vec{B}\right)=\frac{1}{4\pi}\times\mu_{\rm B}^{-4}\times\\
\int d\vec{u}_{ij}\left|\bra{+-} 3 \left(\vec{\mu}_{i}\cdot\vec{u}_{ij}\right) \left(\vec{\mu}_{j} \cdot \vec{u}_{ij} \right)-\vec{\mu}_{i} \cdot \vec{\mu}_{j} \ket{-+}\right|^2,
\label{eq:Xi_factor}
\end{multline}
does not refer to a specific interacting ion pair but depends on the $\bar{\bar{g}}$-tensor, and on $\vec{B}$ direction, which determines the eigenvectors $\left(\ket{+},\ket{-}\right)$. The detailed calculation of $\Xi\left(\bar{\bar{g}},\vec{B}\right)$ is deferred to Appendix~\ref{sec:Vanalytical}.

To evaluate the sum over $1/r_{ij}^6$ in Eq.~\ref{eq:sum_over_ions}, we replace the impurity-ion random distribution by an evenly spaced arrangement of $j$ spins in state $\ket{-}$ on a fictitious cubic lattice. Each $i$ spin in state $\ket{+}$ is thus surrounded by layers of those $j$ spins (see Appendix \ref{sec:sum_distance}).  Since the weighting factor $1/r_{ij}^6 $ rapidly drops with distance, one may truncate the sum to a few layers. With truncation to the first 10 layers, the sum reduces to $8.4n_s^2$, where $n_s$ represents the spatial density of spins in state $\ket{-}$, which coincides with half the impurity-ion density. Finally, the FF decay rate reads as:
\begin{equation}
R = \frac{2 \pi }{ \hbar}\left(\frac{\mu_0}{4\pi}\:n_s\:\mu_{\rm B}^2\right)^2\times 8.4 \times\Xi\left(\bar{\bar{g}},\vec{B}\right)\times\frac{1}{\hbar\Gamma^{\rm spin}_{\rm inh}} ,
\label{eq:FGR_rate_final}
\end{equation}

A magnetic dipole moment $\vec{\mu}_i$ can be expressed in terms of the effective spin operator $\vec{\hat{S}}_i$ and of the $\bar{\bar{g}}$-tensor as:
\begin{equation}
\vec{\mu}_i=\mu_{\rm B}\bar{\bar{g}}\cdot\vec{\hat{S}}_i
\label{eq:mag_dipole}
\end{equation}
In the frame $(x,y,z)$ where the $\bar{\bar{g}}$-tensor is diagonal, $\vec{\mu}_i$ reads as:
\begin{equation}
\vec{\mu}_i=\mu_{\rm B}\left(\hat{S}_{i,x}g_x\vec{u}_x+\hat{S}_{i,y}g_y\vec{u}_y+\hat{S}_{i,z}g_z\vec{u}_z\right)
\label{eq:mag_dipole_expansion}
\end{equation}
where:
\begin{equation}
\bar{\bar{g}} = \left( \begin{matrix} g_x & 0 & 0 \\ 0 & g_y & 0 \\ 0 & 0 & g_z \end{matrix} \right) \;
\end{equation}
and $(\vec{u}_x,\vec{u}_y,\vec{u}_z)$ represent the unit vectors of the reference frame.
When $\vec{B}$ is directed along $Oz$, $\hat{S}_z$ and the Zeeman Hamiltonian share the same  eigenvectors. Hence, $\bra{+}\hat{S}_{i,z}\ket{-}$ vanishes and $\Xi\left(\bar{\bar{g}},\vec{B}\right)$ does not depend on $g_z$. Conversely, $\Xi\left(\bar{\bar{g}},\vec{B}\right)$ does not depend on $g_x$ or $g_y$ when $\vec{B}$ is directed along $Ox$ or $Oy$ respectively. That feature will prove to be important in the following.

\subsection{Application to \eryso}
\label{sec:theory_eryso}

In order to apply the above model to \eryso, we only keep the zero nuclear spin isotopes (78\% of the doping ions), ignoring the $^{167}$Er isotope (22\% abundance with non-zero nuclear spin). The required smallness of the spin-spin interaction with respect to Zeeman splitting is satisfied at external magnetic field values as small as $1$~mT. In a 10 ppm \eryso crystal, at $B = 1$~mT, the Zeeman interaction varies between 5 and 50 MHz depending on the crystal orientation, whereas the Er-Er interaction ranges within 1 and 100 kHz.
As we will discuss in \ref{discussion}, we use a very weak magnetic field to minimize the inhomongeous spin broadening. We choose to keep it constant to $B = 0.3$~mT and vary the field orientation to reveal the anisotropy.
The example of \eryso is particularly interesting because it shows strongly anisotropic properties. Indeed, in the frame $(x,y,z)$ where the $\bar{\bar{g}}$-tensor is diagonal,
\begin{equation}
\; g_z \approx 10 \, g_y \approx 30 \, g_x \; .
\end{equation}

As pointed out above, $ \Xi\left(\bar{\bar{g}},\vec{B}\right)$ only depends on the $g_x$ and $g_y$ components, much smaller than $g_z$ in \eryso, when the applied magnetic field is directed along $Oz$. According to the general expression, derived in Appendix~\ref{sec:Vanalytical}:
\begin{equation}
 \Xi\left(\bar{\bar{g}},\vec{B}\right) = \frac{1}{20}   \left( g_x^4 + g_y^4 - g_x^2 g_y^2 \right) \; .
 \label{eq:Vmoy_130deg}
 \end{equation}
This situation is almost reached when $\vec{B}$ is directed at about $\phi = 135^\circ$ from $D_1$ in the $D_1-D_2$ plane, or in other words, the principal axis almost lies in the $D_1-D_2$ plane

The $g_z$ contribution raises as $\vec{B}$ direction departs from $Oz$, reaching a maximum when $\vec{B}$ lies in the $xOy$ plane. Then $\Xi\left(\bar{\bar{g}},\vec{B}\right)$ reads as:
\begin{equation}
  \Xi\left(\bar{\bar{g}},\vec{B}\right) =
\frac{1}{20}\left(g_z^4+g_\bot^4-g_z^2g_\bot^2\right)
  \label{eq:Vmoy_30deg}
  \end{equation}
where:
\begin{equation}
\frac{1}{g_\bot^2}=\frac{1}{B^2}\left(\frac{B_x^2}{g_y^2}+\frac{B_y^2}{g_x^2}\right)
\; .
\end{equation}

This situation occurs for $\phi \approx 30^\circ$. According to Eqs.\ref{eq:Vmoy_130deg} and \ref{eq:Vmoy_30deg}, the relaxation rate $R$ should vary by a factor of $\left( g_z / g_y \right)^4 \simeq 10^4$ between those two situations (neglecting the small $g_x$ contribution). Note that in a scenario where the $g$-tensor is isotropic, $\Xi\left(\bar{\bar{g}},\vec{B}\right)$ would be totally independent of the external magnetic field direction.

We calculate $\Xi\left(\bar{\bar{g}},\vec{B}\right)$ variations with the external magnetic field $\vec{B}$ direction in the crystalline frame $(D_1,D_2,b)$. Fig.~\ref{fig:T1_map_10ppm_ErYSO} shows the map of $\displaystyle T_{ff}=1/R$ following \eqref{eq:FGR_rate_final} as a function of $\phi$ and $\theta$ the angles of $\vec{B}$ in the $(D_1,D_2,b)$ frame. We observe strong variations within several planes, including the $D_1-D_2$ plane that we considered above.

\begin{figure}[t]
\centering
\includegraphics[width=0.9\columnwidth]{T1_map_10ppm_ErYSO_v3.pdf}
\caption{Map of $\displaystyle T_{ff}=1/R$ as a function of the orientation of the external magnetic field in \eryso (10 ppm). The angles $\phi$ and $\theta$ are the angles of the external magnetic field $\vec{B}$ in the crystalline frame $(D_1,D_2,b)$. When $\theta = 90^\circ$, the field lies in the $D_1-D_2$ plane. Note the logarithmic colormap.}
\label{fig:T1_map_10ppm_ErYSO}
\end{figure}

\section{flip-flop rate in erbium doped \yso}

Because of the strong variation of the expected lifetime, different techniques must be used to cover the measurement range. We first discuss the methodology and then detail the experimental results.

\subsection{Spin lifetime optical measurement techniques}\label{techniques}

The method of choice to measure optically the spin lifetime is spectral hole burning (SHB). With modern laser stabilization techniques, one can reach SHB spectral resolution as small as the optical homogeneous linewidth, which gives access to structures below one MHz splitting. Zeeman and hyperfine structure for both Kramers and non-Kramers ions have been revealed and studied in pioneering works \cite{Macfarlane:81, PhysRevB.38.11061}. Stimulated by the recent development of diode lasers and the contextual interest in quantum information for which rare-earth crystals represent a promising quantum memory support, SHB has reappeared as a versatile tool to study the population dynamics and prepare the memory in the initial step of many storage protocols \cite{Seze_PhysRevB.73.085112, Hastings-Simon_PhysRevB.77.125111, Hastings, AFZELIUS20101566, Lauritzen_PhysRevB.85.115111, Cruzeiro}.

\begin{figure}[!ht]
\centering
\includegraphics[width=0.9\columnwidth]{shb_vs_pi.pdf}
\caption{Optical measurement techniques. Left - Accumulated spectral holeburning: the spin polarization is obtained by pumping through the optically excited state. Right - Optical inversion recovery: A spin population imbalance can also be created by transient optical excitation ($\pi$-pulse typically) when the spin lifetime is much shorter than the optical decay times. See text for details.}
\label{fig:shb_vs_pi}
\end{figure}

In a simplified picture, SHB operates as optical pumping in an equivalent three-level system, as represented in Fig.\ref{fig:shb_vs_pi} (left). This ideal scheme corresponds well to the case of \eryso. The $^{4} I_{15/2} \rightarrow \, ^{4}I_{13/2}$ optical transition of Er$^{3+}$\xspace splits into ground and excited Kramers doublets leading to a 4-level structure. The laser frequency determines the spin levels that are addressed by the optical excitation ($\ket{-}_g \rightarrow \ket{-}_e$).  We can reduce the structure to a 3-level by just neglecting $\ket{+}_e$. This simplification is made possible because the direct spin relaxation in the excited state $\ket{-}_e \leftrightarrow \ket{+}_e$ can be neglected. More precisely, in the FF regime, the optically excited spin density is so weak that the  $\ket{-}_e \leftrightarrow \ket{+}_e$ cross-relaxation is extremely unlikely \cite{welinski2018electron, rakonjac2018spin}. The SHB structure consists of 3 levels, namely $\ket{-}_g$, $\ket{-}_e$ and $\ket{+}_g$. Instead, if the laser excites the crossed transition $\ket{-}_g \rightarrow \ket{+}_e$, which corresponds to a different subset of ions (frequency class) \cite{nilsson2004hole}, the $\ket{-}_e$ level can be neglected.

In this archetypal SHB 3-level structure, the continuous optical excitation is redistributed between the ground state sublevels, in accordance with the branching ratio $\mathcal{R}$ from the excited state. The hole decay reflects the spin relaxation. Ideally, the cumulative SHB procedure can fully polarize the spins, provided spin relaxation is slow as compared to the optical pumping rate $\mathcal{R}/T_1^\mathrm{opt}$ where $T_1^\mathrm{opt}$ is the optical lifetime. This condition may not be satisfied in Kramers ions. This is precisely the case for Er$^{3+}$\xspace. As a consequence, the spin state preparation can be quite challenging \cite{Lauritzen_PhysRevA.78.043402}.

As we discussed in \ref{sec:theory_eryso}, the expected FF decay time varies by four orders of magnitude, and can be faster or slower than optical pumping, which is always limited by the exceptionally long value of $T_1^\mathrm{opt}=11$~ms among the rare-earth elements. When spin relaxation is too fast, one cannot accumulate spin population imbalance, and an alternative to SHB is needed. Instead of using a continuous laser for optical pumping, we perform transient optical excitation as illustrated in \figref{fig:shb_vs_pi} (right). The $\pi$-pulse achieves transient spin imbalance between $\ket{-}_g$ and $\ket{+}_g$ by promoting the $\ket{-}_g$ population into the optical excited state. The rapid recovery, caused by spin cross-relaxation, can be probed on the optical $\ket{-}_g \rightarrow \ket{-}_e$ transition. Having no time to scan the probe frequency, as we did in SHB, we instead only measure the transmission of a weak and short probe pulse. The transmission decay gives the spin lifetime. These two complementary techniques, namely accumulated SHB and optical inversion recovery ($\pi$-pulse excitation) as distinguished in \figref{fig:shb_vs_pi}, are used alternatively depending whether the expected FF lifetime is shorter or larger than $T_1^\mathrm{opt}=11$~ms. Rigorously, the relevant parameter to compare with the FF rate is the pumping rate $\mathcal{R}/T_1^\mathrm{opt}$. Nonetheless, the branching ratio  $\mathcal{R}$ also exhibits a strong magnetic field orientation dependency \cite{lauritzen}. So for simplicity, we have decided to keep $T_1^\mathrm{opt}=11$~ms as a  reference value to distinguish the different dynamical regimes.

In practice, as calculated for a 10 ppm \eryso sample (see section \ref{sec:theory_eryso}), the FF lifetime depends strongly on the magnetic field angle $\phi$ in the $D_1$-$D_2$ plane and varies between $30$~$\mu$s and $200$~ms. For $\phi \leq 70^\circ$, the FF rate is fast and so we use the optical inversion recovery measurement technique (blue shaded area in \figref{fig:Tz10ppm}). Conversely, SHB can be accumulated for $\phi > 70^\circ$ and longer spin lifetimes can be extracted from the decay curve (yellow shaded area in \figref{fig:Tz10ppm}).

\subsection{Experimental results}\label{experiment}

The experimental validation of the analytical formula \eqref{eq:FGR_rate_final} is quite challenging because of the four orders of magnitude covered by the FF rates and therefore, the spin lifetimes we need to measure. We also propose to investigate different doping concentrations to validate our formula.

The two crystals under study are monoclinic \yso, grown by Scientific Materials Corporation. Erbium substitutes for yttrium at dopant concentrations of respectively $10$ ppm and $50$ ppm. Among the two crystallographic substitution sites, we consider the one at $1536.48$ nm, referenced as site 1\cite{bottger2006spectroscopy}. We use a liquid helium cryostat to cool down the crystal at $2$ K and a superconducting coil, monitored by a low voltage supply, to generate a weak $0.3$ mT magnetic field. We use the 3 optical extinction axis ($D_1$-$D_2$-$b$) as a reference frame. Light propagates along $b$ whereas the crystal is placed on a rotating mount (Attocube ANRv51/LT) so that the external magnetic field can be rotated in the whole $D_1$-$D_2$ plane.

As discussed in \ref{techniques}, to measure the spin lifetimes of the $10$ ppm crystal, we implement the optical inversion recovery technique for angles $\phi$ from $5^\circ$ to $60^\circ$, where spin lifetimes are expected to be shorter than the optical lifetime, and a SHB experiment from $80^\circ$ to $180^\circ$ for longer lifetimes. For the $50$ ppm sample, as the rate is proportional to the density squared [see \eqref{eq:FGR_rate_final}], the lifetimes are expected to be shorter than $T_1^\mathrm{opt}=11$~ms at any angle. Therefore, the optical inversion recovery is implemented exclusively.

For both measurement techniques, the optical setup is the same. We split a Koheras fiber laser in two beams, the pump and the probe. The pump is amplified with a ManLight erbium-doped fiber amplifier. These two beams are temporally shaped by acousto-optic modulators, controlled by an arbitrary waveform generator (Tektronix AWG520). Then, pump and probe beams counterpropagate in the crystal with typical power of $20$ mW and $100$ $\mu$W respectively. Finally, the probe is measured by an avalanche photodetector.
The time sequences are different for the two techniques. For SHB experiment, we use a $100$ ms pump and $500$ $\mu$s weak probe. We sweep the probe frequency to record the central hole area decays. For the optical inversion recovery, a $200$ ns rms-duration pulse has a typical area of $\pi$ (transient excitation, see \figref{fig:shb_vs_pi}). The mean transmission of a $10$ $\mu$s probe pulse gives us the temporal evolution of the absorption at different waiting times after the $\pi$-pulse (inversion recovery). \figref{fig:2exp10ppm} gives two examples of the technique for the 10 ppm amd 50 ppm samples respectively.

\begin{figure}[!t]
\centering
\includegraphics[width=0.9\columnwidth]{decayAll_1.pdf}
\caption{a : Optical inversion recovery decay of the 10 ppm sample for $\vec{B}$ at $45^\circ$. We fit data with $2$ exponential decay curves ($N=1$). b : Optical inversion recovery decay of the 50 ppm sample for $\vec{B}$ at $115^\circ$. We fit data with 2 or 3 exponential decay curves (respectively black for $N=1$ and red for $N=2$).}
\label{fig:2exp10ppm}
\end{figure}

Both techniques exhibit decay curves characterizing the population dynamics.
In a simplified $3$-level system, the exponential decays should have two characteristic times (see \figref{fig:2exp10ppm}), $T_1^\mathrm{opt}$ and $\tau$ where $T_1^\mathrm{opt}$ is the optical lifetime and $\tau$ the spin lifetime to be compared to $T_{ff}$. We should be able to extract the spin lifetime by fitting data with the following formula
\begin{equation}\label{eq_exp1}
\Delta_{\alpha L}(t)=\alpha L((1-a)e^{-\frac{t}{T_1^\mathrm{opt}}}+ae^{-\frac{t}{\tau}})
\end{equation}
where $T_1^\mathrm{opt}$ is fixed to $11$ ms. We choose to let the coefficients $\alpha L$ and $a$ as free parameters even thought they could be measured and calculated respectively. We actually observe that their values depend on the efficiency of the $\pi$-pulse excitation which varies significantly when rotating the sample.

However in practice, \eqref{eq_exp1} is not sufficient to fit the very different situtations encountered when we change the magnetic field and the concentration. To obtain a good agreement between the data and the fit, we use a multi exponential decay formula as
\begin{equation}\label{eq_expN}
\Delta_{\alpha L}(t)=\alpha L\left((1-\sum\limits_{i=1}^Na_i)e^{-\frac{t}{T_1^\mathrm{opt}}}+\sum\limits_{i=1}^Na_ie^{-\frac{t}{\tau_i}}\right)
\end{equation}
and increase $N$ until $\chi^2 < 10^{-3}$ (good agreement threshold).
For the 10 ppm crystal, we reach a good agreement ($\chi^2 < 10^{-3}$) with $N=1$ for the optical inversion recovery technique (blue shaded area in \figref{fig:Tz10ppm}) as illustrated by the example in \figref{fig:2exp10ppm}). For the same sample using the SHB technique, it is necessary to increase to  $N=3$ to properly fit data, thus revealing longer timescale (typically from $1$ s to $100$ s) in this configuration as we will discuss in \ref{discussion}. For the 50 ppm sample, $N=2$ is sufficient for the different angles using optical inversion recovery as illustrated in \figref{fig:2exp10ppm}.b.

To evaluate the FF rates, we extract the characteristic time $\tau_1$ given by the dominant exponential term where $a_1$ is the largest $a_i$. The measurements for a varying magnetic field are represented in \figref{fig:Tz10ppm} and \figref{fig:Tz50ppm} for 10 ppm and 50 ppm respectively and compared with the theoretical prediction.

\begin{figure}[!ht]
\centering
\includegraphics[width=0.9\columnwidth]{Tz10ppm.pdf}
\caption{Variation of the spin lifetime with the angle $\phi$ of the external magnetic field in the $D_1 - D_2$ plane for the $10$ ppm crystal. Blue shaded area, optical inversion recovery technique when the expected lifetime is shorter than $T_1^\mathrm{opt}=11$~ms (green dashed line). Yellow shaded area, accumulated SHB technique for longer lifetimes. The markers are the experimental measurements of $\tau_1$ and the red curves the theoretical calculation of $T_{ff}$ (see text for details)}
\label{fig:Tz10ppm}
\end{figure}

\begin{figure}[!ht]
\centering
\includegraphics[width=0.9\columnwidth]{Tz50ppm.pdf}
\caption{Variation of the spin lifetime with the angle $\phi$ of the external magnetic field in the $D_1 - D_2$ plane for the $50$ ppm crystal. The experimental measurements corresponds to the optical inversion recovery technique only. We cannot measure lifetime shorter than $10$ $\mu$s (represented as a black dashed line) corresponding to our probe pulse duration.}
\label{fig:Tz50ppm}
\end{figure}

For the theoretical calculation, we plot a confidence interval (light-red area in Figs.\ref{fig:Tz10ppm} and \ref{fig:Tz50ppm}) corresponding to our uncertainty in the spin inhomogeneous broadening  $\Gamma_\mathrm{inh}^\mathrm{spin}$, which appears in \eqref{eq:FGR_rate_final}. $\Gamma_\mathrm{inh}^\mathrm{spin}$ varies with the magnetic field strength and orientation as will be discussed in section \ref{discussion}. We measure it for a $0.3$ mT external fields and with different orientations via SHB technique. The width of the associated antihole (positionned at the excited state splitting) corresponds to the inhomogeneous broadening of spin levels \cite{Lauritzen_PhysRevA.78.043402}. The red solid curves bounding the confidence interval have been calculated with the extremal values we observe for $\Gamma_\mathrm{inh}^\mathrm{spin}$ ($\Gamma_\mathrm{inh}^\mathrm{spin}\in \left[2.3,6.3\right]$ MHz for $10$ ppm and $\Gamma_\mathrm{inh}^\mathrm{spin}\in \left[2.8,6\right]$ MHz for $50$ ppm).

\subsection{Discussion}\label{discussion}

The agreement between the experimental data and the theoretical prediction is overall satisfying. First, we observe a strong dependence of the spin lifetime at very low magnetic field, reflecting the anisotropic character of the FF rate, as predicted by our theoretical model. From our measurements in the 10 ppm crystal, the anisotropy factor on the spin lifetime is $1.2 \times 10^3$ ($1.4 \times 10^4$ theoretically) as the magnetic field is rotated in the $D_1 - D_2$ plane.
When the lifetimes are minimum (for $\phi \approx 30^{\circ}$), the experimental values are about 10 times larger than the predictions (fact which is enhanced for the short lifetimes by the logarithmic scale). This could be explained by an uncertainty on the inhomogeneous broadening of the spin transition for instance, as discussed below. In the case of the 50 ppm crystal, the agreement is less satisfying because the ratio between the maximum and the minimum lifetimes is only $10^2$ instead of $10^4$. As we will discuss below, we cannot measure lifetimes shorter than 10 $\mu$s, still the exponential fit converges to well defined decay times. It is nevertheless difficult to track down the different decays because of the shorter orders of magnitude and the technical lower 10 $\mu$s limit. It should be also noted that this angular region ($\phi \approx 30^{\circ}$) is particularly singular because of the low $g$-factor values for ground and excited states \cite{sun2008magnetic}. In other words, the optically probed transitions $\ket{-}_g \rightarrow \ket{-}_e$ and $\ket{+}_g \rightarrow \ket{+}_e$ (see fig.\ref{fig:shb_vs_pi}) have a comparable frequency. So a mixing between different mechanisms is possible. Precautions in the analysis are needed in this angular region.
However, the overall variation caused by the change of concentration of a factor of $5^2 = 25$ is clearly observed. This validates the fact that we are in presence of a concentration dependent phenomenon, i.e. a process involving two distinct Er$^{3+}$ ions. The FF process indeed theoretically shows a quadratic dependence with the concentration. In the following we expose our experimental sources of uncertainties, as well as fundamental reasons that could explain the differences between experimental and theoretical results.

One of the main experimental sources of error is the uncertainty in measuring the spin inhomogeneous broadening $\Gamma_\mathrm{inh}^\mathrm{spin}$, which induces a significant spread of the theoretical values. The exact modeling of the spin broadening mechanism is beyond the scope of this paper but it should be noted that $\Gamma_\mathrm{inh}^\mathrm{spin}$ exhibits a specific dependency when the magnetic field is rotated. This latter can be evaluated quantitatively from the $g$-tensor \cite{Welinski2016} and would induce an extra angular dependency of the FF rate. In other words, both terms in \eqref{eq:FGR_rate_final}, $\Xi\left(\bar{\bar{g}},\vec{B}\right)$ and $\Gamma_\mathrm{inh}^\mathrm{spin}$, would depend on the magnetic field orientation. For the sake of simplicity, we plot the rate $R$ for the two measured extremal values of $\Gamma_\mathrm{inh}^\mathrm{spin}$ and keep an angular dependency of $\Xi\left(\bar{\bar{g}},\vec{B}\right)$ solely. This gives the red solid curves bounding the confidence interval in Figures~\ref{fig:Tz10ppm} and \ref{fig:Tz50ppm} (light red area). Without going much deeper into the analysis, it should be nevertheless noted than we observe at 0.3~mT approximately the same spin linewidths ($\sim 5$~MHz) than the ones measured via EPR at 60~mT \cite{Welinski2016}. Assuming a linear broadening with magnetic field, this difference corresponds to a factor of 15 in the broadening coefficients. This questions the origin of the inhomogeneous broadening, which may vary from low to high magnetic fields at different orientations. Even if we intentionally use a very weak magnetic to reduce the spin inhomogeneous broadening and measure it conveniently with the SHB technique, further investigations are needed to fully understand its origin.

Another experimental limitation imposed by our measurement techniques is the shortest measurable spin lifetimes of about 10 $\mu$s. In the optical inversion recovery technique (see Section~\ref{techniques}), we use a strong $\pi$-pulse with a duration of less than 1 $\mu$s. However, this strong pulse blinds the photo-detector for approximately 10 $\mu$s. Then, we send a 10 $\mu$s probe pulse whose integration gives the average absorption. As a consequence, the shortest spin lifetime that we can measure is 10 $\mu$s whereas the shortest expected lifetime for the 50 ppm crystal is 1 $\mu$s.

To discuss the possible influences that are not included in our model, we can first mention the presence of Er$^{3+}$\xspace spins in the crystallographic site 2. They should not modify the FF rate of Er$^{3+}$\xspace ions in site 1 that we probe optically. Indeed, both sites have very different spin transitions (different $g$-factors), so cross-relaxation between site 1 and 2 is off-resonant and therefore very unlikely. In other words, site 1 and 2 should behave as two independent groups. This distinction is generally true except if the sites have the same effective $g$-factor. This singular situation occurs close to $\phi \sim 135 ^\circ$ \cite{sun2008magnetic}. In that case, cross-relaxation between sites is resonant and may happen as soon as the spin transition difference (between sites) is smaller than $\Gamma_\mathrm{inh}^\mathrm{spin}$. This corresponds to a few degrees around $135 ^\circ$ and may locally increase the FF rate because the density of resonantly interacting ions is doubled. Around this region, two points of \figref{fig:Tz10ppm} seem to have a lower $T_{ff}$ value and confirm this interaction. This would deserve more investigations near $\phi \sim 135 ^\circ$.

Another source of magnetic interaction is given by the so-called superhyperfine interaction \cite{guillot-noel_direct_2007} with the ligand (yttrium in that case) that we neglect in a first approach. The electro-nuclear spin coupling makes the electronic spin description significantly more complex, as recently discussed \cite{Car18}. Despite the weakness of the interaction between Er$^{3+}$\xspace and the Y$^{3+}$\xspace nuclear spin, the Er$^{3+}$\xspace-Y$^{3+}$\xspace proximity makes the superhyperfine interaction comparable with the Er$^{3+}$\xspace-Er$^{3+}$\xspace electronic coupling in diluted samples (about 100~kHz in a 10 ppm crystal). Before qualitatively incriminating any interaction that is not included in the model as the responsible of the partial disagreement with the theoretical prediction, it should be kept in mind that these interactions, including the superhyperfine coupling, are weak compared to the spin inhomogeneous broadening $\Gamma_\mathrm{inh}^\mathrm{spin}$. So they should not drastically influence the FF dynamics even if the exact coupling between Er$^{3+}$\xspace-Y$^{3+}$\xspace electro-nuclear mixed states still remains to be evaluated, first in a perturbative approach.

Finally, we wanted to discuss the presence of multiple characteristic decay times for the spin relaxation. As described in section \ref{experiment}, depending on the concentration and the field orientation, the population dynamics cannot be described by two timescales only, namely the optical radiative relaxation and the FF rate. We also recognize long lived structures as observed early in similar conditions by Hastings-Simon {\it at al.} \cite{Hastings}. Even if these extra mechanisms are not dominant, in the sense that they represent a small fraction of the ions, they still reveal a complex dynamics. In our SHB measurements, those long decays are visible on the anti-holes dynamics (enhanced absorption shifted by $\Delta_g - \Delta_e$ where $\Delta_g$ and $\Delta_e$ are the ground and excited state Zeeman splitting respectively). This means that the ions with the same Zeeman coefficient and sensing the same magnetic field can relax differently. This excludes for instance an explanation based on the $^{167}$Er isotope, which represents 22\% of the ions and possesses a nuclear spin. They may exhibit slower FF rates because of their much weaker effective concentration, not only because of the 22\% abundance, but also because their population is distributed among different hyperfine states making the resonant FF process much less likely. However, because of the hyperfine coupling, the shift of $\Delta_g - \Delta_e$ will be different than for the other isotopes, so $^{167}$Er do not participate to the anti-holes that we observed. Nevertheless, a comparable analysis with a 10 ppm or 50 ppm isotopically pure $^{167}$Er crystal would be interesting to explore the FF rate in the presence of hyperfine coupling.

\section{Extrapolation to other materials}

After an experimental validation in \eryso, our study can be extended to other combinations of dopants and hosts to predict the orders of magnitude of the FF rate and its angular dependence (magnetic field orientation). In this section, we compare the spin lifetime limited by FF in other materials that are investigated  in quantum information namely \ercawo \cite{bertaina_rare-earth_2007}, Er$^{3+}$:LiNbO$_3$\xspace \cite{0953-4075-45-12-124013, Hastings-Simon2006} and \ndyso \cite{Cruzeiro}.

We first consider Er$^{3+}$\xspace -doped materials to reveal the influence of the $g$-tensor anisotropy. The comparison between \ercawo and Er$^{3+}$:LiNbO$_3$\xspace is particularly enlightening in that sense. For both, the Er$^{3+}$\xspace site has a higher symmetry than in \yso.
We also analyze \ndyso to compare two Kramers ions (Er$^{3+}$\xspace and Nd$^{3+}$\xspace) in the anisotropic low symmetry \yso matrix. Finally we present the spin lifetime dependency as a function of the doping concentration for all these materials.

\subsection{Er$^{3+}$\xspace ions in different host matrices}

In the previous section, we studied experimentally the FF rate in \eryso and showed its strong anisotropic behavior, as predicted by the theoretical model. \eryso is a very interesting material because of its strong anisotropic properties and long optical coherence lifetimes \cite{bottger_optical_2006}, but other host matrices are good alternatives if the FF is a limiting factor. We here focus on \ercawo and Er$^{3+}$:LiNbO$_3$\xspace, which both present simpler crystalline structures and higher symmetries (S$_4$ and C$_3$ respectively for the substitution site). \ercawo exhibits long spin coherence time \cite{mims_phase_1968}. Er$^{3+}$:LiNbO$_3$\xspace would benefit from all the knowledge in opto-electronics where the \linbo matrix is a reference photonic platform. As already mentioned, both materials find application in quantum information \cite{bertaina_rare-earth_2007, 0953-4075-45-12-124013, Hastings-Simon2006}.

\begin{figure}[t!]
\centering
\includegraphics[width=0.95\columnwidth]{TzCWO_LNO.pdf}
\caption{Dependence of the spin lifetime on the orientation of the external magnetic field for \ercawo (red line) and Er$^{3+}$:LiNbO$_3$\xspace (blue line) crystals. We choose a concentration of 10 ppm and a inhomogeneous spin broadening $\Gamma_{\rm inh}^{\rm spin}$ of 5 MHz in both cases.}
\label{fig:Tz_10ppm_CaWO4_LiNbO3}
\end{figure}

\ercawo and Er$^{3+}$:LiNbO$_3$\xspace show a diagonal $g$-tensor in the crystalline frame $(a,b,c)$ where some properties are invariant by rotation around the $c$-axis ($a$ and $b$ are equivalent). Thus, the $g$-tensor has the following form
\begin{equation}
\bar{\bar{g}} = \left( \begin{matrix} g_\perp & 0 & 0 \\ 0 & g_\perp & 0 \\ 0 & 0 & g_\parallel \end{matrix} \right) \; .
\end{equation}
The main difference between the two materials is that $g_\perp > g_\parallel$ for \ercawo, whereas $g_\parallel > g_\perp$ for Er$^{3+}$:LiNbO$_3$\xspace. This distinction has a major impact on the FF mechanism.

More precisely, the coupling term $\Xi\left(\bar{\bar{g}},\vec{B}\right)$ is almost isotropic in the case of \ercawo, i.e. almost constant when we rotate the magnetic field in the $(a,c)$ plane ($\theta$ is the angle between $\vec{B}$ and the $c$-axis). Indeed, there is always one perpendicular component that is large because $g_\perp > g_\parallel$ (see section~\ref{sec:theory_eryso}) thus the coupling $\Xi\left(\bar{\bar{g}},\vec{B}\right)$ is dominated by $g_\perp^4$ in any case.

Conversely, Er$^{3+}$:LiNbO$_3$\xspace shows a very anisotropic $\Xi\left(\bar{\bar{g}},\vec{B}\right)$ because when $\vec{B} \parallel c$, the two perpendicular components are negligible. Figure~\ref{fig:Tz_10ppm_CaWO4_LiNbO3} shows the spin lifetime in these two materials as a function of the magnetic field angle $\theta$ from the $c$-axis for an Er$^{3+}$\xspace concentration of 10 ppm. We clearly observe the quasi isotropic behavior of FF rate in \ercawo, whereas Er$^{3+}$:LiNbO$_3$\xspace is strongly anisotropic. This leads to a much longer optimum spin lifetime (for $\theta = 0$) in Er$^{3+}$:LiNbO$_3$\xspace than for \ercawo. We also report the maximum and minimum values of $T_{z}$ in Table~\ref{tab:comparison_mat} for a concise material comparison.

\begin{widetext}
\begin{center}
\begin{table}[t]
\centering
\begin{tabular}{>{\columncolor{Plum}}c|>{\columncolor{LightCyan}}c|>{\columncolor[gray]{0.9}}c|>{\columncolor{LightCyan}}c|>{\columncolor[gray]{0.9}}c|>{\columncolor{LightCyan}}c|>{\columncolor[gray]{0.9}}c}
 Crystal & $g_{\rm eff}^{\rm max}$  & $g_{\rm eff}^{\rm min}$ & $T_{ff}^{\rm max}$ & $T_{ff}^{\rm min}$ & $\vec{B}$ orientation for  $T_{ff}^{\rm max}$ & $\vec{B}$ orientation for  $T_{ff}^{\rm min}$   \\
\hline
\hline
\ercawo & 8.38	& 1.25 &	199 $\mu$s & 149 $\mu$s &	$\theta = 8 ^\circ$ & $\theta = 0 ^\circ$ \\
\hline
Er$^{3+}$:LiNbO$_3$\xspace & 15.14 &	2.15 &	15.8 ms &	6.55 $\mu$s &	$\theta = 0 ^\circ$ & $\theta = 90 ^\circ$ \\
\hline
\eryso & 11.7  &	1.7	& 165 ms & 	32 $\mu$s &	$\phi = 133 ^\circ$ &	$\phi = 27 ^\circ$ \\
\hline
\ndyso & 2.87 &	0.98 &	299 ms &	5.2 ms & $\phi = 104 ^\circ$ & $\phi = 12 ^\circ$ \\
\end{tabular}
\caption{Comparison between materials. For \ercawo and Er$^{3+}$:LiNbO$_3$\xspace, $\theta$ is the angle of $\vec{B}$ with respect to the $c$-axis (see main text). For \eryso and \ndyso, we restrict our study to magnetic field in the $(D_1,D_2)$ plane for simplicity. The angle $\phi$ is then the angle of $\vec{B}$ in the $(D_1,D_2)$ plane (see \ref{sec:theory_eryso}). The $g$-factor value $g_{\rm eff}^{\rm max}$ (resp. $g_{\rm eff}^{\rm min}$) corresponds to the maximum (resp. minimum) $T_{ff}^{\rm max}$ lifetime (resp. $T_{ff}^{\rm min}$) }
\label{tab:comparison_mat}
\end{table}
\end{center}
\end{widetext}

\subsection{Extrapolation to Nd$^{3+}$\xspace in \yso}

In a comparative approach, we propose to investigate the change of rare-earth ion dopant in the same matrix. We keep \yso because of its remarkable properties, (such as its strong anisotropy) and its wide use in the community. Concerning the dopant, Nd$^{3+}$\xspace is also used in several applications. This latter is also a Kramers ions, but its $g$-tensor exhibits smaller components. Indeed, the maximum component in the principal axis is 4.17 for \ndyso \cite{wolfowicz_coherent_2015}, as compared to 14.64 for \eryso, whereas the two other components have the same order of magnitude. Figure \ref{fig:Tz_YSO_Er_Nd} shows the spin lifetime as a function of the orientation of the magnetic field in the $(D_1,D_2)$ plane for \ndyso and \eryso. This latter serves as a reference with the experimental investigation reported in Fig.\ref{fig:Tz10ppm}. We observe that the maximum is about the same, but that the minimum is larger for \ndyso. Thus, the anisotropy is enhanced for \eryso because of the very large $g$-tensor component.

\begin{figure}[t!]
\centering
\includegraphics[width=0.95\columnwidth]{Tz_YSO_Er_Nd.pdf}
\caption{Dependence of the spin lifetime as a function of the orientation of the external magnetic field in the $(D_1,D_2)$ plane for \ndyso (red line) and \eryso (blue line) crystals. We choose a concentration of 10 ppm and a inhomogeneous spin broadening $\Gamma_{\rm inh}^{\rm spin}$ of 5 MHz in both cases.}
\label{fig:Tz_YSO_Er_Nd}
\end{figure}

\subsection{Expected concentration dependence}

We finally present the expected concentration dependence of the spin lifetime for all the materials discussed in this work. Figure \ref{fig:Tz_vs_C} shows the lifetimes for the optimized magnetic field orientation (leading to the longest time $T_{ff}^{\rm max}$) in \eryso, \ndyso, Er$^{3+}$:LiNbO$_3$\xspace and \ercawo as a function of the dopant concentration. We retrieve the orders of magnitude difference now expressed in terms of concentration. For example, if we aim at a FF-limited spin lifetime of 10~s (with the best orientation) to favor the hole-burning dynamics, this can be only achieved at concentration lower than 1.8 ppm, 1.3 ppm, 0.4 ppm and 0.05 ppm for \ndyso, \eryso, Er$^{3+}$:LiNbO$_3$\xspace and \ercawo respectively. The lowest concentration level can totally prevent the optical measurements due to the lack of absorption. This aspect should be kept in mind when different combination of dopant and matrices are considered for SHB application at low magnetic field.

\begin{figure}[t!]
\centering
\includegraphics[width=0.95\columnwidth]{TzvsCwithNd.pdf}
\caption{Dependence of the spin lifetime on the doping concentration for \eryso, \ndyso, Er$^{3+}$:LiNbO$_3$\xspace and \ercawo. We choose an inhomogeneous spin broadening $\Gamma_{\rm inh}^{\rm spin}$ of $5$~MHz and the best magnetic field orientation, i.e. the one which gives the longest $T_{ff}$, in all cases.}
\label{fig:Tz_vs_C}
\end{figure}

\section{Conclusion}

We have shown experimentally that the erbium flip-flop dynamics is extremely anisotropic in \yso for different orientations of the magnetic field. We have verified that an order of magnitude anisotropy in the $g$-factor translates into four orders of magnitude for the flip-flop rate \cite{Cruzeiro}. In this regime, typically at low magnetic field where many experiments are operated, this scaling severly limits the spin lifetime and may directly impact the optical pumping of ensembles. Depending on the sample under consideration, a proper choice of magnetic field orientation and dopant concentration has to be made.

It should be noted that the flip-flop reduction may impose a severe additional constraint on EPR experiments when for example a large microwave coupling is targeted. For example in \eryso, a bias magnetic field minimizing the flip-flops (large $g$-factor) irremediably produces a weak coupling to the oscillating excitation field in the orthogonal direction. In other words, a proper trade-off has to be considered between cross-relaxation and an strong microwave coupling.

In any case, our study can be used as a guideline to understand the spin dymanics of other Kramers ions in different matrices, as soon as the $g$-tensor is well caracterized.

\section*{Acknowledgements}
We have received funding from the Investissements d'Avenir du LabEx PALM ExciMol and OptoRF-Er (ANR-10-LABX-0039-PALM) and from the ITMO Cancer: AVIESAN (National Alliance for Life Sciences \& Health) within the framework of the Cancer Plan.


\section{Introduction}
\subsection{Real-time systems and automata based modeling}
Real-time systems maintain timely interaction with the environment. The timing constraints are crucial to such systems and it may lead to dramatic consequences when these constraints are not met. We therefore would like to model such systems using automata theory and model check the timed properties of real-time systems. However, we don't have the necessary resources to model time in finite automata. In 1994, Alur and Dill published the results of their study as the theory of timed automata [1], thereby solving this problem. Timed automata theory is an extension to finite automata theory that describes a way of measuring real-time event occurrences in such systems. It equips us with the required resources in verifying these real-time systems and obtaining provably correct guarantees of timed properties in them. We will study various properties of timed automata and prove some of them in section II. We will also look at a class of systems known as Hybrid Systems that have mixed discrete-continuous behavior. We will look at the hybrid automata theory for modeling this behavior in section III. In section IV we will look at some of the applications of systems known as Cyber-Physical systems that have these hybrid properties. In section V we look at a case-study analysis of the peterson's mutual exclusion protocol in a tool for modeling and verifying real-time systems, and discuss open problems and conclusions in section VI.

\subsubsection{Real-time systems} Real-time systems are encountered in many instances, and we interact with them more often than we realize. Some examples of real-time systems are event response systems like airbag systems, closed-loop control like cruise control system in a car, aircraft control systems, cardiac pacemakers. Before looking at automata for modeling real-time systems, let us first take a look at a class of words that have an embedded time component within them.

\subsubsection{Timed words} An alphabet $\Sigma$ is defined over a finite set of letters. A timed word is a tuple (w, t) where w is a word over the alphabet $\Sigma$, w=$a_1$$a_2$$a_3$...$a_n$ and t is a time sequence, t=$t_1$$t_2$$t_3$...$t_n$ where $t_1{\leq}t_2{\leq}t_3{\leq}$...$t_n$ and $t_i$ $\in$ $R_{\geq0}$, the set of non-negative real numbers. A timed language is a set of words in L such that $L$ $\in$ $T{\Sigma}^*$ is a property over timed words [2].

\subsubsection{Automaton modeling timed properties} Fig. 1. shows an automaton modeling a lamp. It also describes the timed properties of the lamp. The lamp starts in off state. When the switch is pressed in off state, it takes the transition to low state while resetting a clock, y to zero. From low state if the switch is pressed again before 5 seconds have elapsed, it transitions to bright state. From low state if the switch is pressed after 5 seconds have elapsed, it takes a transition back to the off state and accepts. From bright state, if the switch is pressed, it transitions to off state and accepts.

\begin{figure}
\begin{center}
\begin{tikzpicture}[scale=0.2]
\tikzstyle{every node}+=[inner sep=0pt]
\draw [black] (22.1,-26.6) circle (3);
\draw (22.1,-26.6) node {$off$};
\draw [black] (22.1,-26.6) circle (2.4);
\draw [black] (39.9,-26.6) circle (3);
\draw (39.9,-26.6) node {$low$};
\draw [black] (56.8,-26.6) circle (3);
\draw (56.8,-26.6) node {$bright$};
\draw [black] (22.463,-23.636) arc (-196.62661:-343.37339:8.91);
\fill [black] (39.54,-23.64) -- (39.79,-22.73) -- (38.83,-23.01);
\draw (31,-16.78) node [above] {$press?\mbox{ }y:=0$};
\draw [black] (38.36,-29.161) arc (-39.96724:-140.03276:9.604);
\fill [black] (23.64,-29.16) -- (23.77,-30.09) -- (24.54,-29.45);
\draw (31,-33.1) node [below] {$press?\mbox{ }y>=5$};
\draw [black] (56.344,-29.561) arc (-13.69418:-166.30582:17.388);
\fill [black] (22.56,-29.56) -- (22.26,-30.46) -- (23.23,-30.22);
\draw (39.45,-43.33) node [below] {$press?$};
\draw [black] (42.9,-26.6) -- (53.8,-26.6);
\fill [black] (53.8,-26.6) -- (53,-26.1) -- (53,-27.1);
\draw (48.35,-27.1) node [below] {$press?\mbox{ }y<5$};
\end{tikzpicture}
\caption{Automaton modeling timed properties of a lamp} \label{fig1}
\end{center}
\end{figure}

\subsection{Timed Automata} A timed automaton is composed of a finite automaton and a finite set of real-valued clocks. All the clocks values increase at the same rate. Guards can be placed on the transitions of the automaton using which we can enable or disable that transition, thereby constraining the behavior of the automaton and describing timed properties through the language of the automaton. Guards are comparisons of clock values with constants that are non-negative rational numbers and can either evaluate to true or false, g: $x{\leq}c$ $\mid$ $x{\geq}c$ $\mid$ ${\neg}g$ $\mid$ $g{\land}g$ where $x$ $\in$ $Clocks$, $c$ $\in$ $Q_{\geq0}$. The clocks can be reset.

\subsubsection{Formal defn.} A timed automaton M is defined as a tuple, M = (Q, $\Sigma$, C, Inv, $\delta$, $q_0$, F). Q is the finite set of states. $\Sigma$ is a finite set of actions. C is a finite set of clocks. Inv associates each location with an invariant. $\delta$ is a set of transitions. (q,a,g,r,q') is a transition from \textit{q} to \textit{q'}, executing an action \textit{a}, satisfying a guard \textit{g} and clock resets in \textit{r}. There are two types of transitions namely a location switch and a time switch. By elapsing time, and satisfying the location invariant, the automaton can stay in the same state thus making a time switch. By satisfying the guard and taking a transition to another location, the automaton makes a location switch.

\section{Timed Automata properties}
\subsection{Non-deterministic timed automata}
There can be non-determinism concerning location and non-determinism concerning time. $\epsilon$ transitions may also introduce non-determinism w.r.t location. Fig. 2. shows an example of location non-determinism. Upon reading a \textit{b} from the state $s_1$ we cannot determine deterministically if the automaton is going to stay in the state $s_1$ or if it is going to transition to the state $s_2$. Fig. 3. shows an example of time non-determinism. Upon reading an \textit{a} from the state $s_1$ we cannot determine deterministically the time of occurrence of the event (reading an a) as long as it satisfies the guard condition t $\leq$ 4 and location invariant t $\leq$ 5.

\begin{figure}
\begin{center}
\begin{tikzpicture}[scale=0.2]
\tikzstyle{every node}+=[inner sep=0pt]
\draw [black] (28.5,-25.4) circle (3);
\draw (28.5,-25.4) node {$s_1$};
\draw [black] (54.5,-25.4) circle (3);
\draw (54.5,-25.4) node {$s_2$};
\draw [black] (54.5,-25.4) circle (2.4);
\draw [black] (31.5,-25.4) -- (51.5,-25.4);
\fill [black] (51.5,-25.4) -- (50.7,-24.9) -- (50.7,-25.9);
\draw (41.5,-25.9) node [below] {$b,\mbox{ }x:=0$};
\draw [black] (27.177,-22.72) arc (234:-54:2.25);
\draw (28.5,-18.15) node [above] {$b$};
\fill [black] (29.82,-22.72) -- (30.7,-22.37) -- (29.89,-21.78);
\draw [black] (53.177,-22.72) arc (234:-54:2.25);
\draw (54.5,-18.15) node [above] {$a,b,\mbox{ }x\neq1$};
\fill [black] (55.82,-22.72) -- (56.7,-22.37) -- (55.89,-21.78);
\end{tikzpicture}
\caption{Non-deterministic Timed Automaton (location non-determinism)} \label{fig2}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\begin{tikzpicture}[scale=0.2]
\tikzstyle{every node}+=[inner sep=0pt]
\draw [black] (28.5,-25.4) circle (3);
\draw (28.5,-25.4) node {$y\leq5$};
\draw [black] (54.5,-25.4) circle (3);
\draw (54.5,-25.4) node {$s_2$};
\draw [black] (54.5,-25.4) circle (2.4);
\draw [black] (31.5,-25.4) -- (51.5,-25.4);
\fill [black] (51.5,-25.4) -- (50.7,-24.9) -- (50.7,-25.9);
\draw (41.5,-25.9) node [below] {$b,\mbox{ }y:=0$};
\draw [black] (27.177,-22.72) arc (234:-54:2.25);
\draw (28.5,-18.15) node [above] {$a,\mbox{ }y\leq4$};
\fill [black] (29.82,-22.72) -- (30.7,-22.37) -- (29.89,-21.78);
\draw [black] (53.177,-22.72) arc (234:-54:2.25);
\draw (54.5,-18.15) node [above] {$a,b\mbox{ }y\leq10$};
\fill [black] (55.82,-22.72) -- (56.7,-22.37) -- (55.89,-21.78);
\end{tikzpicture}
\caption{Non-deterministic Timed Automaton (time non-determinism)} \label{fig3}
\end{center}
\end{figure}

\subsection{Deterministic timed automata}
A timed automaton is deterministic if a) It has only one initial location, b) It doesn't have $\epsilon$ transitions, c) Event determinism: Two edges with same source and same label have disjoint guards ($g_1$ $\cap$ $g_2$ = $\emptyset$), d) Time determinism: For every transition, the intersection of g with $I_q$ is at most a singleton.

\subsection{Expressiveness of $\epsilon$ transitions}
$\epsilon$ transitions add to the expressiveness of a timed automaton, i.e. the language recognized by a timed automaton with $\epsilon$ transitions may not be recognized by a timed automaton without the $\epsilon$ transitions. Take the automaton in Fig. 4. as an example. This automaton accepts timed words over \textit{a} such that every occurrence time is an integer and no two \textit{a} events occur at the same time. This language cannot be accepted by a timed automaton if $\epsilon$ switches are not allowed. If the largest constant in such timed automaton is \textit{c}, then a timed automaton without $\epsilon$ transitions cannot distinguish between the words (a, c+1) and (a, c+1.1). In this case, it can't distinguish between (a, 2) and (a, 2.1).

\begin{figure}
\begin{center}
\begin{tikzpicture}[scale=0.2]
\tikzstyle{every node}+=[inner sep=0pt]
\draw [black] (37.6,-32.3) circle (3);
\draw [black] (37.6,-32.3) circle (2.4);
\draw [black] (36.277,-29.62) arc (234:-54:2.25);
\draw (37.6,-25.05) node [above] {$a,\mbox{ }x=1,\mbox{ }x:=0$};
\fill [black] (38.92,-29.62) -- (39.8,-29.27) -- (38.99,-28.68);
\draw [black] (38.923,-34.98) arc (54:-234:2.25);
\draw (37.6,-39.55) node [below] {$\epsilon,\mbox{ }x=1,\mbox{ }x:=0$};
\fill [black] (36.28,-34.98) -- (35.4,-35.33) -- (36.21,-35.92);
\end{tikzpicture}
\caption{Expressiveness of $\epsilon$ transitions} \label{fig4}
\end{center}
\end{figure}

\subsection{Timed regular languages}
A timed language is called timed regular if it can be accepted by a timed automaton.

\subsection{Closure properties}
\begin{theorem}Timed regular languages are closed under the following operations.\end{theorem} Let us consider two timed regular languages, $L_1$, $L_2$.
\begin{enumerate}
   \item \textit{Union}: Union of two timed regular languages, ${L_1}\cup{L_2}$ is timed regular.
   \item \textit{Intersection}: Intersection of them, ${L_1}\cap{L_2}$ is timed regular.
   \item \textit{Projection}: Projection of a timed regular language is timed regular.
   \item \textit{Untime}: If L is timed regular, then untime(L) is $\omega$-regular [3].
\end{enumerate}

Closure under Union and Intersection are established by constructing product of the timed automata. Projection can be proved by labeling transitions with $\epsilon$. The proof for untime(L) being $\omega$-regular is established by region construction as shown in [1].

\subsection{Closure under complementation}
\begin{theorem}Timed regular languages are not closed under complementation [3].\end{theorem}
\begin{proof}
Let $\Sigma$ = \{a, b\}. L is a timed language consisting of timed words \textit{w}, containing an $\textit{a}$ event at some time $\textit{t}$ such that no event occurs at time $\textit{t + 1}$. \textit{L} is accepted by the timed automaton in Fig. 5. We will show that $\overline{L}$ is not timed regular. \textit{untime(L)} accepts \textit{(a+b)*a(a+b)*}.

\begin{figure}
\begin{center}
\begin{tikzpicture}[scale=0.2]
\tikzstyle{every node}+=[inner sep=0pt]
\draw [black] (28.5,-25.4) circle (3);
\draw (28.5,-25.4) node {$s_1$};
\draw [black] (54.5,-25.4) circle (3);
\draw (54.5,-25.4) node {$s_2$};
\draw [black] (54.5,-25.4) circle (2.4);
\draw [black] (31.5,-25.4) -- (51.5,-25.4);
\fill [black] (51.5,-25.4) -- (50.7,-24.9) -- (50.7,-25.9);
\draw (41.5,-25.9) node [below] {$a,\mbox{ }x:=0$};
\draw [black] (27.177,-22.72) arc (234:-54:2.25);
\draw (28.5,-18.15) node [above] {$a,b$};
\fill [black] (29.82,-22.72) -- (30.7,-22.37) -- (29.89,-21.78);
\draw [black] (53.177,-22.72) arc (234:-54:2.25);
\draw (54.5,-18.15) node [above] {$a,b\mbox{ }x\neq1$};
\fill [black] (55.82,-22.72) -- (56.7,-22.37) -- (55.89,-21.78);
\end{tikzpicture}
\caption{Timed Automaton for disproving closure on complementation} \label{fig6}
\end{center}
\end{figure}

Consider timed regular language L', consisting of timed words \textit{w}, such that untimed word of \textit{w} is in \textit{a*b*}. All \textit{a} events happen before time 1 and no two \textit{a} events happen at same time. Observe the word $a^{m}b^{n}$ $\in$ untime($\overline{L}{\cap}L'$) iff m $\geq$ n. Timed regular languages are closed under intersection, \textit{untime(L)} of the timed regular language \textit{L} is regular, but since we also know that $a^{m}b^{n}$ m $\geq$ n is not regular, we proved that $\overline{L}$ is not timed regular. \end{proof}

\subsection{Emptiness}
Given a timed automaton over finite words, is the language accepted by it empty? This problem is known as the \textbf{emptiness} problem for timed automata. It can be analyzed by checking if there is a run of the automaton from an initial state to a final state.
\begin{proof}
We have a problem in performing this check since the number of configurations is uncountably infinite due to real-valued clock times. The solution to this problem was proposed in [1] by constructing a finite region graph using region equivalence technique. The caveat here is that although there are infinite configurations, some clock valuations are equivalent and hence a timed automaton cannot distinguish between them. By grouping such configurations using region equivalence technique, we construct a finite region graph. The problem thus reduces to a reachability problem over this finite region graph, which is decidable. Hence it is proved that emptiness problem for timed automata is decidable.
\end{proof}

\subsection{Universality}
Given a timed automaton over finite words, is the language accepted by it the set consisting of every timed word? This problem is known as the \textbf{universality} problem for timed automata. This problem can be proved undecidable by reducing the halting problem over two-counter machines, which is known to be undecidable to the universality problem [3].

\subsection{Language inclusion}
Given two timed automata \textit{A} and \textit{B} over finite words, checking if L(A) $\subseteq$ L(B) is known as the \textbf{language inclusion} problem. This problem can be proved undecidable by reducing the halting problem over two-counter machines, which is known to be undecidable to the language inclusion problem [3].

\subsection{Decidability in special cases}
Let us consider the universality problem on a given timed automaton having at most one clock. In this special case, this problem has been proved to be decidable by Abdulla et al. in [4]. Given two timed automata A and B, such that the timed automaton B only has one clock, and the only constant appearing in the clock constraints of B is 0, the language inclusion problem on A and B, formulated as checking if L(A) $\subseteq$ L(B) has been proved to be decidable by Ouaknine et. al. in [5].

\subsection{Decidability of determinizability}
In [6], E. Asarin has posed an open question concerning timed automata as follows. \textit{Given a timed automaton A, is it possible to decide whether it is equivalent to a deterministic one?} This question remained open until O. Finkel has proved in [7] that this problem is undecidable.

\section{Hybrid Systems}
\subsection{Background}
Hybrid Systems are dynamical systems with interacting continuous-time dynamics and discrete-time dynamics. The evolution of the state of a continuous time system is described by an \textbf{ordinary differential equation} (ODE), $\dot{x} = Ax$, whereas the evolution of the state of a discrete-time system is described by a \textbf{difference equation}, $x_{k+1}$ = A$x_k$. The continuous time dynamics are referred to as flows and the discrete-time dynamics are referred to as jumps. The flows cause the system's state to make a smooth continuous transition and jumps cause the system to transition to a different set of flow equations by making a discrete jump to the new continuous-time dynamics. The hybrid behavior arrives in various contexts, a) Continuous systems with a phased operation, like the bouncing ball, biological cell growth to name a few, b) Continuous systems controlled by discrete logic, like control modes for complex systems, c) Coordinating processes, like multi-agent systems. When formally modeling the hybrid system behavior, the modeling paradigm we choose needs to have both continuous and discrete parts to it [8].

\subsection{Bouncing ball}
Let us consider an example of a ball of mass \textit{m} freely falling from height \textit{x} $\geq$ 0 that bounces after hitting the ground. The analysis of this system has two parts to it.
\subsubsection{Part-I Free Fall}
In free fall while x $\geq$ 0, the ball is under the influence of gravity. The equation of motion satisfies $\ddot{x}$ = -g, where \textit{x} is the position of the ball and \textit{g} is the gravitational constant. We reduce the order of the ODE by substituting $\dot{x}$ = v, where \textit{v} is the velocity of the ball. We have $\dot{v}$ = -g. This describes the continuous time dynamics of the system.
\subsubsection{Part-II Bouncing}
When the ball is at x = 0, and has a velocity downwards ($v<0$), the ball bounces on the ground. There is a loss of velocity due to deformation and friction. The velocity discretely jumps to a different value as per the equation v := -cv, where $c<1$. The velocity's magnitude is reduced and it's direction flips.
\subsubsection{Hybrid Automaton}
The hybrid automaton for bouncing ball is illustrated in Fig. 6. $x_1$ represents the vertical position and $x_2$ represents the velocity of the ball [9]. In Fig. 7 we illustrate the bouncing ball simulation in MATLAB, where the yellow legend corresponds to position of the ball, and blue legend corresponds to velocity of the ball.

\begin{figure}
\begin{minipage}{.5\textwidth}
\begin{center}
\begin{tikzpicture}[scale=0.175]
\tikzstyle{every node}+=[inner sep=0pt]
\draw [black] (38.2,-33.4) circle (5.9);
\draw (38.2,-31.4) node {$\dot{x_1}\mbox{ }=\mbox{ }x_2$};
\draw (38.2,-34.4) node {$\dot{x_2}\mbox{ }=\mbox{ }-g$};
\draw (38.2,-36.4) node {$m$};
\draw [black] (32.877,-30.72) arc (259:-12:4.25);
\draw (32.89,-17.96) node [above] {$x_1=0 \wedge x_2\leq0$};
\draw (32.89,-19.96) node [above] {impact};
\draw (32.89,-21.96) node [above] {${x_1}' = {x_2} \wedge {x_2}' = -c{x_2}$};
\draw (35.89,-40.96) node [below] {${x_1} \geq 0$};
\fill [black] (37.88,-27.66) -- (38.38,-26.45) -- (37.38,-26.46);
\end{tikzpicture}
\caption{Hybrid automaton} \label{fig6}
\end{center}
\end{minipage}
\begin{minipage}{.5\textwidth}
\begin{center}
\includegraphics[width=6.3cm,height=4.9cm,keepaspectratio]{bouncing_ball_simulations.PNG}
\caption{Bouncing ball simulation} \label{fig7}
\end{center}
\end{minipage}
\end{figure}

\subsubsection{Formal defn.} A hybrid automaton H is defined as a tuple, H = (M, $\Sigma$, $M_0$, X, $\delta$, I, F, J). M is the finite set of control modes, also known as locations. $\Sigma$ is a finite set of actions. Actions are given as sets of differential equations. $M_0$ is a finite initial set of control modes, ${M_0}$ ${\subseteq}$ $M$. X is a finite set of real-valued variables. $\delta$ is a set of transitions. I is the mode-invariant function. F is the mode dependent flow function, characterizing flow at each mode. J is the jump function. The configuration of a hybrid automaton is $(m,v)$ where $m$ ${\in}$ $M$  is a mode/location and v is a variable valuation.

\subsubsection{Reachability}
Linear Hybrid Automata are a class of hybrid automata where its activities, invariants, and transition relations are defined by linear expressions. Given a linear hybrid automaton, is there a run reaching a particular state from the initial state? This problem is known as the \textbf{reachability} problem for linear hybrid automata. In [10], Alur et al. have reduced the halting problem of two-counter machines, which is known to be undecidable to the reachability problem for linear hybrid automata. Hence, it is proved to be undecidable.\\
Further, in [11], S.N. Krishna et al. have reduced the halting problem for two-counter machines to the reachability problem on recursive hybrid automata, thus showing that it is undecidable. However, bounded reachability is still decidable. There exist some incomplete algorithms for reachability. In [12], E. Abraham presents one such incomplete algorithm for linear hybrid automata based on fixed-point computation. Termination of this algorithm corresponds to finding the least fixed-point for the one-step forward reachability starting from the initial set.

\section{Cyber-Physical Systems}
\subsection{Background}
As defined in [13], a cyber-physical system consists of a collection of computing devices communicating with one another and interacting with the physical world through sensors and actuators in a feedback loop. There are quite a few applications of such systems everywhere, whether it is smart buildings or medical devices or even automobiles. One of the distinguishing characteristics of cyber-physical systems is that they are reactive in nature, meaning these systems interact with the environment in an ongoing manner through inputs and outputs, like a program for a cruise controller in a car just to name one. They are concurrent in nature, meaning as opposed to sequential computation, there are multiple threads known as processes executing concurrently. They have mixed discrete-continuous behavior such as the hybrid systems. They often have real-time system characteristics and are in the most safety-critical areas where errors could lead to catastrophic eventualities. We will study some crucial applications of cyber-physical systems and look at automata based modeling and analysis of such systems.

\subsection{Cardiac Pacemaker}
Implantable cardiac pacemakers are life-saving devices that help patients with heart arrhythmia conditions. Ironically, there are bugs even in such safety-critical devices. To guarantee the correct operation of such devices, we need a model of the heart that captures the physiological conditions of the heart and respond to pacemaker outputs, using which we can verify the safety properties of the pacemaker. The pacemaker functions autonomously according to the device algorithm that is implemented onto it, interacts with its environment through sensors, actuators and other devices. It is, therefore a perfect example of a cyber-physical system as defined in the previous section. Timed automata is an appropriate formalism for such a heart model since most timing behaviors of heart can be captured by timed automata. In [14] Z. Jiang et. al. have proposed a real-time heart model based on timed automata and capture such crucial timing properties of the heart to then verify the cardiac pacemaker.

\section{Uppaal case study}
\subsection{Background}
Uppaal is a model checker for real-time systems. We can conveniently model a real-time system as a network of timed automata running concurrently. Uppaal is used to verify system properties specified in CTL over the timed automata system model [15]. The CTL formulae can be specified over paths, and are classified broadly as \textit{safety}, \textit{liveness} and \textit{reachability} properties. Uppaal runs on a client-server architecture and is split into the backend model checking engine and frontend GUI communicating via TCP/IP protocol. It has a simulator where the user can run the system manually, and a verifier where the user may specify the properties and run the model checker on the system.

\subsection{Peterson's mutual exclusion}
The idea of mutual exclusion is that two processes which have critical sections cannot enter those sections at the same time. It is a technique for avoiding race conditions and keeping the result deterministic while running the processes concurrently. The Peterson's mutual exclusion algorithm is illustrated below.
\begin{algorithm}
\caption{Peterson's Mutex}\label{mutex}
\begin{minipage}{.5\textwidth}
\begin{algorithmic}[1]
\Procedure{Process 1}{}\\
req1 = 1;\\
turn = 2;\\
while(turn!=1 \&\& req2!=0);\\
\textit{//critical section}\\
job1();\\
req1 = 0;
\EndProcedure
\end{algorithmic}
\end{minipage}
\begin{minipage}{.5\textwidth}
\begin{algorithmic}[1]
\Procedure{Process 2}{}\\
req2 = 1;\\
turn = 1;\\
while(turn!=2 \&\& req1!=0);\\
\textit{//critical section}\\
job2();\\
req2 = 0;
\EndProcedure
\end{algorithmic}
\end{minipage}
\end{algorithm}

Fig. 8. illustrates the results of Uppaal model checker on the correctly implemented mutex algorithm as above. First property species that there exists a run where critical section is reachable as $E<> (P1.CS)$. Second property captures the mutual exclusion property, where on all paths both P1 and P2 cannot be in the critical section at the same time as A$[]$ not (P1.CS and P2.CS). We see that both the properties are satisfied in this case.

\begin{figure}
\begin{center}
\includegraphics[scale=0.35]{uppaal_2.png}
\caption{Model checking mutex} \label{fig8}
\end{center}
\end{figure}

Fig. 9. illustrates incorrectly implemented mutex, where two processes are seen to be in the \textit{criticial section} at the same time in the simulator. Fig. 10 illustrates the results where Uppaal says that mutual exclusion property is not satisfied on the incorrectly implemented mutex. Mutual exclusion property has failed to satisfy the model.

\begin{figure}
\begin{minipage}{.5\textwidth}
\begin{center}
\includegraphics[scale=0.25]{uppaal_3.png}
\caption{Mutex protocol with bug} \label{fig9}
\end{center}
\end{minipage}
\begin{minipage}{.5\textwidth}
\begin{center}
\includegraphics[scale=0.31]{uppaal_4.png}
\caption{Model checking mutex with bug} \label{fig10}
\end{center}
\end{minipage}
\end{figure}

\section{Conclusion}
The motive of this work is to study extensions to finite automata theory used to model systems of the 21st century that are becoming increasingly complex to verify, due to the confluence of sensors, actuators, real-time constraints, and their concurrent and reactive nature. The automata theory based modeling paradigms discussed above are helping us design computer-based models of such systems that can be used to analyze their behavior mathematically and ultimately build safety-critical systems that are more dependable and secure although complex.\\

As we use model checking to verify the properties of such systems, we also see that some seemingly simple problems like reachability don't have a decidable algorithm that works for any general system. Researchers in this area are working towards devising algorithms to model check these systems. These algorithms although while being incomplete algorithms, terminate in most cases that are very relevant to the system analysis. Future work also lies ahead in researching methods to synthesize correct-by-construct designs for controllers and other components involved in such real-time hybrid systems. Some open problems in hybrid systems are concerning non-linear hybrid systems. Scalability is a challenge for most of the modeling tools available to us currently. Furthermore, high-dimensional systems pose a unique challenge to reachability analysis in these tools. Recent works such as [16], show promise in this direction. Schupp et al. [17] present an analysis of the current challenges in verification of hybrid systems.