# Unsafe Probabilities and Risk Contours for Stochastic Processes using Convex
Optimization

Jared Miller1, Matteo Tacchi2, Didier Henrion3, Mario Sznaier4 1 J. Miller is
with the Automatic Control Laboratory (IfA), Department of Information
Technology and Electrical Engineering (D-ITET), ETH Zürich, Physikstrasse 3,
8092, Zürich, Switzerland (e-mail: jarmiller@control.ee.ethz.ch).2 M. Tacchi
is with Univ. Grenoble Alpes, CNRS, Grenoble INP (Institute of Engineering
Univ. Grenoble Alpes), GIPSA-lab, 38000 Grenoble, France. (e-mail:
matteo.tacchi@gipsa-lab.fr)3 D. Henrion is with the Polynomial Optimization
group of LAAS-CNRS, Université de Toulouse, CNRS, Toulouse, France; and the
Faculty of Electrical Engineering, Czech Technical University in Prague, Czech
Republic. (henrion@laas.fr)4 M. Sznaier is with the Robust Systems Lab, ECE
Department, Northeastern University, Boston, MA 02115. (e-mail:
msznaier@coe.neu.edu).J. Miller and M. Sznaier were partially supported by NSF
grants CNS–1646121, ECCS–1808381 and CNS–2038493, AFOSR grant
FA9550-19-1-0005, and ONR grant N00014-21-1-2431. J. Miller was partially
supported by the Swiss National Science Foundation under NCCR Automation,
grant agreement 51NF40_180545.

###### Abstract

This paper proposes an algorithm to calculate the maximal probability of
unsafety with respect to trajectories of a stochastic process and a hazard
set. The unsafe probability estimation problem is cast as a primal-dual pair
of infinite-dimensional linear programs in occupation measures and continuous
functions. This convex relaxation is nonconservative (to the true probability
of unsafety) under compactness and regularity conditions in dynamics. The
continuous-function linear program is linked to existing probability-
certifying barrier certificates of safety. Risk contours for initial
conditions of the stochastic process may be generated by suitably modifying
the objective of the continuous-function program, forming an interpretable and
visual representation of stochastic safety for test initial conditions. All
infinite-dimensional linear programs are truncated to finite dimension by the
Moment-Sum-of-Squares hierarchy of semidefinite programs. Unsafe-probability
estimation and risk contours are generated for example stochastic processes.

##  1 Introduction

This work performs risk analysis of stochastic processes by computing the
probability that trajectories will reach an unsafe set $X_{u}$. Trajectories
evolve over a maximal time horizon of $T$ in a state space
$X\subseteq\mathbb{R}^{n}$, beginning in an initial set $X_{0}\subseteq X$ and
possibly entering the unsafe set $X_{u}\subseteq X$. We consider stochastic
processes described by a generator $\mathcal{L}$ (e.g., SDE Stochastic
Differential Equation (SDE)), for which $X_{t}$ is the time-dependent state
probability distribution of the process starting from $X_{0}$ and terminating
upon exit from $X$ (with exit time distribution
$\tau_{X}=\inf\\{t:X_{t}\in\partial X\\}$). The worst-case probability of
unsafety for $\mathcal{L}$-trajectories starting at a point $x_{0}\in X$
within a time horizon of $[t_{0},T]$ is:

$\displaystyle P^{*}(t_{0},x_{0})=$
$\displaystyle\sup_{t^{*}\in[t_{0},T]}\quad\textrm{Prob}_{X_{t^{*}}}[x\in
X_{u}]$ (1a) $\displaystyle x(t)\ \text{follows}\ \mathcal{L}\quad\forall
t\in[t_{0},\min(t^{*},\tau_{X})]$ (1b) $\displaystyle x(0)=x_{0}.$ (1c)

The stopping time $\min(t^{*},\tau_{X})$ in (1b) also ensures that
trajectories $x(t)$ remain in $X$ for all relevant $t$. The worst-case
probability of unsafety for $\mathcal{L}$-trajectories over an initial set
$X_{0}\subseteq X$ is:

$\displaystyle P^{*}(t_{0},X_{0})=$ $\displaystyle\sup_{x_{0}\in
X_{0}}P^{*}(t_{0},x_{0})$ (2a) $\displaystyle=$
$\displaystyle\sup_{t^{*}\in[0,T],\ x_{0}\in
X_{0}}\quad\textrm{Prob}_{X_{t^{*}}}[x\in X_{u}]$ (2b) $\displaystyle x(t)\
\text{follows}\ \mathcal{L}\quad\forall t\in[t_{0},\min(t^{*},\tau_{X})]$ (2c)
$\displaystyle x(0)=x_{0}.$ (2d)

The initial set of $X_{0}$ is safe if $P^{*}(t_{0},X_{0})=0$ and is unsafe if
$P^{*}(t_{0},X_{0})=1$. Any other value $P^{*}(t_{0},X_{0})\in(0,1)$ returns a
maximal probability of unsafety of the SDE.

The map $(t_{0},x_{0})\mapsto P^{*}(t_{0},x_{0})$ can be interpreted as a risk
function in trajectory planning, in which level-sets of $P^{*}$ are risk
contours with constant unsafe probability. This paper will also develop a
piecewise-polynomial sequence of risk functions that will converge from above
in an $L_{1}$ sense to the optimal map $P^{*}$. Superlevel sets of any one of
these risk-function-proxies can be used to identify unsafe regions.

Safety verification of stochastic systems includes the tasks of stochastic
reachability [1], reach-avoid analysis [2, 3], ruin [4], and performance
checking [5].

Infinite-dimensional convex Linear Program (LP) formulations of generically
nonconvex and nonlinear stochastic trajectory problems can be used to analyze
and control system executions, with no conservatism added under appropriate
compactness and regularity conditions. The infinite-dimensional LP must be
truncated into finite-dimensional programs in order to obtain numerical
solutions. Algorithms to perform this truncation for generic stochastic
processes include gridding [6, 7], radial basis function selection [8], and
random sampling [9].

Barrier functions [10] provide a level-set certificates of invariance with
probability $(1-\epsilon)$. Control barrier functions [11, 12, 13] may also be
generated to form controllers with stochastic guarantees of safety.

In the specific case where the stochastic process and all problem data have a
rational representation, the moment- Sum of Squares (SOS) hierarchy [14] can
be employed to truncate the LP into a hierarchy of Semidefinite Programs in
increasing size, parameterized by the polynomial degree. Application of the
moment-SOS hierarchy for stochastic analysis and control includes regional
verification [15], reach-avoid estimation [16], reach-avoid control [17],
region of attraction estimation [18], chance-constrained trajectory planning
via moment propagation [19], exit-time estimation [20], and (conditional)
value-at-risk upper-bounding [5].

This paper continues a research direction involving safety quantification of
nonlinear systems. The safety of a set $X_{0}$ in (2) is quantified by its
worst-case probability of unsafety. Other such quantities include the expected
time a trajectory spends in the unsafe set [21], constraint violations
representing an expanded unsafe set [22], the distance of closest approach
[23], value-at-risk of state functions [5], and the perturbation intensity
necessary to crash [24].

The two most relevant prior works to our paper are [25] and [26]. The work in
[25] presents infinite-dimensional LPs to generate superlevel-based outer
approximations of probability-$p$-safe sets of an SDE for a fixed initial
distribution. Over the course of their derivations, they present a
nonconservative formulation in Lemma 2 that is similar to our approach. Our
problem evaluates safety of a specific initial set $X_{0}$ with respect to a
free initial distribution, performs risk contour analysis, and extends the
framework towards more general stochastic processes. The work in [26]
synthesizes risk contours for optimization problems with distributional
uncertainty (with an emphasis on obstacles and path planning), and our paper
extends this work by analyzing safety of stochastic dynamical systems.

This paper has the following structure: Section 2 introduces preliminaries
such as notation, stochastic processes, and SOS proofs of nonnegativity.
Section 3 poses infinite-dimensional convex LPs to compute the worst-case
probability of unsafety of a stochastic process w.r.t. an initial set $X_{0}$.
Section 4 modifies the worst-case unsafety LPs to create risk contours that
upper bound the true risk map $P^{*}(t_{0},x_{0})$. Section 5 truncates the
unsafety LPs into a converging sequence of finite-dimensional SDPs using the
moment-SOS hierarchy. Section 6 demonstrates the computation of unsafe
probabilities and risk contours on example stochastic process with polynomial
dynamics. Section 7 concludes the paper.

##  2 Preliminaries

BSA

    Basic Semialgebraic
LP

    Linear Program
ODE

    Ordinary Differential Equation
PSD

    Positive Semidefinite
SDE

    Stochastic Differential Equation
SDP

    Semidefinite Program
SOS

    Sum of Squares
WSOS

    Weighted Sum of Squares

###  2.1 Notation

The subset of natural numbers between $a$ and $b$ is $a..b\subset\mathbb{N}$.
The minimum of two quantities will be denoted as $a\wedge b=\min(a,b)$. The
set of polynomials with real-valued coefficients in an indeterminate
$x\in\mathbb{R}^{n}$ is $\mathbb{R}[x]$. Every polynomial $c\in\mathbb{R}[x]$
may be uniquely described by a finite sum over multi-indices
$\alpha\in\mathbb{N}^{n}$ by $c(x)=\sum_{\alpha}c_{\alpha}x^{\alpha}$. The
degree of a polynomial $\deg{c}$ is equal to the largest exponent-sum
$\sum_{i=1}^{n}\alpha_{i}$ such that $c_{\alpha}\neq 0$. The set of
polynomials in $x$ with degree at most $2k$ is $\mathbb{R}[x]_{\leq 2k}$.

The set of continuous functions over a set $S$ is $C(S)$. Its subcone of
nonnegative continuous functions is $C_{+}(S)$. Given a product set $S\times
H$, the set of functions that are once-continuously differentiable in the
first variable $s$ and are twice-continuously differentiable in the second
variable $h$ is $C^{1,2}(S\times H)$. The set of signed Borel measures
supported in $S$ is $\mathcal{M}(S)$, and its subset of nonnegative Borel
measures supported in $S$ is $\mathcal{M}_{+}(S)$. The support of a measure
$\mu\in\mathcal{M}_{+}(\mathbb{R}^{n})$ is the smallest closed set of points
$s$ such that every open neighborhood $N_{\epsilon}(s)$ has positive measure
$\mu(N_{\epsilon}(s))>0$. The sets $C_{+}(S)$ and $\mathcal{M}_{+}(S)$ are in
topological duality when $S$ is compact, and they admit a duality pairing
$\langle\cdot,\cdot\rangle$ by Lebesgue integration: $\forall f\in C_{+}(S),\
\mu\in\mathcal{M}_{+}(S):\langle f,\mu\rangle=\int_{S}f(s)d\mu(s)$. This
duality pairing will be extended into a bilinear pairing between $C(S)$ and
$\mathcal{M}_{+}(S)$. The $\mu$-measure of a set $A\subseteq S$ may also be
expressed as the pairing of $\mu$ with the indicator function $I_{A}:\
\mu(A)=\langle I_{A},\mu\rangle=\inf\\{\langle w,\mu\rangle\;|\;w\in
C(S),w\geq I_{A}\\}$. The mass of $\mu$ is $\mu(S)=\langle 1,\mu\rangle$, and
$\mu$ is a probability distribution if this mass is 1. A vital example of a
probability distribution is the Dirac delta $\delta_{s^{\prime}}$ supported
only at $s^{\prime}$, obeying the point-evaluation pairing rule $\forall f\in
C(S):\ \langle f,\delta_{s^{\prime}}\rangle=f(s^{\prime}).$ The Lebesgue
(volume) distribution over a space $S$ is $\lambda_{S}$ with $\forall f\in
C(S):\ \langle f,\lambda_{S}\rangle=\int_{S}f(s)ds.$ The product measure
between $\mu\in\mathcal{M}_{+}(S_{1})$ and $\nu\in\mathcal{M}_{+}(S_{2})$ is
the unique measure $\mu\otimes\nu$ satisfying $\forall A_{1}\times
A_{2}\subseteq S_{1}\times S_{2}:\ \mu\otimes\nu(A_{1}\times
A_{2})=\mu(A_{1})\nu(A_{2}).$

###  2.2 Stochastic Processes

A stochastic process is a time-indexed set of random variables $\\{X_{t}\\}$
that are related together through system dynamics [27] (pushforward of an
initial distribution along flow maps). Properties of the stochastic process
will be analyzed in terms of test functions in an appropriate set
$\mathcal{C}$ (such as $\mathcal{C}=C([0,T]\times X)$. The expectation of a
test function $v(s,x)$ at the time $t$ according to the distribution $X_{t}$
is $\mathbb{E}[v(t,x)\mid X_{t}]$. Letting ${\Delta t}>0$ be a time step, the
generator $\mathcal{L}$ of a stochastic process satisfies (for all appropriate
test functions $v\in\textrm{dom}{\mathcal{L}}$)

$\mathcal{L}_{{\Delta t}}v(t,x)=\lim_{{{\Delta t}}^{\prime}\rightarrow{{\Delta
t}}}\frac{\mathbb{E}[v(t+{{\Delta t}}^{\prime},x)\mid X_{t+{{\Delta
t}}^{\prime}}]-v(t,x)}{{{\Delta t}}^{\prime}}.$ (3)

We will express the domain of $\mathcal{L}$ as
$\mathcal{C}=\textrm{dom}{\mathcal{L}}$, such that $\mathcal{C}$ is a subset
of the preimage of continuous functions under $\mathcal{L}$. The generator for
a discrete-time stochastic process is $\mathcal{L}_{{{\Delta t}}}$ with
${{\Delta t}}>0$, which is defined w.r.t. the test function class
$\mathcal{C}=C([0,T]\times X)$. For a discrete-time law following
$x_{t+{{\Delta t}}}=f(t,x_{t},\lambda_{t})$ in which the time-varying
parameter $\lambda_{t}\in\Lambda$ has a probability distribution of
$\xi(\lambda_{t})$, the associated generator satisfies

$\displaystyle\mathcal{L}_{{\Delta t}}v(t,x)=({{\Delta
t}}^{-1})\int_{\Lambda}v\left(t+{{\Delta
t}},f(t,x,\lambda)\right)d\xi(\lambda)-v(t,x).$ (4)

The generator for a continuous-time (Feller) stochastic process is
$\mathcal{L}_{{{\Delta t}}=0}$ for the class $\mathcal{C}=C^{1,2}([0,T]\times
X)$ [28]. The dynamical behavior of an Itô SDE has a description

$\displaystyle dx=f(t,x)dt+g(t,x)dW,$ (5) in which $W$ is the Wiener process.
The generator of the SDE in (5) is $\displaystyle\mathcal{L}_{0}v(t,x)=$
$\displaystyle\partial_{t}v(t,x)+f(t,x)\cdot\nabla_{x}v(t,x)$
$\displaystyle+g(t,x)^{T}\nabla^{2}_{xx}v(t,x)g(t,x)/2.$ (6)

A sequence of random variables $\\{Y_{t}\\}$ is a martingale if
$\mathbb{E}[Y_{t+\delta t}\mid\\{Y_{t^{\prime}}\\}_{t^{\prime}\leq t}]=Y_{t}$
[27] (in the sense of generalized conditional expectation). The generator
$\mathcal{L}$ solves a martingale problem for all possible time steps $s>0$
and test functions $v\in\mathcal{C}$ [29]. The martingale problem in discrete-
time (with $s\in\mathbb{N}$) is

$\displaystyle 0=\mathbb{E}[v(t+s,x)\mid X_{t+s}]$
$\displaystyle-\mathbb{E}[v(t,x)\mid X_{t}]$ (7)
$\displaystyle\textstyle\quad-\sum_{s^{\prime}=t}^{t+s}\mathbb{E}[\mathcal{L}_{{\Delta
t}}v(t,x)\mid X_{s^{\prime}}],$ and in continuous time (with
$s\in\mathbb{R}_{\geq 0}$) is $\displaystyle 0=\mathbb{E}[v(t+s,x)\mid
X_{t+s}]$ $\displaystyle-\mathbb{E}[v(t,x)\mid X_{t}]$ (8)
$\displaystyle\textstyle\quad-\int_{s^{\prime}=t}^{t+s}\mathbb{E}[\mathcal{L}_{0}v(t,x)\mid
X_{s^{\prime}}].$

For the rest of this paper, we will refer to the generator $\mathcal{L}$
interchangeably with its stochastic process.

The martingale problem can be described using the theory of occupation
measures. Given an initial condition (point) $x_{0}\in X_{0}\subseteq X$ and a
time $t\in[0,T]$, we refer to $x(t\mid x_{0})$ as the trajectory of
$\\{X_{t}\\}$ with generator $\mathcal{L}$ starting at the initial
distribution $\delta_{x=x_{0}}$. The occupation measure
$\mu\in\mathcal{M}_{+}([0,T]\times X)$ and terminal measure
$\mu_{\tau}\in\mathcal{M}_{+}([0,T]\times X)$ of the stochastic process
$\\{X_{t}\\}$ from a starting time $t_{0}$ up to a stopping time $t^{*}$
satisfies $\forall A\subseteq[t_{0},T],B\subseteq X$:

$\displaystyle\mu(A\times B)$
$\displaystyle=\int_{X_{0}}\int_{t=t_{0}}^{t^{*}}I_{A\times B}\left(t,x(t\mid
x_{0})\right)dt\,d\mu_{0}(x_{0})$ (9a) $\displaystyle\mu_{\tau}(A\times B)$
$\displaystyle=\int_{X_{0}}I_{A\times B}\left(t^{*},x(t^{*}\mid
x_{0})\right)\,d\mu_{0}(x_{0}).$ (9b)

The measures in (9) may also be defined with respect to a probability
distribution of stopping times $t^{*}\in\mathcal{M}_{+}([0,T])$.

A tuple of measures
$(\mu_{0},\mu,\mu_{\tau})\in\mathcal{M}_{+}(X_{0})\times\mathcal{M}_{+}([0,T]\times
X)^{2}$ obeys the martingale relations (7) or (8) of the stochastic process
with generator $\mathcal{L}$

$\displaystyle\forall v\in\mathcal{C}:$ $\displaystyle\langle
v,\mu_{\tau}\rangle$ $\displaystyle=\langle
v(t_{0},\cdot),\mu_{0}\rangle+\langle\mathcal{L}v,\mu\rangle.$ (10)

Relation (10) is called Liouville equation for Ordinary Differential Equations
and is referred to as Dynkin’s equation in the context of SDEs. We will also
equivalently express (10) in a shorthand form by

$\displaystyle\mu_{\tau}=\delta_{t_{0}}\otimes\mu_{0}+\mathcal{L}^{\dagger}\mu.$
(11)

###  2.3 Sum-of-Squares Methods

Let $c\in\mathbb{R}[x]_{\leq 2k}$ be a polynomial. The polynomial $c$ is
nonnegative if $\forall x\in\mathbb{R}^{n}:c(x)\geq 0$. Verification of
polynomial nonnegativity is generically an NP-hard task. One tractable method
to construct nonnegative polynomials is through SOS representations. The
polynomial $c$ is SOS if there exists $N$ polynomials $q_{j}\in\mathbb{R}[x]$
such that $c(x)=\sum_{j=1}^{N}q_{j}^{2}(x)$. The cone of SOS polynomials is
$\Sigma[x]$, and its subset of SOS polynomials with degree $\leq 2k$ is
$\Sigma[x]_{\leq 2k}.$ The set of SOS polynomials equals the cone of
nonnegative polynomials only when $n=1,2k=2,$ or $(n=2,2k=4)$ [30]. SOS
polynomials become vanishingly small in the set of nonnegative polynomials as
$n$ and $d$ rise [31].

SOS representations up to fixed degree can be conducted by formulating SDPs.
To each polynomial $c\in\mathbb{R}[x]_{\leq 2k}$, there exists a (typically
nonunique) associated polynomial vector $v_{k}(x)\in(\mathbb{R}[x]_{\leq
d})^{s}$ and a symmetric Gram matrix $Q\in\mathbb{S}^{s}$ such that
$c(x)=v_{k}(x)^{T}Qv_{k}(x)$. One such choice of vectors $v_{k}(x)$ is the set
of monomials in $x$ between degrees $0..k$, in which case $s=\binom{n+k}{k}$.
A polynomial $c$ is SOS if its Gram matrix $Q$ is also Positive Semidefinite
(PSD) [32].

BSA Basic Semialgebraic (BSA) set $\mathbb{K}$ is a set described by a finite
number of bounded-degree polynomial inequality $(N_{g})$ and equality
constraints $(N_{h})$:

$\displaystyle\mathbb{K}=\\{x\in\mathbb{R}^{n}\mid g_{i}(x)\geq 0,\
h_{j}(x)=0,\forall i\in 1..N_{g},j\in 1..N_{h}\\}.$ (12)

The class of Weighted Sum of Squares (WSOS) polynomials over $\mathbb{K}$ in
(12) is the set $\Sigma[\mathbb{K}]$ of polynomials $c$ that admit a
description of

$\displaystyle c(x)=$
$\displaystyle\sigma_{0}(x)+\textstyle\sum_{i}{\sigma_{i}(x)g_{i}(x)}+\textstyle\sum_{j}{\phi_{j}(x)h_{j}(x)}$
(13a)
$\displaystyle\sigma_{0},\sigma_{i}\in\Sigma[x],\quad\phi_{j}\in\mathbb{R}[x].$
(13b)

The set $\mathbb{K}$ is compact if there exists a finite $R>0$ such that
$\\{x\mid\lVert x\rVert_{2}^{2}\leq R\\}\supseteq\mathbb{K}$. The set is
additionally Archimedean if $R-\lVert x\rVert_{2}^{2}\in\Sigma[\mathbb{K}].$
Not every compact set is Archimedean [33], but if an $R$ verifying compactness
of $\mathbb{K}$ is known, then the constraint $R-\lVert x\rVert_{2}^{2}\geq 0$
can be added to the description of $\mathbb{K}$ to render $\mathbb{K}$
Archimedian. Every compact polytope and ellipsoid is also Archimedean [14].

The Putinar Positivestellensatz states that every positive polynomial over
$\mathbb{K}$ is also a member of $\Sigma[\mathbb{K}]$ when $\mathbb{K}$ is
Archimedean [34]. However, the polynomial degrees of the multipliers
$(\sigma,\phi)$ needed to certify nonnegativity of such a positive $c$ could
be exponential in $n$ and $k$ [35].

##  3 Unsafe Linear Programs

This section will present LPs to analyze safety of a stochastic process with
generator $\mathcal{L}$.

###  3.1 Assumptions

We will posit the following assumptions throughout this paper:

  * A1

The state-sets $X_{0},X_{u},X$ are all compact.

  * A2

The time interval $[t_{0},T]$ is compact.

  * A3

Trajectories stop upon the first exit from $X$ ($\tau_{X}\wedge T$).

  * A4

The test function set $\mathcal{C}=\textrm{dom}(\mathcal{L})$ satisfies
$\mathcal{C}\subseteq C([t_{0},T]\times X)$ with $1\in\mathcal{C}$ and
$\mathcal{L}1=0$.

  * A5

The set $\mathcal{C}$ separates points and is multiplicatively closed.

  * A6

There exists a countable set $\\{v_{k}\\}\in\mathcal{C}$ such that $\forall
v\in\mathcal{C}:$ $(v,\mathcal{L}v)$ is contained in the bounded pointwise
closure of the linear span of $\\{(v_{k},\mathcal{L}v_{k})\\}$.

Assumptions A4-A6 are the same as in Condition 1 of [6]. We note that
generators for Markov, SDE, and Lévy processes evolving in a compact set
satisfy A1-A6.

###  3.2 Measure Program

We will pose a convex but infinite-dimensional LP in order to upper-bound
program (2). Such an LP will be posed in terms of a relaxed occupation measure
$(\mu_{0},\mu_{\tau},\mu)$.

######  Theorem 3.1.

The following LP will upper-bound (2) with $p^{*}(t_{0},X_{0})\geq
P^{*}(t_{0},X_{0})$ for each $X_{0}\subseteq X$

$\displaystyle p^{*}(t_{0},X_{0})=$ $\displaystyle\sup\quad\langle
I_{X_{u}},\mu_{\tau}\rangle$ (14a)
$\displaystyle\mu_{\tau}=\delta_{t_{0}}\otimes\mu_{0}+\mathcal{L}^{\dagger}\mu$
(14b) $\displaystyle\langle 1,\mu_{0}\rangle=1$ (14c)
$\displaystyle\mu_{0}\in\mathcal{M}_{+}(X_{0})$ (14d) $\displaystyle\mu,\
\mu_{\tau}\in\mathcal{M}_{+}([t_{0},T]\times X).$ (14e)

###### Proof.

Let $x_{0}\in X_{0}$ be an initial condition and $t^{*}\in[t_{0},T]$ be a
stopping time. A set of measures $(\mu_{0},\mu,\mu_{\tau})$ may be constructed
to fulfill constraints (14b)-(14e). To begin, the initial measure $\mu_{0}$
may be chosen as the Dirac delta $\delta_{x=x_{0}}$. The measure $\mu$ can be
set to the occupation measure of the stochastic trajectory $x(t\mid x_{0})$
between the initial time of ${t_{0}}$ and the stopping time of
$t^{*}\wedge\tau_{X}$. The terminal measure $\mu_{\tau}(t,x)$ may then be
picked as the the time-state distribution at time $t^{*}$ (in which
trajectories are stopped according to $t^{*}\wedge\tau_{X}$). Because there
exists a mapping from every feasible point $(t^{*},x_{0})$ of (2) to a set of
measures in (16b)-(16f), it holds that $p^{*}_{\tau}({t_{0}},X_{0})\geq
P^{*}({t_{0}},X_{0})$. ∎

The function $I_{X_{u}}(x)$ in the objective (14a) is a discontinuous function
of $x$. We can reformulate (14) as an LP where all terms are continuous in
$(t,x)$ (under A1-A6) through the use of the following restriction lemma:

######  Lemma 3.2 (Theorem 3.1 of [36]).

Let $S_{1},S_{2}$ be spaces with $S_{1}\subseteq S_{2}$, and let
$\xi\in\mathcal{M}_{+}(S_{2})$ be a measure. The restriction of $\xi$ to
$S_{1}$ is uniquely given by the measure $\nu^{*}$ solving the following
optimization problem

$\displaystyle Q^{*}=\ $ $\displaystyle\sup_{\nu,\hat{\nu}}\ \langle
1,\nu\rangle$ (15a) $\displaystyle\xi=\nu+\hat{\nu}$ (15b)
$\displaystyle\nu\in\mathcal{M}_{+}(S_{1})$ (15c)
$\displaystyle\hat{\nu}\in\mathcal{M}_{+}(S_{2}).$ (15d)

The objective $Q^{*}$ will satisfy $Q^{*}=\langle 1,\nu^{*}\rangle=\langle
I_{S_{1}},\xi\rangle$. The measure $\hat{\nu}^{*}$ will be supported on
$\textbf{closure}(S_{2}\setminus S_{1})$.

We will define a peak measure $\mu_{p}\in\mathcal{M}_{+}([t_{0},T]\times
X_{u})$ and a complement measure $\mu_{c}\in\mathcal{M}_{+}([t_{0},T]\times
X)$ in order to split the terminal measure $\mu_{\tau}$ as
$\mu_{\tau}=\mu_{p}+\mu_{c}$. This decomposition (through Lemma 3.2) will
allow for the definition of an LP with continuous coefficients in $(t,x)$.
This continuity provides strong duality (under A1-A6) and convergence of the
moment-SOS hierarchy to $P^{*}(t_{0},X_{0})$ from (2) (under A1-A6 and
assumptions on polynomial structure). The LP in which $\mu_{\tau}$ is
decomposed as $\mu_{p}+\mu_{c}$ is:

######  Theorem 3.3.

The following LP has the same objective value $p^{*}({t_{0}},X_{0})$ as in
(2):

$\displaystyle p^{*}({t_{0}},X_{0})=$ $\displaystyle\sup\quad\langle
1,\mu_{p}\rangle$ (16a)
$\displaystyle\mu_{p}+\mu_{c}=\delta_{t_{0}}\otimes\mu_{0}+\mathcal{L}^{\dagger}\mu$
(16b) $\displaystyle\langle 1,\mu_{0}\rangle=1$ (16c)
$\displaystyle\mu_{0}\in\mathcal{M}_{+}(X_{0})$ (16d) $\displaystyle\mu,\
\mu_{c}\in\mathcal{M}_{+}([t_{0},T]\times X)$ (16e)
$\displaystyle\mu_{p}\in\mathcal{M}_{+}([t_{0},T]\times X_{u}).$ (16f)

###### Proof.

1) Let $(\mu_{0},\mu_{p},\mu_{c},\mu)$ be a feasible point to (16b)-(16f).
Then $(\mu_{0},\mu_{p}+\mu_{c},\mu)$ is a feasible point to (14b)-(14e).

2) Let $(\mu_{0},\mu_{\tau},\mu)$ be a feasible point to the constraints
(14b)-(14e). Unique measures $\mu_{p},\ \mu_{c}$ may be defined according to
the restriction Lemma 3.2 with $S_{1}=[t_{0},T]\times X_{u}$ and
$S_{2}=[t_{0},T]\times X$ by solving

$\displaystyle\sup\ \langle 1,\mu_{\tau}\rangle$ (17a)
$\displaystyle\mu_{\tau}=\mu_{p}+\mu_{c}$ (17b)
$\displaystyle\mu_{c}\in\mathcal{M}_{+}([t_{0},T]\times X)$ (17c)
$\displaystyle\mu_{p}\in\mathcal{M}_{+}([t_{0},T]\times X_{u}),$ (17d)

with $\langle I_{X_{u}},\mu_{\tau}\rangle=\langle 1,\mu_{p}\rangle$.

There exists a mapping relating each feasible point of (14) to (16) and vice
versa while keeping the same objective, which proves that the two objectives
are equal.

∎

######  Theorem 3.4.

Under assumptions A1-A6, the objectives of (16) and (1b) are equal
$(p^{*}({t_{0}},X_{0})=P^{*}({t_{0}},X_{0}))$.

###### Proof.

By Theorem 3.1 of [6] (and under Assumptions A1-A6), every tuple
$(\mu_{0},\mu,\mu_{\tau})$ with $\mu_{0}\in\mathcal{M}_{+}([t_{0},T]\times X)$
(with $\langle 1,\mu_{0}\rangle=1$) and
$\mu,\mu_{\tau}\in\mathcal{M}_{+}([t_{0},T]\times X)$ satisfying a Martingale
relation (11) is supported on the graph of a stochastic process (1b).

The tuple $(\mu_{0},\mu,\mu_{p}+\mu_{c})$ from the Martingale relation (11) is
supported on the graph of a stochastic process. The measures $\mu_{p}$ and
$\mu_{c}$ can be separated according to Lemma 3.2. This concludes the no-
relaxation-gap proof. ∎

###  3.3 Function Program

######  Theorem 3.5.

LP LP in terms of an auxiliary function $v$ that forms a weak dual
($p^{*}({t_{0}},X_{0})\geq d^{*}({t_{0}},X_{0})$) to (16) is

$\displaystyle d^{*}({t_{0}},X_{0})=$
$\displaystyle\inf_{\gamma\in\mathbb{R}}\quad\gamma$ (18a)
$\displaystyle\gamma\geq v(t_{0},x)$ $\displaystyle\forall x\in X_{0}$ (18b)
$\displaystyle\mathcal{L}v(t,x)\leq 0$
$\displaystyle\forall(t,x)\in[t_{0},T]\times X$ (18c) $\displaystyle
v(t,x)\geq 0$ $\displaystyle\forall(t,x)\in[t_{0},T]\times X$ (18d)
$\displaystyle v(t,x)\geq 1$ $\displaystyle\forall(t,x)\in[t_{0},T]\times
X_{u}$ (18e) $\displaystyle v\in\mathcal{C}.$ (18f)

Strong duality ($p^{*}(t_{0},X_{0})=d^{*}(t_{0},X_{0})$) will hold under
assumptions A1-A4.

###### Proof.

See Appendix A. ∎

######  Remark 1.

Constraints (18b)-(18f) have the same form as a stochastic (time-dependent)
barrier function at a fixed probability level $\gamma$ (from [37]). Our work
involves application of $\gamma$ as an optimization variable, as well as a
proof of no relaxation gap in Theorem 3.4 under A1-A6.

##  4 Risk Contour Linear Programs

The LPs in (16) and (18) determine the worst-case probability of unsafety
among stochastic trajectories starting from an initial set $X_{0}$. The
optimal value of (16) and (18) is therefore equal to 1 when $X_{0}=X$ and
$X_{u}$ is nonempty (because then $\varnothing\neq X_{u}\subset X_{0}$ and
there exists trajectories starting in the unsafe set). This section will form
a continuous-function LP similar to (18), whose function values evaluated at a
given initial point $X_{0}$ form a upper-bound on the probability of unsafety
when starting at $x_{0}$. These upper-bounds will be proven to converge in an
$L_{1}$ sense over $X$ to the true risk function $P^{*}$. Such upper-bounds on
unsafety allow for interpretation of risk when performing path planning tasks.

Let $\mu_{0}\in\mathcal{M}_{+}(X)$ be a fixed probability distribution
$(\langle 1,\mu_{0}\rangle=1)$. The $\mu_{0}$-averaged probability of unsafety
is

$\displaystyle M^{*}({t_{0}},\mu_{0})=$
$\displaystyle\int_{X}P^{*}({t_{0}},x_{0})\ d\mu_{0}(x_{0}).$ (19)

One specific case we will focus on in the examples is where $\mu_{0}$ is the
the uniform distribution over $X$ $(\mu_{0}=\lambda_{X}/\textrm{vol}(X)).$

###  4.1 Measure Program

######  Corollary 1.

LP LP to upper-bound (2) with $m^{*}({t_{0}},\mu_{0})\geq
M^{*}({t_{0}},\mu_{0})$ in general, and
$m^{*}({t_{0}},\mu_{0})=M^{*}({t_{0}},\mu_{0})$ under A1-A6, is

$\displaystyle m^{*}({t_{0}},\mu_{0})=$ $\displaystyle\sup\quad\langle
1,\mu_{p}\rangle$ (20a)
$\displaystyle\mu_{p}+\mu_{c}=\delta_{t_{0}}\otimes\mu_{0}+\mathcal{L}^{\dagger}\mu$
(20b) $\displaystyle\mu,\ \mu_{c}\in\mathcal{M}_{+}([t_{0},T]\times X)$ (20c)
$\displaystyle\mu_{p}\in\mathcal{M}_{+}([t_{0},T]\times X_{u}).$ (20d)

###### Proof.

These optimality properties follow from arguments used in Theorem 3.3 (upper-
bound by a measure construction starting from the distribution $\mu_{0}$) and
Theorem 3.4 (restriction measures). The measures $(\mu,\mu_{p}+\mu_{c})$ are
supported on the graph of the stochastic process $\mathcal{L}$ starting from
$\mu_{0}$ for every feasible $(\mu,\mu_{p},\mu_{c})$. ∎

###  4.2 Function Program

For given $v\in\mathcal{C}$ and $\mu_{0}\in\mathcal{M}_{+}(X)$, let
$v(t_{0},\bullet)$ be shorthand notation for the function $x\mapsto
v(t_{0},x)$ such that $\langle v(t_{0},\bullet),\mu_{0}\rangle=\langle
v(t_{0},x),\mu_{0}(x)\rangle$. We will modify terms (18a)-(18b) in order to
form LP LP penalized over the whole space $X$ with an objective integrating
against the weighting distribution $\mu_{0}$:

$\displaystyle J^{*}(t_{0},\mu_{0})=$ $\displaystyle\inf\langle
v(t_{0},\bullet),\mu_{0}\rangle\\!\\!\\!\\!$ (21a)
$\displaystyle\mathcal{L}v(t,x)\leq 0$
$\displaystyle\forall(t,x)\in[t_{0},T]\times X$ (21b) $\displaystyle
v(t,x)\geq 0$ $\displaystyle\forall(t,x)\in[t_{0},T]\times X$ (21c)
$\displaystyle v(t,x)\geq 1$ $\displaystyle\forall(t,x)\in[t_{0},T]\times
X_{u}$ (21d) $\displaystyle v\in\mathcal{C}.$ (21e)

######  Lemma 4.1.

The objective (21a) from (21) is upper-bounded by $J^{*}(t_{0},\mu_{0})\leq
1$.

###### Proof.

The function $v(t,x)=1$ is feasible for all constraints (21b)-(21e). Given
that $\mu_{0}$ is a probability measure, the objective in (21a) will satisfy
$J^{*}\leq\langle 1,\mu_{0}\rangle=1$. ∎

We now show that any solution of (21b)-(21e) returns an upper-bound of the
probability of unsafety.

######  Theorem 4.2.

Let $v$ be feasible for (21b)-(21e) and let assumptions A1-A6 hold. Then

$v(t_{0},x_{0})\geq P(t_{0},x_{0})$ (22)

for every initial condition $x_{0}\in X$.

###### Proof.

See Appendix B. ∎

######  Remark 2.

Level sets of any $v(t,x)$ solving (21b)-(21e) can therefore be used to upper-
bound the risk of stochastic execution for $\mathcal{L}$ when starting at any
$(t,x_{0})\in[t_{0},T]\times X$.

The notation $p(t_{0},x_{0})$ will be subsequently used to denote the function
obtained by solving (16) as $x_{0}\mapsto p(t_{0},X_{0}=\\{x_{0}\\})$.

######  Theorem 4.3.

Programs (20) and (21) obey strong duality
$m^{*}(t_{0},\mu_{0})=J^{*}(t_{0},\mu_{0})$ under A1-A6.

###### Proof.

This strong duality holds by slight modification of the method used in Theorem
3.5 and Appendix A. ∎

######  Theorem 4.4.

Let $\\{v_{k}\\}_{k\geq 1}$ be a sequence of solutions to (21b)-(21e) with
$\langle v_{k}(t_{0},\bullet),\mu_{0}\rangle\rightarrow J^{*}(t_{0},\mu_{0})$
(under A1-A6). Then $v_{k}$ converges in $L_{1}(\mu_{0})$ to $p^{*}$ as in

$\displaystyle\int_{X}\left(\vphantom{\sum}v_{k}(t_{0},x_{0})-p^{*}(t_{0},x_{0})\right)d\mu_{0}(x_{0})\underset{k\to\infty}{\longrightarrow}0.$
(23)

###### Proof.

See Appendix C.

∎

##  5 Unsafe Semidefinite Programs

We will approximate programs (18) and (21) through the moment-SOS hierarchy of
SDPs, as reviewed in Section 2.3. In order to apply SOS methods towards
convergent approximation, we require additional assumptions:

  * A7

The sets $X_{0},\ X_{u},\ X$ are all Archimedean BSA sets.

  * A8

The generator $\mathcal{L}$ is closed under polynomials
($v(t,x)\in\mathbb{R}[t,x]\implies\mathcal{L}v(t,x)\in\mathbb{R}[t,x]$).

Given a degree $k$ and a generator $\mathcal{L}$ obeying Assumption A8, we
define the dynamics degree $\tilde{k}$ as

$\displaystyle\tilde{k}=\lceil\max_{v\in\mathbb{R}_{\leq
2k}[t,x]}\deg(\mathcal{L}v(t,x))/2\rceil.$ (24)

###  5.1 Unsafe-Probability SDP

The degree-$k$ SOS tightening of program (18) is

$\displaystyle d^{*}_{k}(t_{0},X_{0})=$
$\displaystyle\inf_{\gamma\in\mathbb{R}}\quad\gamma$ (25a)
$\displaystyle\gamma-v(t_{0},x)\in\Sigma[X_{0}]_{\leq 2k}$ (25b)
$\displaystyle-\mathcal{L}v(t,x)\in\Sigma[[0,T]\times X]_{\leq 2\tilde{k}}$
(25c) $\displaystyle v(t,x)\in\Sigma[[0,T]\times X]_{\leq 2k}$ (25d)
$\displaystyle v(t,x)-1\in\Sigma[[0,T]\times X_{u}]_{\leq 2k}.$ (25e)

######  Remark 3.

Constraint (25d) absorbs the polynomial restriction of
$v\in\mathbb{R}[t,x]_{\leq 2k}$.

We require the following Lemma to ensure convergence of (25) to (18) as
$k\rightarrow\infty$.

######  Lemma 5.1.

Under Assumptions A1-A6, all feasible measures
$\boldsymbol{\mu}=(\mu_{0},\mu,\mu_{c},\mu_{\tau})$ in (16) are bounded.

###### Proof.

Boundedness of nonnegative measures will be proven by the sufficient condition
of compact support and finite mass. Compact support of measures in
$\boldsymbol{\mu}$ is ensured by A1. The initial measure $\mu_{0}$ is a
probability distribution $(\langle 1,\mu_{0}\rangle=1)$ by constraint (16c).
The sum $\mu_{c}+\mu_{\tau}$ likewise has mass 1 by the Liouville constraint
(16b), when passing in the test function $v=1$. Given that $\mu_{c}$ and
$\mu_{\tau}$ are each nonnegative measures, it holds that they both have
finite mass (upper-bounded by 1). Lastly, assignment of $v=t$ to (16b) with
$\mathcal{L}t=1$ results in $\langle 1,\mu\rangle=\langle
t,\mu_{c}+\mu_{\tau}\rangle\leq T$. All measures have bounded masses and
compact supports, and therefore are bounded. ∎

######  Theorem 5.2.

Under assumptions A1-A8, program (25) will converge to (18) with
$\lim_{k\rightarrow\infty}d^{*}_{k}(t_{0},X_{0})=d^{*}(t_{0},X_{0}).$

###### Proof.

This convergence will hold by Corollary 8 of [38], because all sets are
Archimedean (A7), dynamics are polynomial (A8), measures are bounded (Lemma
5.1), and the objective value of (25) is finite (bounded below by 0). ∎

###  5.2 Risk SDP

The degree-$k$ SOS tightening for the risk-contour program in (21) is

$\displaystyle J^{*}_{k}(t_{0},\mu_{0})=$
$\displaystyle\inf\int_{X}v(t_{0},x)d\mu_{0}(x)$ (26a)
$\displaystyle-\mathcal{L}v(t,x)\in\Sigma[[0,T]\times X]_{\leq 2\tilde{k}}$
(26b) $\displaystyle v(t,x)\in\Sigma[[0,T]\times X]_{\leq 2k}$ (26c)
$\displaystyle v(t,x)-1\in\Sigma[[0,T]\times X_{u}]_{\leq 2k}.$ (26d)

######  Remark 4.

The function $v(t,x)=1$ is feasible for constraints (26b)-(26d) at every
degree $k\in\mathbb{N}$. As a result, the objective in (26) is always upper-
bounded by $J_{k}^{*}(t_{0},\mu_{0})\leq 1$.

######  Theorem 5.3.

Under assumptions A1-A8, the risk contour program (26) will converge to (21)
with (18) with $\lim_{k\rightarrow\infty}J^{*}_{k}=J^{*}.$

###### Proof.

The proof of this theorem follows the same steps as in Theorem 5.2 with
respect to the measure program in (20). ∎

Let $v_{k}(t,x)\in\Sigma[[0,T]\times X]_{\leq 2k}$ be the solution to (26) at
degree-$k$, and let $I_{u}$ be the $0/1$ indicator function on the unsafe set
$X_{u}$. Then the probability of unsafety when starting at a point $x_{0}\in
X$ is upper-bounded by

$\displaystyle q_{1:k}(x)=\min(1,\min_{k^{\prime}\in
1..k}v_{k^{\prime}}(t_{0},x)).$ (27)

######  Corollary 2.

The sequence of functions $q_{1:k}(x)$ in increasing $k$ will converge in
measure $\mu_{0}$ to $v(t_{0},x)$.

###### Proof.

This corollary follows from Theorem 4.4, in which the sequence
$\\{v_{k^{\prime}}\\}$ is used to approximate $v^{*}$. The minimization among
all $k^{\prime}$ in (27) further sharpens the estimate of $v^{*}$. ∎

###  5.3 Computational Complexity

We will quantify computational complexity of the degree-$k$ tightenings of
(25) and (26) by the size of the maximal Gram matrices involved in their SOS
program. In the typical case where $\tilde{k}>k$ (only violated under A1-A8
when $\mathcal{L}$ maps every polynomial to a constant), the largest Gram
matrix will occur in the Lie constraints (25c) and (26b). The Lie constraints
each have $n+1$ variables $(t,x)$, so the Gram matrix size when using the
monomial basis is $\binom{n+1+\tilde{k}}{\tilde{k}}.$ All other constraints
have a lower degree ($k$ rather than $\tilde{k}$), or are posed only over the
$n$ variables $x$. The complexity of using an interior-point method to solve
the SOS programs will therefore scale based on $O((n+1)^{6d})$ for fixed $d$
or $O(d^{(n+1)})$ for fixed $n$ [14, 39].

##  6 Numerical Examples

MATLAB (2021a) code to reproduce all examples is available at
https://github.com/jarmill/prob_unsafe. All programs are modeled using YALMIP
[40] and solved using Mosek 10 [41]. All examples will involve an initial time
of $t_{0}=0$.

###  6.1 Two-Dimensional Cubic SDE

Our first demonstration analyzes safety of a cubic polynomial SDE from Example
1 of [10]:

$dx=\begin{bmatrix}x_{2}\\\
-x_{1}-x_{2}-\frac{1}{2}x^{3}_{1}\end{bmatrix}dt+\begin{bmatrix}0\\\
0.1\end{bmatrix}dW.$ (28)

Safety of (28) is evaluated within the state space of $X=[-2,2]^{2}$ until a
time horizon of $T=5$. The unsafe set is a moon-shaped region
$X_{u}=\\{x\in\mathbb{R}^{2}\mid
0.9083^{2}\leq(x_{1}+0.5006)^{2}+(x_{2}+0.2902)^{2},0.5^{2}\geq(x_{1}-0.2)^{2}+x_{2}^{2}\\}.$
The initial set $X_{0}$ is a circle of radius $R_{0}$ and center
$[0.85;-0.75]$. Figure 2(a) plots trajectories of (28) starting from $X_{0}$
(magenta region with $R_{0}=0.2$). Note how some sampled trajectories touch
the right corner of the red half-circle $X_{u}$ (and are therefore unsafe).

Figure 1:  Trajectories of SDE (28) (cyan) initialized in a disk $X_{0}$
(purple boundary), hitting an unsafe region $X_{u}$ (red moon).

Table 1 reports the probability of unsafety computed by program (25) for a
circular initial set of radius $R_{0}$.

Table 1: Unsafe probability upper-bounds for system (28)  order | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---
$R_{0}=0$ | 1 | 0.9442 | 0.6072 | 0.4805 | 0.4068 | 0.3696
$R_{0}=0.2$ | 1 | 0.9736 | 0.7943 | 0.7091 | 0.6390 | 0.6132

Figure 2(a) plots the unsafety upper-bounding function $\min(1,v(0,x))$ at
$T=5$ acquired by solving (25) at order $k=$6 with $R_{0}=0.2$. The unsafe set
$X_{u}$ is the red moon. The magenta circle is the boundary of the initial set
$X_{0}$. Note how the probability estimate is sharper in the region
surrounding $X_{0}$, as compared to the sea of $\textrm{Prob}=1$ saturation
away from $X_{0}$.

In contrast, Figure 2(b) solves the order-6 SOS tightening of program (21) to
produce a risk map $v$ valid in $X=[-2,2]^{2}$. This risk map results in
$v(0,[0.85;-0.75])=0.4366,$ which is looser than the order-6 probability bound
of $0.3696$ from Table 1. However, Figure 2(b) produces an interpretable
visualization of risk across the entire set $X$.

(a) solution of (25)

(b) solution of (26)

Figure 2: Risk levels of the unsafety upper-bound function $v(0,x)$ at $t=T=5$
and order $k=$6 for SDE (28) initialized in a disk $X_{0}$ (purple boundary)
with unsafe region $X_{u}$ (red moon).

###  6.2 Adjustment of Time Horizons

In the second experiment, we consider a linear SDE with dynamics

$\displaystyle dx=[-x_{2};x_{1}]dt+[0;0.1]dW$ (29)

The unsafe set for this example is a half-circular region with $X_{u}=\\{x\mid
0.5^{2}\geq(x+0.5)^{2}+(x_{2}+0.75)^{2},x_{2}\geq-0.75\\}$. The plots in
Figure 3 are obtained by solving (26) at order $k=6$ for time horizons of
$T\in\\{1,3,5\\}$ and a state set of $X=[-2,2]^{2}$. As the time horizon $T$
increases, so too does the area of $X$ with high probabilities of unsafety.
The growth of the high-probability regions occurs approximately in a clockwise
pattern from $X_{u}$, matching the interpretation of (29) as simple harmonic
motion perturbed by diffusion in the vertical coordinate.

(a) $T=1$

(b) $T=3$

(c) $T=5$

Figure 3: Risk contours for the linear SDE in (29) in increasing time horizon
$T$ w.r.t. unsafe region $X_{u}$ (red half-disk)

###  6.3 Discrete-time example

This subsection will focus on a discrete-time stochastic process involving a
time-step of $\tau=1$. The law of this stochastic process with parameter
$\lambda\in\mathbb{R}$ is

$\displaystyle x_{+}=\begin{bmatrix}-0.3x_{1}+0.8x_{2}+x_{1}x_{2}\lambda/4\\\
-0.9x_{1}--0.1x_{2}-0.2x_{1}^{2}+\lambda/40\end{bmatrix},$ (30)

in which $\lambda$ is i.i.d. sampled according to a unit-normal distribution
at each time $(\lambda_{t}\sim\mathcal{N}(0,1))$.

We evaluate the safety of (30) with respect to the state set
$X=[-1.5,1.5]^{2}$, the time horizon of $T=10$, and the half-circle unsafe set
of $X_{u}=\\{x\mid 0.4^{2}\geq(x_{1}-0.8)^{2}+(x_{2}-0.2)^{2},\
x_{1}+x_{2}\geq 1\\}$. Stochastic trajectories of (30) evolve starting at an
circular initial set with radius $R_{0}$ and center $[-1;0]$. Figure 4 plots
5,000 sampled trajectories (blue dots) with respect to the magenta circular
initial set and the red half-circle unsafe set. Some of the sampled trajectory
points fall inside the red half-circle, corresponding to unsafety when
beginning in $X_{0}$.

Figure 4:  Trajectories of discrete-time (30) (cyan) initialized in a disk
$X_{0}$ (purple boundary), hitting an unsafe region $X_{u}$ (red half disk).

Table 2 reports probabilities of unsafety for (30) found by solving (25) at
initial radii $R_{0}=0$ and $R_{0}=0.4$. To improve numerical conditioning, we
normalize the time-steps from $(\tau,T)=(1,10)$ to $(\tau,T)=(0.1,1)$ without
affecting the autonomous dynamics in (30).

Table 2: Unsafe probability upper-bounds for system (30)  order | 1 | 2 | 3 | 4 | 5 | 6
---|---|---|---|---|---|---
$R_{0}=0$ | 1.00 | 1.00 | 0.1569 | 0.0103 | 1.871$\mathbf{e}$-3 | 7.052$\mathbf{e}$-4
$R_{0}=0.4$ | 1.00 | 1.00 | 0.9801 | 0.7054 | 0.5225 | 0.4017

Figure 5(a) plots risk contours found by solving (25) at order 6 for
$R_{0}=0.4$. Figure 5(b) plots risk contours $v(0,x)$ acquired from (26) at
degree $k=$6 and time $t=T=10$. The returned risk map has an evaluation of
$v(0,[-1;0])=0.4915$, as compared to the (25) estimate of $7.052\times
10^{-4}$ at $R_{0}=0$.

(a) solution of (25)

(b) solution of (26)

Figure 5: Risk levels of unsafety upper-bound function $v(0,x)$ at $t=T=10$
and order $k=$6 for SDE (30) initialized in a disk $X_{0}$ (purple boundary)
with unsafe region $X_{u}$ (red moon).

##  7 Conclusion

This paper presents a method to analyze the probability of unsafety for
stochastic processes by forming LPs in occupation measures. These LPs are
upper-bounded by the moment-SOS hierarchy of SDPs, yielding a convergent
sequence of bounds to the true probability of unsafety. Modification of the
objective leads to the development of visually interpretable risk contours for
use in analysis and motion planning. The presented method can be used for any
stochastic process satisfying assumptions A1-A8 (compactness, described by
generators, having a polynomial description).

Future work involves finding redundant constraints to refine and improve the
unsafe-probability bounds. Other work includes extension to Conditional Value-
at-Risk-minimizing stochastic control methods.

## Acknowledgements

The authors would like to thank Roy S. Smith, Niklas Schmid, and the Automatic
Control Laboratory at ETH Zürich for their advice and support.

##  Appendix A Strong Duality

This appendix will prove strong duality between (16) and (18) by using
arguments from Theorem 2.6 of [42].

We collect the variables in (16) and (18) respectively into,

$\displaystyle\boldsymbol{\mu}$ $\displaystyle=(\mu_{0},\mu,\mu_{c},\mu_{p})$
$\displaystyle\boldsymbol{\ell}$ $\displaystyle=(\gamma,v).$ (31)

Variable spaces related to $\boldsymbol{\mu}$ are,

$\displaystyle\mathcal{X}^{\prime}$ $\displaystyle=C(X_{0})\times
C([0,T]\times X)^{2}\times C([0,T]\times X_{u})$ (32)
$\displaystyle\mathcal{X}$
$\displaystyle=\mathcal{M}(X_{0})\times\mathcal{M}([0,T]\times
X)^{2}\times\mathcal{M}([0,T]\times X_{u}),$

with nonnegative subcones of

$\displaystyle\mathcal{X}_{+}^{\prime}$ $\displaystyle=C_{+}(X_{0})\times
C_{+}([0,T]\times X)^{2}\times C_{+}([0,T]\times X_{u})$ (33)
$\displaystyle\mathcal{X}_{+}$
$\displaystyle=\mathcal{M}_{+}(X_{0})\times\mathcal{M}_{+}([0,T]\times
X)^{2}\times\mathcal{M}_{+}([0,T]\times X_{u}).$

Corresponding spaces to $\boldsymbol{\ell}$ are

$\displaystyle\mathcal{Y}^{\prime}$
$\displaystyle={\mathcal{C}\times\mathbb{R}}$ (34) $\displaystyle\mathcal{Y}$
$\displaystyle={\mathcal{C}^{\prime}\times 0}.$ (35)

Following the notation of [42], we write that
$\mathcal{Y}_{+}=\\{0_{\mathcal{Y}}\\}$ and
$\mathcal{Y}_{+}^{\prime}=\mathcal{Y}^{\prime}$. We note the containments of
$\boldsymbol{\mu}\in\mathcal{X}_{+},\
\boldsymbol{\ell}\in\mathcal{Y}^{\prime},\ $ and also that
$(\mathcal{X},\mathcal{X}^{\prime})$ form a pair of topological dual spaces
under Assumption A1. The topologies for the spaces $\mathcal{X}$ and
$\mathcal{Y}^{\prime}$ are the weak-* topology and sup-norm-bounded weak
topology respectively.

An adjoint pair of affine maps $\mathcal{A}$ and $\mathcal{A}^{\prime}$ may be
defined as

$\displaystyle\mathcal{A}(\boldsymbol{\mu})$
$\displaystyle=[\mu_{p}+\mu_{c}-\mathcal{L}^{\dagger}\mu-\delta_{t_{0}}\otimes\mu_{0},\
\langle 1,\mu_{0}\rangle]$ (36)
$\displaystyle\mathcal{A}^{*}(\boldsymbol{\ell})$
$\displaystyle=[\gamma-v(t_{0},\bullet),\ -\mathcal{L}_{v},\ v,\ v],$

with vectors describing the cost and constraint terms as

$\displaystyle\mathbf{b}$ $\displaystyle=[0,\ 1]$ (37a)
$\displaystyle\mathbf{c}$ $\displaystyle=[0,\ 0,\ 0,\ 1].$ (37b)

Pairings with the vectors in (37) satisfy

$\displaystyle\langle\mathbf{c},\boldsymbol{\mu}\rangle$
$\displaystyle=\langle 1,\mu\rangle$ (38a)
$\displaystyle\langle\boldsymbol{\ell},\mathbf{b}\rangle$
$\displaystyle=\gamma.$ (38b)

Problem (16) may be expressed as

$\displaystyle p^{*}(t_{0},X_{0})=$
$\displaystyle\sup_{\boldsymbol{\mu}\in\mathcal{X}_{+}}\langle\mathbf{c},\boldsymbol{\mu}\rangle$
$\displaystyle\mathbf{b}-\mathcal{A}(\boldsymbol{\mu})\in\mathcal{Y}_{+}.$
(39) Similarly, the function LP in (18) can be expressed as $\displaystyle
d^{*}(t_{0},X_{0})=$
$\displaystyle\inf_{\boldsymbol{\ell}\in\mathcal{Y}^{\prime}_{+}}\langle\boldsymbol{\ell},\mathbf{b}\rangle$
$\displaystyle\mathcal{A}^{\prime}(\boldsymbol{\ell})-\mathbf{c}\in\mathcal{X}_{+}.$
(40)

Sufficient conditions to prove strong duality between (39) and (40) (by
Theorem 2.6 of [42]) are:

  * R1

All feasible measures $\boldsymbol{\mu}\in\mathcal{X}_{+}$ with
$b-\mathcal{A}(\boldsymbol{\mu})\in\mathcal{Y}_{+}$ are bounded, and there
exists such a feasible $\boldsymbol{\mu}$.

  * R2

All functions used to define $\mathbf{c},\ \mathbf{b},\ \mathcal{A}$ are
continuous.

Boundedness of measures in R1 is proven in Lemma 5.1, and feasibility of a
measure solution is demonstrated by the construction process used in the proof
of Theorem 3.3. For R2, we note that both $\mathbf{c}$ and $\mathbf{b}$ in
(37) are constant and are therefore continuous. Additionally, the set
$\mathcal{C}$ is specifically chosen as the preimage of continuous functions
under $\mathcal{L}$. R1 and R2 are both satisfied under assumptions A1-A6,
thus proving strong duality.

##  Appendix B Superlevel Approximation of Risk Map

This appendix provides the proof of Theorem 4.2.

###  B.1 Supermartingale Property

We first introduce the notion of supermartingales, and then note that any
feasible $v$ for (21b)-(21e) is a supermartingale.

######  Definition B.1.

A process $\\{Y_{t}\\}$ is a supermartingale if $\mathbb{E}[Y_{t+\Delta
t}\mid\\{Y_{t}\\}]\leq Y_{t}$ [29].

######  Proposition B.1.

Any $v$ that satisfies (21b) and (21e) is a supermartingale for the process
$\mathcal{L}$ [28].

###### Proof.

Let $\\{X_{t}\\}$ be the state-dependent distribution for a trajectory of the
stochastic process $\mathcal{L}$ in (1b), respecting exit time
$t^{*}\wedge\tau_{X}$. This set of probability distributions satisfies the
Martingale property of (7) (discrete-time) or (8) (continuous-time) as
appropriate. Given that $\mathcal{L}v\leq 0$ from (21), it holds that
$\mathbb{E}[\mathcal{L}_{0}v(t,x)\mid X_{s^{\prime}}]\leq 0$ for every
stopping time $s^{\prime}$ adapted to $t^{*}\wedge\tau_{X}$. The
supermartingale relation $\mathbb{E}[v(t+s,x)\mid X_{t+s}]\leq v(t,x)$
therefore holds for $v$. ∎

The following lemma uses the previous supermartingale criterion from
Proposition B.1 in combination with Doob’s supermartingale inequality to
provide a probability bound on $v$.

######  Lemma B.2 (based on Lemma 6 of [10]).

Let $v(t,x)\in\mathcal{C}$ be a nonnegative function over $[t_{0},T]\times X$
and also form a supermartingale with respect to $\mathcal{L}$ and $X_{t}$. For
a value $\lambda\geq 0$, an initial point $x_{0}\in X$, and a
$\mathcal{L}$-trajectory $x(t)$ starting from $x_{0}$ at time $t_{0}$, the
following inequality holds:

$\displaystyle\textrm{Prob}_{\\{X_{t}\\}}\left[\sup_{t\in[t_{0},T]}v(t,x(t))\geq\lambda\right]\leq
v(t_{0},x_{0})/\lambda.$ (41)

Lemma B.2 modifies Lemma 6 of [10] by stopping at time $T$ and allowing for
time-dependent functions $v(t,x)$.

###  B.2 Superlevel Property

Theorem 4.2 can now be proven based on arguments from Theorem 7 of [10]. By
constraint (21d), the unsafe set $X_{u}$ is inside the 1-superlevel set of
$v(t,x)$ at every time $t\in[t_{0},T]$. This superlevel relation implies that

$\displaystyle P(t_{0},x_{0})=$
$\displaystyle\sup_{t^{*}\in[t_{0},T]}\textrm{Prob}_{\\{X_{t^{*}}\\}}\left[x(t^{*})\in
X_{u}\mid x_{0}\right]$ (42)
$\displaystyle\leq\textrm{Prob}_{\\{X_{t}\\}}\left[\sup_{t\in[t_{0},T]}v(t,x(t))\geq
1\mid x_{0}\right].$

Through Lemma B.2 with $\lambda=1$, it holds that

$\textrm{Prob}_{\\{X_{t}\\}}\left[\sup_{t\in[t_{0},T]}v(t,x(t))\geq 1\mid
x_{0}\right]\leq v(t_{0},x_{0}),$ (43)

thus proving the theorem.

##  Appendix C Risk Contour Convergence

This appendix will provide a proof of Theorem 4.4.

We first notice that

$\displaystyle 0\leq p^{*}(t_{0},x_{0})$ $\displaystyle=\sup\\{\langle
1,\mu_{p}\rangle\;|\;\mu_{p}+\mu_{c}=\delta_{(t_{0},x_{0})}+\mathcal{L}^{\dagger}\mu\\}$
$\displaystyle=\inf\\{v(t_{0},x_{0})\;|\;\mathcal{L}v\leq 0\;;\;v\geq
I_{X_{u}}\\}$

by applying Theorem 4.2 (using A1-A6) to the case where
$\mu_{0}=\delta_{x_{0}}$ (to have strong duality and thus access the function
interpretation of the problem). We remark that our $v_{k}$ is feasible for the
dual problem on functions (21), so by optimality of $p^{*}(t_{0},x_{0})$ for
each $x_{0}\in X_{0}$, one has $v_{k}(t_{0},x_{0})\geq p^{*}(t_{0},x_{0})$.
Moreover, as $J^{*}(t_{0},\mu_{0})\leq 1<\infty$ by Lemma 4.1, each
$v_{k}(t_{0},\bullet)$ is $\mu_{0}$-integrable, and thus the function
$v^{\star}(t_{0},\bullet)$ is $\mu_{0}$-integrable as well. It remains to
prove that their difference $v^{k}(t_{0},\bullet)-p^{*}(t_{0},\bullet)$
ultimately vanishes in a $\mu_{0}$-sense. This vanishing difference is proven
by contradiction: suppose that there exists a $\eta>0$ such that for all
$k\in\mathbb{N}$, there is a $\varphi(k)\geq k$ with $\langle
v_{\varphi(k)}(t_{0},\bullet)-p^{*}(t_{0},\bullet),\mu_{0}\rangle\geq\eta$
(which is the negation of (23)). Then, one has

$\langle p^{*}(t_{0},\bullet),\mu_{0}\rangle\leq\langle
v_{\varphi(k)}(t_{0},\bullet),\mu_{0}\rangle-\eta\underset{k\to\infty}{\longrightarrow}J^{*}(t_{0},\mu_{0})-\eta.$

However, as we already noticed, $p^{*}(t_{0},x_{0})$ is feasible for (21), so
this is a contradiction with optimality of the value $J^{*}(t_{0},\mu_{0})$
(given that $J^{*}-\eta<J^{*}$). Such a contradiction proves (23), and
therefore proves Theorem 4.4.

## References

  * [1] M. Prandini, J. Hu, C. Cassandras, and J. Lygeros, “Stochastic reachability: Theory and numerical approximation,” _Stochastic hybrid systems, Automation and Control Engineering Series_ , vol. 24, pp. 107–138, 2006.
  * [2] P. Mohajerin Esfahani, D. Chatterjee, and J. Lygeros, “The stochastic reach-avoid problem and set characterization for diffusions,” _Automatica_ , vol. 70, pp. 43–56, 2016.
  * [3] N. Schmid and J. Lygeros, “Probabilistic reachability and invariance computation of stochastic systems using linear programming,” _IFAC-PapersOnLine_ , vol. 56, no. 2, pp. 11 229–11 234, 2023, 22nd IFAC World Congress.
  * [4] S. Asmussen and H. Albrecher, _Ruin Probabilities_. World scientific, 2010, vol. 14.
  * [5] J. Miller, M. Tacchi, M. Sznaier, and A. Jasour, “Peak Value-at-Risk Estimation for Stochastic Processes using Occupation Measures,” 2023, arXiv:2303.16064.
  * [6] M. J. Cho and R. H. Stockbridge, “Linear Programming Formulation for Optimal Stopping Problems,” _SIAM J. Control Optim._ , vol. 40, no. 6, pp. 1965–1982, 2002.
  * [7] A. Abate, M. Prandini, J. Lygeros, and S. Sastry, “Probabilistic reachability and safety for controlled discrete time stochastic hybrid systems,” _Automatica_ , vol. 44, no. 11, pp. 2724–2734, 2008.
  * [8] N. Kariotoglou, M. Kamgarpour, T. H. Summers, and J. Lygeros, “The Linear Programming Approach to Reach-Avoid Problems for Markov Decision Processes,” _Journal of Artificial Intelligence Research_ , vol. 60, pp. 263–285, 2017.
  * [9] V. A. Huynh, S. Karaman, and E. Frazzoli, “An incremental sampling-based algorithm for stochastic optimal control,” in _2012 IEEE International Conference on Robotics and Automation_. IEEE, 2012, pp. 2865–2872.
  * [10] S. Prajna, A. Jadbabaie, and G. J. Pappas, “Stochastic Safety Verification Using Barrier Certificates,” in _2004 43rd IEEE conference on decision and control (CDC)(IEEE Cat. No. 04CH37601)_ , vol. 1. IEEE, 2004, pp. 929–934.
  * [11] A. Clark, “Control barrier functions for stochastic systems,” _Automatica_ , vol. 130, p. 109688, 2021.
  * [12] A. Salamati, A. Lavaei, S. Soudjani, and M. Zamani, “Data-driven safety verification of stochastic systems via barrier certificates,” _IFAC-PapersOnLine_ , vol. 54, no. 5, pp. 7–12, 2021.
  * [13] C. Santoyo, M. Dutreix, and S. Coogan, “A barrier function approach to finite-time stochastic system verification and control,” _Automatica_ , vol. 125, p. 109439, 2021.
  * [14] J. B. Lasserre, _Moments, Positive Polynomials And Their Applications_ , ser. Imperial College Press Optimization Series. World Scientific Publishing Company, 2009.
  * [15] J. Steinhardt and R. Tedrake, “Finite-time Regional Verification of Stochastic Nonlinear System,” _The International Journal of Robotics Research_ , vol. 31, no. 7, pp. 901–923, 2012.
  * [16] B. Xue, N. Zhan, and M. Fränzle, “Reach-avoid analysis for stochastic differential equations,” _arXiv:2208.10752_ , 2022.
  * [17] D. Drzajic, N. Kariotoglou, M. Kamgarpour, and J. Lygeros, “A Semidefinite Programming Approach to Control Synthesis for Stochastic Reach-Avoid Problems.” in _ARCH@ CPSWeek_ , 2016, pp. 134–143.
  * [18] E. Ahbe, A. Iannelli, and R. S. Smith, “Region of attraction analysis of nonlinear stochastic systems using polynomial chaos expansion,” _Automatica_ , vol. 122, p. 109187, 2020.
  * [19] A. Wang, A. Jasour, and B. C. Williams, “Non-Gaussian Chance-Constrained Trajectory Planning for Autonomous Vehicles Under Agent Uncertainty,” _IEEE Robotics and Automation Letters_ , vol. 5, no. 4, pp. 6041–6048, 2020\.
  * [20] D. Henrion, M. Junca, and M. Velasco, “Moment-SOS hierarchy and exit time of stochastic processes,” _arXiv preprint arXiv:2101.06009_ , 2021.
  * [21] X. Chen, S. Chen, and V. M. Preciado, “Safety Verification of Nonlinear Polynomial System via Occupation Measures,” in _2019 IEEE 58th Conference on Decision and Control (CDC)_ , 2019, pp. 1159–1164.
  * [22] J. Miller, D. Henrion, and M. Sznaier, “Peak Estimation Recovery and Safety Analysis,” _IEEE Control Systems Letters_ , vol. 5, no. 6, pp. 1982–1987, 2021.
  * [23] J. Miller and M. Sznaier, “Bounding the Distance to Unsafe Sets with Convex Optimization,” _IEEE Transactions on Automatic Control_ , pp. 1–15, 2023\.
  * [24] ——, “Quantifying the Safety of Trajectories using Peak-Minimizing Control,” 2023, arXiv:2303.11896.
  * [25] C. Sloth and R. Wisniewski, “Safety Analysis of Stochastic Dynamical Systems,” _IFAC-PapersOnLine_ , vol. 48, no. 27, pp. 62–67, 2015, analysis and Design of Hybrid Systems ADHS.
  * [26] A. M. Jasour and B. C. Williams, “Risk Contours Map for Risk Bounded Motion Planning under Perception Uncertainties,” in _Robotics: Science and Systems_ , 2019.
  * [27] D. W. Stroock and S. S. Varadhan, _Multidimensional Diffusion Processes_. Springer Science & Business Media, 1997, vol. 233.
  * [28] L. C. Rogers and D. Williams, _Diffusions, Markov Processes, and Martingales: Volume 1, foundations_. Cambridge university press, 2000, vol. 1.
  * [29] D. Williams, _Probability with Martingales_. Cambridge university press, 1991.
  * [30] D. Hilbert, “Über die Darstellung definiter Formen als Summe von Formenquadraten,” _Mathematische Annalen_ , vol. 32, no. 3, pp. 342–350, 1888.
  * [31] G. Blekherman, “There are Significantly More Nonnegative Polynomials than Sums of Squares,” _Israel Journal of Mathematics_ , vol. 153, no. 1, pp. 355–380, 2006.
  * [32] M.-D. Choi, T. Y. Lam, and B. Reznick, “Sums of squares of real polynomials,” in _Proceedings of Symposia in Pure mathematics_ , vol. 58. American Mathematical Society, 1995, pp. 103–126.
  * [33] J. Cimprič, M. Marshall, and T. Netzer, “Closures of quadratic modules,” _Israel Journal of Mathematics_ , vol. 183, no. 1, pp. 445–474, 2011.
  * [34] M. Putinar, “Positive Polynomials on Compact Semi-algebraic Sets,” _Indiana University Mathematics Journal_ , vol. 42, no. 3, pp. 969–984, 1993\.
  * [35] J. Nie and M. Schweighofer, “On the complexity of Putinar’s Positivstellensatz,” _Journal of Complexity_ , vol. 23, no. 1, pp. 135–150, 2007.
  * [36] D. Henrion, J. B. Lasserre, and C. Savorgnan, “Approximate Volume and Integration for Basic Semialgebraic Sets,” _SIAM review_ , vol. 51, no. 4, pp. 722–743, 2009.
  * [37] S. Prajna, A. Jadbabaie, and G. J. Pappas, “A Framework for Worst-Case and Stochastic Safety Verification Using Barrier Certificates,” _IEEE Transactions on Automatic Control_ , vol. 52, no. 8, pp. 1415–1428, 2007.
  * [38] M. Tacchi, “Convergence of Lasserre’s hierarchy: the general case,” _Optimization Letters_ , vol. 16, no. 3, pp. 1015–1033, 2022.
  * [39] J. Miller, T. Dai, and M. Sznaier, “Data-Driven Superstabilizing Control of Error-in-Variables Discrete-Time Linear Systems,” in _2022 61st IEEE Conference on Decision and Control (CDC)_ , 2022, pp. 4924–4929.
  * [40] J. Lofberg, “YALMIP : a toolbox for modeling and optimization in MATLAB,” in _ICRA (IEEE Cat. No.04CH37508)_ , 2004, pp. 284–289.
  * [41] M. ApS, _The MOSEK optimization toolbox for MATLAB manual. Version 10.1._ , 2023\.
  * [42] M. Tacchi, “Moment-sos hierarchy for large scale set approximation. application to power systems transient stability analysis,” Ph.D. dissertation, Toulouse, INSA, 2021.

  *[SDE]: Stochastic Differential Equation
  *[LP]: Linear Program
  *[SOS]: Sum of Squares
  *[LPs]: Linear Program
  *[SDPs]: Semidefinite Program
  *[SDEs]: Stochastic Differential Equation
  *[PSD]: Positive Semidefinite
  *[BSA]: Basic Semialgebraic
  *[WSOS]: Weighted Sum of Squares


# Temperature dependence of divacancy spin coherence in implanted silicon
carbide

Wu-Xi Lin  CAS Key Laboratory of Quantum Information, University of Science
and Technology of China, Hefei 230026, People’s Republic of China  CAS Center
for Excellence in Quantum Information and Quantum Physics, University of
Science and Technology of China, Hefei 230026, People’s Republic of China
Fei-Fei Yan  CAS Key Laboratory of Quantum Information, University of Science
and Technology of China, Hefei 230026, People’s Republic of China  CAS Center
for Excellence in Quantum Information and Quantum Physics, University of
Science and Technology of China, Hefei 230026, People’s Republic of China
Qiang Li<EMAIL_ADDRESS>CAS Key Laboratory of Quantum Information,
University of Science and Technology of China, Hefei 230026, People’s Republic
of China  CAS Center for Excellence in Quantum Information and Quantum
Physics, University of Science and Technology of China, Hefei 230026, People’s
Republic of China  Jun-feng Wang  CAS Key Laboratory of Quantum Information,
University of Science and Technology of China, Hefei 230026, People’s Republic
of China  CAS Center for Excellence in Quantum Information and Quantum
Physics, University of Science and Technology of China, Hefei 230026, People’s
Republic of China  Zhi-He Hao  CAS Key Laboratory of Quantum Information,
University of Science and Technology of China, Hefei 230026, People’s Republic
of China  CAS Center for Excellence in Quantum Information and Quantum
Physics, University of Science and Technology of China, Hefei 230026, People’s
Republic of China  Ji-Yang Zhou  CAS Key Laboratory of Quantum Information,
University of Science and Technology of China, Hefei 230026, People’s Republic
of China  CAS Center for Excellence in Quantum Information and Quantum
Physics, University of Science and Technology of China, Hefei 230026, People’s
Republic of China  Hao Li  Li-Xing You  State Key Laboratory of Functional
Materials for Informatics, Shanghai Institute of Microsystem and Information
Technology, Chinese Academy of Sciences(CAS), Shanghai 200050, People’s
Republic of China  Jin-Shi Xu<EMAIL_ADDRESS>CAS Key Laboratory of Quantum
Information, University of Science and Technology of China, Hefei 230026,
People’s Republic of China  CAS Center for Excellence in Quantum Information
and Quantum Physics, University of Science and Technology of China, Hefei
230026, People’s Republic of China  Chuan-Feng Li<EMAIL_ADDRESS>CAS Key
Laboratory of Quantum Information, University of Science and Technology of
China, Hefei 230026, People’s Republic of China  CAS Center for Excellence in
Quantum Information and Quantum Physics, University of Science and Technology
of China, Hefei 230026, People’s Republic of China  Guang-Can Guo  CAS Key
Laboratory of Quantum Information, University of Science and Technology of
China, Hefei 230026, People’s Republic of China  CAS Center for Excellence in
Quantum Information and Quantum Physics, University of Science and Technology
of China, Hefei 230026, People’s Republic of China

###### Abstract

Spin defects in silicon carbide (SiC) have attracted increasing interest due
to their excellent optical and spin properties, which are useful in quantum
information processing. In this paper, we systematically investigate the
temperature dependence of the spin properties of divacancy defects in
implanted 4 _H_ -SiC. The zero-field splitting parameter $D$, the
inhomogeneous dephasing time $T_{2}^{*}$, the coherence time $T_{2}$, and the
depolarization time $T_{1}$ are extensively explored in a temperature range
from 5 to 300 K. Two samples implanted with different nitrogen molecule ion
fluences ($\rm{N_{2}}^{+}$, $1\times 10^{14}/\rm cm^{2}$ and $1\times
10^{13}/\rm cm^{2}$) are investigated, whose spin properties are shown to have
similar temperature-dependent behaviors. Still, the sample implanted with a
lower ion fluence has longer $T_{2}$ and $T_{1}$. We provide possible
theoretical explanations for the observed temperature-dependent dynamics. Our
work promotes the understanding of the temperature dependence of spin
properties in solid-state systems, which can be helpful for constructing wide
temperature-range thermometers based on the mature semiconductor material.

##  I introduction

Spin defects in solid-state materials have become an important candidate for
qubits during the last two decades for their possible applications in quantum
communication, quantum computation and quantum sensing (AwschalomSum, ;
AwschalomCompute, ; WrachtrupDiamond, ; AwschalomRoom, ; WrachtrupSingle, ;
Dimond2, ; Dimond3, ; AwschalomPolytype, ; Astakhov, ; Castelletto, ; SiCAdd1,
; SiCAdd2, ; Lienhard, ; SiCAdd3, ; SiCAddLab, ). Nitrogen-vacancy (NV)
centers in diamond have been most frequently explored (WrachtrupDiamond, ;
Dimond2, ; Dimond3, ). However, because the fabrication of suitable diamond-
based devices for large-scale practical applications is still met with great
challenges, much effort has been put into the search for similar spin defects
in more mature solid-state materials. Silicon carbide is an outstanding
semiconductor material that has been widely used in many electronic devices
(Material1, ; Material2, ) due to its excellent properties such as high
thermal conductivity, chemical stability, and controllable
electroconductivity. There are sophisticated technologies for the large-scale
crystal growth of SiC and mature processes for the industrial production of
SiC-based devices. So, the defects in SiC, which have been explored as bright
single-photon sources (Astakhov, ; Castelletto, ; Lienhard, ) and potential
qubits (AwschalomRoom, ; WrachtrupSingle, ; AwschalomPolytype, ; SiCAdd1, ;
SiCAdd2, ; SiCAdd3, ; SiCAddLab, ), have attracted great interest.

SiC has many polytypes and one of the most popular polytypes is 4 _H_ -SiC,
which has a high crystal quality and a layered hexagonal crystal structure
(AwschalomRoom, ). One of the most frequently studied defects in 4 _H_ -SiC is
the divacancy defect consisting of a silicon vacancy and an adjacent carbon
vacancy, which has spin $S=1$ (AwschalomRoom, ; 4HSiC, ). There are two kinds
of divacancies in 4 _H_ -SiC, the _c_ -axis ones with $C_{3v}$ symmetry and
the _basal_ -type ones with $C_{1h}$ symmetry (AwschalomPolytype, ). For a _c_
-axis divacancy, we mark its axis of symmetry as the _c_ axis and the
direction of its spin is parallel to the _c_ axis. Recently, divacancies
inside stacking faults have been shown to have excellent optical and spin
properties at room temperature. Especially, the _c_ -axis divacancy named PL6
is shown to have bright photoluminescence fluorescence and a high spin readout
contrast (Li Qiang's paper, ), which are comparable to those of NV centers in
diamond.

In this paper, we systematically investigate the temperature dependence of the
spin properties of the PL6 defects from 5 to 300 K in 4 _H_ -SiC samples
implanted by nitrogen molecule ion ($\rm{N_{2}}^{+}$) beams with two different
ion fluences (the two samples are marked as sample A and sample B, where
sample A is implanted with an ion fluence of $1\times 10^{14}/\rm cm^{2}$, and
sample B is implanted with an ion fluence of $1\times 10^{13}/\rm cm^{2}$).
The temperature-dependent behaviors of four important spin properties, the
zero-field splitting (ZFS) parameter $D$, the inhomogeneous dephasing time
$T_{2}^{*}$, the coherence time $T_{2}$, and the depolarization time $T_{1}$,
are measured. $D$ decreases monotonically as the temperature increases. The
inhomogeneous dephasing time $T_{2}^{*}$ fluctuates in a small time region at
a wide range of temperatures. Although the coherence time $T_{2}$ decreases as
the temperature increases, it is revived in a temperature range between 200
and 300 K. The inverse of the depolarization time $1/T_{1}$ shows an almost
linear temperature-dependent property under 200 K and shows an approximately
polynomial temperature dependence above 250 K. The temperature-dependent spin
dynamics are similar for samples A and B implanted with these two different
fluences of nitrogen molecule ions, but at the same temperature, the low-dose
implanted sample B with a smaller PL6 concentration has longer $T_{2}$ and
$T_{1}$. We discuss possible theoretical explanations for the observed
temperature-dependent dynamics. This paper provides a profound understanding
of the temperature-related behaviors of the divacancies in 4 _H_ -SiC, which
can be used for making wide temperature-range thermometers and designing
quantum devices based on silicon carbide.

##  II Experimental Setup and Results

The experimental setup is shown in Fig. 1. An infrared laser with a central
wavelength of 920 nm modulated by an acousto-optic modulator (AOM) is used to
pump the defect spins in the sample. The laser is focused on the sample by an
near-infrared objective [Olympus, numerical aperture (NA) = 0.65] after being
reflected by a 980-nm long-pass dichroic mirror (DM, Thorlabs). The
fluorescence is collected by the same objective and further filtered by a
1000-nm long-pass interference filter (Thorlabs), then detected by a
superconducting single-photon detector (Photon Technology). The sample is
bonded on a metal circuit board with two electrodes, and we weld a
20-$\rm\upmu m$-diameter copper wire between the electrodes. The wire is
attached to the surface of the sample. A microwave source (Mini-circuits,
SSG6000) modulated by a microwave switch (Mini-circuits, ZASWA-2-50DR+) and
amplified by a microwave amplifier (Mini-circuits, ZHL-30W-252+) is used to
generate the microwave pulses that can manipulate the spin states of PL6
defects. The pulse sequences of the laser and microwave are controlled by a
computer. The experimental results in this paper are obtained from the change
of photoluminescence fluorescence excited by the laser between the two
circumstances which are with or without the action of the microwave, which is
marked as $\Delta\rm PL$ in arbitrary units (a.u.).

Samples are sliced from a commercially available high-purity 4 _H_ -SiC
epitaxial wafer with a 7-$\rm\upmu m$-thick epitaxial layer. According to some
previous works, a hydrogen ion ($\rm H^{+}$) beam (SiCAddLab, ) or carbon ion
($\rm C^{+}$) beam (Li Qiang's paper, ) can be used to generate the PL6
defects in 4 _H_ -SiC. However, nitrogen molecule ion ($\rm{N_{2}}^{+}$)
implantation is rarely used to generate these kinds of defects in 4 _H_ -SiC.
But since $\rm{N_{2}}^{+}$ has a larger collision cross section than those of
$\rm H^{+}$ and $\rm C^{+}$ when hitting the sample, we can reasonably guess
that the $\rm{N_{2}}^{+}$ ion fluence can also effectively generate PL6. So,
in this work, we use $\rm{N_{2}}^{+}$ ions with a kinetic energy of 30 keV to
implant samples A and B with fluences of $1\times 10^{14}/\rm cm^{2}$ and
$1\times 10^{13}/\rm cm^{2}$, respectively. After implantation, the samples
are annealed at 1050∘C for 30 min. Under the above implantation and annealing
conditions, the PL6 defects in 4 _H_ -SiC can be efficiently generated. The
PL6 concentration in sample B implanted with an $\rm{N_{2}}^{+}$ ion fluence
of $1\times 10^{13}/\rm cm^{2}$ is smaller and this low-dose implanted sample
is used for comparative experiments.

The sample is mounted on the stage in a cryogenic vacuum chamber (Montana
Instruments) with the temperature at around 4 K. The temperature is controlled
by a heater whose temperature controlling range is between 5 and 300 K with an
accuracy of 0.001 K. In order to split the approximately degenerate spin
states, we put the sample in a static magnetic field, which is provided by a
permanent magnet placed in this chamber. In the experiment, the magnetic field
is set to be $1.80\times 10^{-2}$ T.

Figure 1: Experimental setup. The sample is mounted on a metal circuit board
(MCB), which is attached to the stage in a cryogenic vacuum chamber (CVC). The
temperature is controlled by a heater controller (HC). The pump laser
modulated by an acousto-optic modulator (AOM) is focused by an objective to
the sample. The fluorescence is collected by the same objective and filtered
by a dichroic mirror (DM) and an interference filter (IF), which is then
coupled by a fiber coupler (FC) and sent through the fiber to a
superconducting single-photon detector (SSPD). The pump laser and microwave
(MW) pulse sequences are controlled by a computer. The magnetic field is
provided by a permanent magnet in the cryogenic chamber.

We first investigate the spin properties of sample A implanted by nitrogen
molecule ions with a fluence of $1\times 10^{14}/\rm cm^{2}$. The ground-state
spin Hamiltonian of PL6 defects can be written as (WrachtrupSingle, )

$H=D[S_{z}^{2}-\frac{1}{3}S(S+1)]+E(S_{x}^{2}-S_{y}^{2})+g\mu_{B}\vec{S}\cdot\vec{B},$
(1)

where $g\approx 2$ is the electron’s Lande-$g$ factor and $\mu_{B}$ is the
Bohr magneton. $D$ and $E$ are both parameters describing the zero-field
splitting which is caused by the electric dipole interaction between
electrons, and they are related to the axisymmetric and non-axisymmetric part
of the electric dipole interaction, respectively. For PL6, $D\gg E$
(AwschalomPolytype, ), and we only focus on $D$. The $S_{x}$, $S_{y}$, and
$S_{z}$ are three components of the spin ($S=1$). There are two possible
transitions between states (denoted as the magnetic quantum numbers)
$\mathinner{|m_{s}=0\rangle}\leftrightarrow\mathinner{|m_{s}=+1\rangle}$ and
$\mathinner{|m_{s}=0\rangle}\leftrightarrow\mathinner{|m_{s}=-1\rangle}$,
respectively. As the applied _c_ -axis static magnetic field $\vec{B}$ is
parallel to $\vec{S}$, the corresponding energy splittings which can be
derived from the Hamiltonian are
$\omega_{1}=D-\sqrt{(g\mu_{B}B/\hbar)^{2}+E^{2}}$ and
$\omega_{2}=D+\sqrt{(g\mu_{B}B/\hbar)^{2}+E^{2}}$, so the ZFS parameter $D$
can be obtained by $D=\frac{1}{2}(\omega_{1}+\omega_{2})$.

Employing optically detected magnetic resonance (ODMR) technology (ODMR, ),
the resonant frequencies can be obtained. Figure 2(a) shows three
representative ODMR spectra at three different temperatures, 5, 140, and 300
K. The values of $\omega_{1}$ and $\omega_{2}$ are obtained from the
Lorentzian fitting lines in the left and right panels, respectively. We will
focus on the resonant frequency $\omega_{1}$ in all the coherent control
experiments described below. Figure 2(b) shows the relationship between the
ZFS parameter $D$ and the temperature. The change of $D$ is about $-13$ MHz
from 5 to 300 K. There are several different formulas that can be used to fit
the curve shown in Fig. 2(b): the Debye-model formula
$D=[1304.1+60.6\exp(-2.7\times 10^{-6}K^{-2}T^{2})]\ \rm MHz$, which is
derived from the simple Debye model theoretically (D1, ); the Varshni-form
formula $D=(1364.7-0.2K^{-1}\frac{T^{2}}{1348.2K+T})\ \rm MHz$ (D2, ); and the
polynomial-form formula $D=(1364.6+3.5\times 10^{-3}K^{-1}T-1.8\times
10^{-4}K^{-2}T^{2}-1.5\times 10^{-7}K^{-3}T^{3}+1.6\times
10^{-9}K^{-4}T^{4}-2.7\times 10^{-12}K^{-5}T^{5})\ \rm MHz$ (D3, ; Important2,
). For all the above fitting functions, the determination coefficients satisfy
$R^{2}>0.99$, so all the functions are suitable for fitting the temperature
dependence of $D$. Because of the comparatively more solid theoretical
foundation (D1, ), we prefer to use the Debye-model formula to fit the
$D$-temperature curve, which is shown in Fig. 2(b). The temperature dependence
of $D$ can be used in the high-sensitivity quantum sensing of temperature
(Important2, ; D2, ; D3, ; Thermal1, ; Thermal2, ).

Because the temperature dependence of the spin properties of the defects can
be affected by their distribution, we also use the stopping and range of ions
in matter (SRIM) (SRIM_1, ; SRIM_2, ) to simulate the depth distribution of
the vacancy defects (including PL6) generated by the 30-keV $\rm{N_{2}}^{+}$
implanted ions. Figure 2(c) shows the simulation result of the normalized
concentration of vacancy defects as a function of the distance to the sample
surface. Deduced from the result, most of the vacancy defects, including the
PL6 defects, are distributed within 50 nm from the surface of the sample.

Figure 2: Temperature dependence of the ODMR spectrum and the ZFS parameter
$D$ of sample A implanted with a $1\times 10^{14}/\rm cm^{2}$ nitrogen
molecule ion fluence along with the simulation result of the vacancy defect’s
depth distribution by using SRIM. In (a) and (b) the dots are experimental
data and the solid lines are fitting lines, while in (c) the solid line is the
simulation curve. (a) The ODMR spectra of PL6 at 5, 140, and 300 K in a
$1.80\times 10^{-2}$ T static magnetic field. The fitting lines are Lorentzian
lines. $\omega_{1}$ and $\omega_{2}$ are the two corresponding resonant
frequencies, respectively, and $D=(\omega_{1}+\omega_{2})/2$. (b) Temperature
dependence of $D$. The fitting line is the Debye-model fitting line with a
determination coefficient of $0.9997$. $D$ decreases monotonically as the
temperature increases. The error bars are deduced from the fitting errors. (c)
The simulation result that shows the normalized vacancy concentration
generated by the $\rm{N_{2}}^{+}$ ion implantation as a function of the depth
from the sample surface. This result demonstrates that most vacancy defects,
including PL6 defects, are distributed within 50 nm from the sample surface.

We then measure the inhomogeneous dephasing time $T_{2}^{*}$ of PL6 at
different temperatures inferred from the Ramsey interference curve.
$T_{2}^{*}$ is an important factor that affects the performance of quantum
sensors based on Ramsey interference (Important2, ; Thermal1, ). The microwave
pulse sequence used in the Ramsey interference experiments is
$\pi/2-\tau-\pi/2$. $\pi/2$ means the duration time of the microwave is
$\pi/(2\Omega_{R})$, where $\Omega_{R}$ is the angular frequency of the Rabi
oscillation (Rabi, ) of the spin system under the action of the microwave with
the same power. $\tau$ is the length of the time interval during which the
system is allowed to evolve freely without the action of the microwave. We set
the detuning frequency of the microwave as $\delta=10\ \rm MHz$. Figure
3(a)-3(c) show three representative Ramsey interference curves at 50, 150, and
300 K, respectively. We fit these curves by the formula
$a\exp[-(\tau/T_{2}^{*})^{2}]\cos(2\pi\delta\tau+\phi)+b$, where $a$, $b$, and
$\phi$ are free parameters, and we take the exponent of $\tau/T_{2}^{*}$ as 2
because of a reasonable hypothesis that the noise which causes the signal to
decay has a Gaussian distribution, as shown in the former theoretical results
(Ramsey, ; Ramsey2, ). $T_{2}^{*}$ can then be derived from the fitting
formula. Figure 3(d) shows the temperature dependence of $T_{2}^{*}$, from
which we can find that $T_{2}^{*}$ remains almost the same. The collapse of
the Ramsey interference oscillation is mainly due to the interaction between
the defect spins and the nuclear spins, the fluctuation of the external
magnetic field, and the instability of the frequency of the microwave. The
increase in temperature mainly influences the intensity of the lattice
vibration, which does not much affect the behavior of inhomogeneous dephasing.

Figure 3: Temperature dependence of Ramsey interference and the inhomogeneous
dephasing time $T_{2}^{*}$ of sample A implanted with a $1\times 10^{14}/\rm
cm^{2}$ nitrogen molecule ion fluence. (a)-(c) Ramsey interference curves at
50, 150, and 300 K with $T_{2}^{*}$ derived from the fitting formula,
respectively. The dots are experimental data and the solid lines are the
fitting lines. (d) Temperature dependence of $T_{2}^{*}$ from 5 to 300 K.
$T_{2}^{*}$ remains stable at different temperatures. Error bars are deduced
from the fitting errors.

Since the coherence time $T_{2}$ is critical in the spin-based quantum
technology (Divincenzo, ), we further measure its temperature-dependent
behavior. The standard method of measuring $T_{2}$ is to measure the spin echo
oscillation, and the microwave pulse sequence used in this experiment is
$\pi/2-\tau/2-\pi-\tau/2-\pi/2$. Figure 4(a) shows a representative spin echo
curve at 5 K. Because of the existence of the static magnetic field, the spin
echo attenuated oscillation is an electron spin echo envelope modulation
(ESEEM) curve (Nuclear1, ; Nuclear2, ; Nuclear3, ). The experimental data are
fitted by the formula $a\exp(-\tau/T_{2})[1-b\sin^{2}(\pi
f_{1}\tau)][1-c\sin^{2}(\pi f_{2}\tau)]+d$, where $a$, $b$, $c$, and $d$ are
free parameters. $f_{1}$ and $f_{2}$ are Larmor precession frequencies of the
spin magnetic moments of the 13C and 29Si nuclei in the static magnetic field.
Figure 4(b) shows the temperature dependence of $T_{2}$. We can find an
unusual phenomenon, that is, in a specific temperature range between 200 and
300 K, $T_{2}$ becomes larger as the temperature increases. But usually the
coherence time of the defects in solid-state systems decreases monotonically
as the temperature increases. Similar behavior was observed from silicon
vacancy spins in 4 _H_ -SiC (Important1, ).

In order to explain this unusual nonmonotonic relationship between $T_{2}$ and
temperature, we assume that it may be attributed to the dynamic Jahn-Teller
effect (Important1, ). In the 4 _H_ -SiC crystal, there are carbon and silicon
atoms adjacent to the divacancy. At low temperatures, orbital pairing occurs
between neighboring atoms’ electrons and leads to the Jahn-Teller distortion,
which defines the structure of the divacancy. As the temperature increases to
a specific value, the energy of the thermal vibration is large enough for
unraveling the orbital pairs, so rapid thermally activated reorientation
occurs and the Jahn-Teller distortion disappears. The restructuring of the
divacancy induced by the orbital pairs’ unraveling can possibly lead the
coherence time to increase and we can find the abnormal phenomenon that
$T_{2}$ is longer at a higher temperature than at a lower temperature.

Figure 4: Temperature dependence of the coherence time $T_{2}$ and the
depolarization time $T_{1}$ of sample A implanted with a $1\times 10^{14}/\rm
cm^{2}$ nitrogen molecule ion fluence. (a) The spin echo oscillation at 5 K
with $T_{2}=30.7\ \rm\upmu s$ derived from the fitting formula (the solid
line). (b) Temperature dependence of $T_{2}$ from 5 to 300 K. An anomalous
revival of $T_{2}$ occurs at around 250 K. (c) The depolarization curve at 5 K
with $T_{1}=567\ \rm\upmu s$ derived from the fitting formula (the solid
line). (d) Temperature dependence of $1/T_{1}$ from 5 to 300 K, which transits
from a linear relationship to a polynomial relationship. The dots are
experimental data. Error bars are deduced from the fitting errors.

The coherence time $T_{2}$ of defects is restricted by the depolarization time
$T_{1}$, which is given by $T_{2}<0.5T_{1}$ (T2T11, ) or $T_{2}<1T_{1}$
(T2T12, ). We measure the temperature-dependent behavior of $T_{1}$. The spin
state is first polarized and then evolves freely for a period of time $\tau$,
from which the photoluminescence fluorescence is detected. A similar detection
is then performed in the case with a $\pi$ microwave pulse after the
polarizing process. The fluorescence difference between these two processes is
measured as the function of $\tau$. Figure 4(c) shows a representative
depolarization curve at 5 K. The experimental data are fitted by the
exponential decay formula $a\exp(-\tau/T_{1})+b$, where $a$ and $b$ are free
parameters. $T_{1}$ can then be derived from the fitting formula. Figure 4(d)
shows the temperature dependence of $1/T_{1}$. In the temperature region from
5 to 200 K, $1/T_{1}$ increases almost linearly as the temperature increases.
Two main factors are involved to cause the spin depolarization, including the
direct energy transition between the defect spins and the phonons, and the
Orbach mechanism. $1/T_{1}$ can be described by $1/T_{1}\approx
a\coth(E_{1}/k_{B}T)+b/[\exp(E_{2}/k_{B}T)-1]+c$ (T1theo1, ; T1theo2, ), where
$k_{B}$ is the Boltzmann’s constant and $c$ is a parameter independent of
temperature. $E_{1}$ is proportional to the energy splitting between the spin
states, and it satisfies $k_{B}T\gg E_{1}$ above 1 K. $E_{2}$ is proportional
to the energy splitting between the orbital ground state and a nearby orbital
excited state, which always satisfies $k_{B}T\ll E_{2}$ below 300 K. So we
have $\coth(E_{1}/k_{B}T)\approx k_{B}T/E_{1}$ and
$b/[\exp(E_{2}/k_{B}T)-1]\approx 0$ and the relationship becomes
$1/T_{1}\approx a^{\prime}T+c^{\prime}$. When the temperature is above 250 K,
the depolarization of the spin state is mainly caused by the two-phonon Raman
process, which can be described by (T1theo1, ; T1theo2, ) $1/T_{1}\approx
aT^{s}+bT^{s+1}+cT^{s+2}$. $s=2d-1$ and $d$ is the dimension of the sample
with the main term being $aT^{s}$. Because PL6 defects are mainly distributed
near the sample surface [as the simulation result displayed in Fig. 2(c)
shows], we can approximately take $d=2$ and have $1/T_{1}\propto T^{3}$. From
the experimental result we find that
$[T_{1}(T=250K)/T_{1}(T=300K)]^{1/3}\approx 1.9^{1/3}\approx 1.24\approx
300K/250K$, which is approximately consistent with the theoretical analysis.
The linear-form temperature dependence of $1/T_{1}$ transits to the
polynomial-form temperature dependence in the temperature range between 200
and 250 K. For this sample A implanted with the $1\times 10^{14}/\rm cm^{2}$
ion fluence, the longest $T_{1}$ is measured to be about 500 $\rm\upmu s$,
while the longest $T_{2}$ is about 30 $\rm\upmu s$. The coherence time can be
increased further by using the dynamical decoupling method (Science 330, ) and
increasing the quality of the crystals, for example, reducing the defects in
the samples and using the isotopic purification samples (NatureMat 19, ;
NatureCom 11, ).

Figure 5: Experimental results of sample B implanted with a nitrogen molecule
ion fluence of $1\times 10^{13}/\rm cm^{2}$. (a) Typical ODMR spectra at 5,
140, and 300 K. The ODMR signals of the low-dose implanted sample B have
smaller intensities, smaller signal-to-noise ratios, and narrower full width
at half maxima (FWHM) than those of the high-dose implanted sample A. (b)
Temperature dependence of $D$. The solid line is the Debye-model fitting curve
with the determination coefficient $R^{2}=0.9986$. (c) Temperature dependence
of $T_{2}^{*}$. $T_{2}^{*}$ remains stable at different temperatures. (d)
Temperature dependence of $T_{2}$. An anomalous revival of $T_{2}$ occurs at
around 250 K. (e) Temperature dependence of $1/T_{1}$, which transits from a
linear relationship to a polynomial relationship. At the same temperature,
sample B has longer $T_{2}$ and $T_{1}$ than those of sample A implanted with
an ion fluence of $1\times 10^{14}/\rm cm^{2}$. Error bars are deduced from
the fitting errors.

Finally, we provide the experimental results of PL6 defects in sample B
implanted with a $1\times 10^{13}/\rm cm^{2}$ nitrogen molecule ion fluence.
Figure 5(a) shows typical ODMR spectra of this low-dose implanted sample B at
5, 140, and 300 K. The ODMR spectra in Figs. 2(a) and 5(a) are normalized by
the same normalization factor. From the two sets of ODMR spectra we can find
that the ODMR signal intensity of sample B is about one-third to one-half of
the ODMR signal intensity of sample A at the same temperature, which shows
that a lower dose of implanted ions can indeed produce a smaller concentration
of PL6 defects. We can also find that the low-dose implanted sample B’s ODMR
signal has a lower signal-to-noise ratio and a narrower full width at half
maximum (FWHM) than those of the high-dose implanted sample A’s ODMR signal at
the same temperature. Figure 5(b) shows the temperature dependence of the ZFS
parameter $D$, which can also be well fitted by the Debye-model formula
$D=[1301.3+64.3\exp(-2.5\times 10^{-6}K^{-2}T^{2})]\ \rm MHz$ (D1, ), Varshni-
form formula $D=(1365.6-0.2K^{-1}\frac{T^{2}}{1268.3K+T})\ \rm MHz$ (D2, ),
and polynomial-form formula $D=(1365.4+1.2\times 10^{-2}K^{-1}T-3.5\times
10^{-4}K^{-2}T^{2}+8.8\times 10^{-7}K^{-3}T^{3}-4.2\times
10^{-10}K^{-4}T^{4}-2.2\times 10^{-12}K^{-5}T^{5})\ \rm MHz$ (D3, ;
Important2, ), respectively. The fitting curve that we show in Fig. 5(b) is
the Debye-model curve. $D$ also changes about $-13$ MHz from 5 to 300 K.
Figure 5(c) shows the temperature dependence of $T_{2}^{*}$, from which we can
find that $T_{2}^{*}$ still remains stable for this low-dose implanted sample
B. The value of $T_{2}^{*}$ is almost not affected by the implanted ion
fluence. An abnormal nonmonotonic behavior is still observed in the
temperature dependence of $T_{2}$ in the temperature range between 200 and 300
K, which is shown in Fig. 5(d). The values of $T_{2}$ are almost two times
those shown in Fig. 4(b), which may due to the smaller total dipole
interaction between the defect spins and environmental nuclear spins in this
low-dose implanted sample B with a smaller defect concentration. Figure 5(e)
shows the temperature dependence of $T_{1}$. The behavior is similar to that
shown in Fig. 4(d), but this low-dose implanted sample B has a longer
depolarization time. The sample with a higher defect concentration exists with
larger resonant magnetic or electric noises, which depolarize the spin states
(T1concen1, ; T1concen2, ; T1concen3, ).

##  III Conclusion

In this paper, we explore the temperature dependence of PL6 defect spin
properties in 4 _H_ -SiC implanted with two different nitrogen molecule ion
fluences between 5 and 300 K. The ZFS parameter is shown to decrease
monotonically as the temperature increases. It changes about $-13$ MHz from 5
to 300 K, which can be used for constructing wide temperature-range
thermometers. The inhomogeneous dephasing time is shown to be stable and
remains almost the same in the two samples implanted with different ion
fluences. Although the evolution trends are similar for both samples, the
coherence time and depolarization time are longer in the sample implanted with
a lower ion fluence. An increase in the coherence time is observed during the
increase of temperature in both samples, which is attributed to the
disappearance of the Jahn-Teller distortion. We discuss possible theoretical
explanations for the observed phenomena. Our work is useful in the
understanding of temperature-related spin dynamics in SiC, and will stimulate
further investigations in this area.

## ACKNOWLEDGMENTS

This work was supported by the National Key Research and Development Program
of China (Grants No. 2016YFA0302700), the National Natural Science Foundation
of China (Grants No. U19A2075, 61725504, 61905233, 11774335, 11821404 and
11975221), the Key Research Program of Frontier Sciences, Chinese Academy of
Sciences (CAS) (Grant No. QYZDY-SSW-SLH003), the Anhui Initiative in Quantum
Information Technologies (Grants No. AHY060300 and No. AHY020100), and the
Fundamental Research Funds for the Central Universities (Grants No.
WK2030380017 and WK2470000026), the National Postdoctoral Program for
Innovative Talents (Grant No. BX20200326). This work was partially performed
at the University of Science and Technology of China Center for Micro and
Nanoscale Research and Fabrication.

W.-X. L. and F.-F. Y. contributed equally to this work.

## References

  * (1) D. D. Awschalom, L.C. Bassett, A. S. Dzurak, E. L. Hu, and J. R. Petta,  Quantum Spintronics: Engineering and Manipulating Atom-Like Spins in Semiconductors,  Science 339, 1174 (2013).
  * (2) J. R. Weber, W. F. Koehl, J. B. Varley, A. Janotti, B. B. Buckley, C. G. Van de Walle, and D. D. Awschalom,  Quantum computing with defects,  MRS bulletin 38, 802 (2013).
  * (3) F. Jelezko and J. Wrachtrup,  Single defect centres in diamond: A review,  Physica Status Solidi 203, 3207 (2010).
  * (4) R. Schirhagl, K. Chang, M. Loretz, and C. L. Degen,  Nitrogen-Vacancy Centers in Diamond: Nanoscale Sensors for Physics and Biology,  Annual Review of Physical Chemistry 65, 83 (2014).
  * (5) T. Gaebel, M. Domhan, I. Popa, C. Wittmann, P. Neumann, F. Jelezko, J. R. Rabeau, N. Stavrias, A. D. Greentree, and S. Prawer,  Room-temperature coherent coupling of single spins in diamond,  Nat. Phys. 2, 408 (2006).
  * (6) F. Fuchs, B. Stender, M. Trupke, D. Simin, J. Pflaum, V. Dyakonov, and G. Astakhov,  Engineering near-infrared single-photon emitters with optically active spins in ultrapure silicon carbide,  Nat. Commun. 6, 7578 (2015).
  * (7) S. Castelletto, B. Johnson, V. Ivády, N. Stavrias, T. Umeda, A. Gali, and T. Ohshima,  A silicon carbide room-temperature single-photon source,  Nat. Mater. 13, 151 (2014).
  * (8) B. Lienhard, T. Schröder, S. Mouradian, F. Dolde, T. T. Tran, I. Aharonovich, and D. Englund,  Bright and photostable single-photon emitter in silicon carbide,  Optica 3, 768 (2016).
  * (9) W. F. Koehl, B. B. Buckley, F.J. Heremans, G. Calusine and D. D. Awschalom,  Room temperature coherent control of defect spin qubits in silicon carbide,  Nature 479, 84 (2011).
  * (10) M. Widmann, S. Y. Lee, T. Rendler, N. T. Son, H. Fedder, S. Paik, L. P. Yang, N. Zhao, S. Yang, I. Booker, A. Denisenko, M. Jamali, S. A. Momenzadeh, I. Gerhardt, T. Ohshima, A. Gali, E. Janzén, and J. Wrachtrup,  Coherent control of single spins in silicon carbide at room temperature,  Nat. Mater. 14, 164 (2015).
  * (11) A. L. Falk, B. B. Buckley, G. Calusine, W. F. Koehl, V. V. Dobrovitski, A. Politi, C. A. Zorman, P. X.-L. Feng, and D. D. Awschalom,  Polytype control of spin qubits in silicon carbide,  Nat. Commun. 4, 1819 (2013).
  * (12) P. G. Baranov, A. P. Bundakova, A. A. Soltamova, S. B. Orlinskii, I. V. Borovykh, R. Zondervan, R. Verberk, and J. Schmidt,  Silicon vacancy in SiC as a promising quantum system for single-defect and single-photon spectroscopy,  Phys. Rev. B 83, 125203 (2011).
  * (13) V. A. Soltamov, A. A. Soltamova, P. G. Baranov, and I. I. Proskuryakov,  Room temperature coherent spin alignment of silicon vacancies in 4 _H_ \- and 6 _H_ -SiC,  Phy. Rev. Lett. 108, 226402 (2012).
  * (14) G. Wolfowicz, C. P. Anderson, B. Diler, O. G. Poluektov, F. J. Heremans, and D. D. Awschalom,  Vanadium spin qubits as telecom quantum emitters in silicon carbide,  Sci. Adv. 6, 1192 (2020).
  * (15) F. F. Yan, A. L. Yi, J. F. Wang, Q. Li, P. Yu, J. X. Zhang, A. Gali, Y. Wang, J. S. Xu, X. Ou, C. F. Li, and G. C. Guo,  Room-temperature coherent control of implanted defect spins in silicon carbide,  npj Quantum Information 6, 38 (2020).
  * (16) H. Matsunami,  Current SiC technology for power electronic devices beyond Si,  Microelectronic Engineering 83, 2 (2006).
  * (17) J. R. Jenny, D. P. Malta, M. R. Calus, S. G. Müller, A. R. Powell, V. F. Tsvetkov, H. M. Hobgood, R. C. Glass, and C. H. Carter,  Development of large diameter high-purity semi-insulating 4 _H_ -SiC wafers for microwave devices,  Mater. Sci. Forum 457-460, 35 (2004).
  * (18) N. T. Son, P. Carlsson, J. ul Hassan, E. Janzén, T. Umeda, J. Isoya, A. Gali, M. Bockstedte, N. Morishita, T. Ohshima, and H. Itoh,  Divacancy in 4 _H_ -SiC,  Phys. Rev. Lett. 96, 055501 (2006).
  * (19) Q. Li, J. F. Wang, F. F. Yan, J. Y. Zhou, H. F. Wang, H. Liu, L. P. Guo, X. Zhou, A. Gali, Z. H. Liu, Z. Q. Wang, K. Sun, G. P. Guo, J. S. Tang, H. Li, L. X. You, J. S. Xu, C. F. Li, and G. C. Guo,  Room temperature coherent manipulation of single-spin qubits in silicon carbide with a high readout contrast,  Nat. Sci. Rev., doi: 10.1093/nsr/nwab122 (2021).
  * (20) D. Carbonera,  Optically detected magnetic resonance (ODMR) of photoexcited triplet states,  Photosynthesis Research 102, 403 (2009).
  * (21) T. Plakhotnik, M. W. Doherty, J. H. Cole, R. Chapman, and N. B. Manson,  All-optical thermometry and thermal properties of the optically detected spin resonances of the nv-center in nanodiamond,  Nano Lett. 14, 4989 (2014).
  * (22) C. C. Li, M. Gong, X. D. Chen, S. Li, B. W. Zhao, Y. Dong, G. C. Guo, and F. W. Sun,  Temperature dependent energy gap shifts of single color center in diamond based on modified Varshni equation,  Diamond Relat. Mater. 74, 119 (2017).
  * (23) D. M. Toyli, D. J. Christle, A. Alkauskas, B. B. Buckley, C. G. Van de Walle, and D. D. Awschalom,  Measurement and control of single nitrogen-vacancy center spins above 600 K,  Phys. Rev. X 2, 031001 (2012).
  * (24) F. F. Yan, J. F. Wang, Q. Li, Z. D. Cheng, J. M. Cui, W. Z. Liu, J. S. Xu, C. F. Li, and G. C. Guo,  Coherent Control of Defect Spins in Silicon Carbide above 550 K,  Phys. Rev. Appl. 10, 044042 (2018).
  * (25) Y. Zhou, J. F. Wang, X. M. Zhang, K. Li, J. M. Cai, and W. B. Gao,  Self-protected thermometry with infrared photons and defect spins in silicon carbide,  Phys. Rev. Appl. 8, 044015 (2017).
  * (26) G. Kucsko, P. C. Maurer, N. Y. Yao, M. Kubo, H. J. Noh, P. K. Lo, H. Park, and M. D. Lukin,  Nanometre-scale thermometry in a living cell,  Nature 500, 54 (2013).
  * (27) J.F. Ziegler, M.D. Ziegler, and J.P. Biersack,  SRIM – the stopping and range of ions in matter,  Nucl. Instrum. Methods Phys. Res., Sect. B 268, 1818 (2010).
  * (28) Q. Li, J. F. Wang, F. F. Yan, Z. D. Cheng, Z. H. Liu, K. Zhou, L. P. Guo, X. Zhou, W. P. Zhang, X. X. Wang, W. Huang, J. S. Xu, C. F. Li, and G. C. Guo,  Nanoscale depth control of implanted shallow silicon vacancies in silicon carbide,  Nanoscale. 11, 20554 (2019).
  * (29) I. I. Rabi, J. R. Zacharias, S. Millman, and P. Kusch,  A new method of measuring nuclear magnetic moment,  Phys. Rev. 53, 318 (1938).
  * (30) X. Rong, J. P. Geng, F. Z. Shi, Y. Liu, K. B. Xu, W. C. Ma, F. Kong, Z. Jiang, Y. Wu and J. F. Du,  Experimental fault-tolerant universal quantum gates with solid-state spins under ambient conditions,  Nat. Commun. 6, 8748 (2015).
  * (31) J. F. Barry, J. M. Schloss, E. Bauch, M. J. Turner, C. A. Hart, L. M. Pham, and R. L. Walsworth,  Sensitivity Optimization for NV-Diamond Magnetometry,  Rev. Mod. Phys. 92, 015004 (2020).
  * (32) D. P. DiVincenzo,  The Physical Implementation of Quantum Computation,  Fortschritte Der Physik 48, 771 (2000).
  * (33) L. P. Yang, C. Burk, M. Widmann, S. Y. Lee, J. Wrachtrup, and N. Zhao,  Electron Spin Decoherence in Silicon Carbide Nuclear Spin Bath,  Phys. Rev. B 90, 241203(R) (2014).
  * (34) J. R. Maze, J. M. Taylor, and M. D. Lukin,  Electron spin decoherence of single nitrogen-vacancy defects in diamond,  Phys. Rev. B 78, 094303 (2008).
  * (35) H. Seo, A. L. Falk, P. V. Klimov, K. C. Miao, G. Galli, and D. D. Awschalom,  Quantum decoherence dynamics of divacancy spins in silicon carbide,  Nat. Commun. 4, 1819 (2013).
  * (36) J. S. Embley, J. S. Colton, K. G. Miller, M. A. Morris, M. Meehan, S. L. Crossen, B. D. Weaver, E. R. Glaser, and S. G. Carter,  Electron spin coherence of silicon vacancies in proton-irradiated 4 _H_ -SiC,  Phys. Rev. B. 95, 045206 (2017).
  * (37) N. BarGill, L. M. Pham, A. Jarmola, D. Budker, and R. L. Walsworth,  Solid-state electronic spin coherence time approaching one second,  Nat. Commun. 4, 1743 (2012).
  * (38) A. M. Tyryshkin, S. A. Lyon, A. V. Astashkin, and A. M. Raitsimring,  Electron spin relaxation times of phosphorus donors in silicon,  Phys. Rev. B 68, 193207 (2003).
  * (39) A. Norambuena, E. Muoz, H. T. Dinani, A. Jarmola, P. Maletinsky, D. Budker, and J. R. Maze,  spin-lattice relaxation of individual solid-state spins,  Phys. Rev. B 97, 094304 (2018).
  * (40) K. N. Shrivastava,  Theory of Spin-Lattice Relaxation,  Phys. status solidi 117, 437 (1983).
  * (41) G. D. Lange, Z. H. Wang, D. Risté, V. V. Dobrovitski, and R. Hanson,  Universal Dynamical Decoupling of a Single Solid-State Spin from a Spin Bath,  Science 330, 60 (2010).
  * (42) A. Bourassa, C. P. Anderson, K. C. Miao, M. Onizhuk, H. Ma, A. L. Crook, H. Abe, J. U.-Hassan, T. Ohshima, N. T. Son, G. Galli, and D. D. Awschalom,  Entanglement and control of single nuclear spins in isotopically engineered silicon carbide,  Nat. Mater. 19, 1319 (2020).
  * (43) N. Morioka, C. Babin, R. Nagy, I. Gediz, E. Hesselmeier, D. Liu, M. Joliffe, M. Niethammer, D. Dasari, V. Vorobyov, R. Kolesov, R. Stöhr, J. U.-Hassan, N. T. Son, T. Ohshima, P. Udvarhelyi, G. Thiering, A. Gali, J. Wrachtrup, and F. Kaiser,  Spin-controlled generation of indistinguishable and distinguishable photons from silicon vacancy centres in silicon carbide,  Nat. Commun. 11, 2516 (2020).
  * (44) S. Y. Paik, S. Y. Lee, W. J. Baker, D. R. Mccamey, and C. Boehme,  $T_{1}$ and $T_{2}$ spin relaxation time limitations of phosphorous donor electrons near crystalline silicon to silicon dioxide interface defects,  Phys. Rev. B 81, 075214 (2010).
  * (45) T. Rosskopf, A. Dussaux, K. Ohashi, M. Loretz, R. Schirhagl, H. Watanabe, S. Shikata, K. M. Itoh, and C. L. Degen,  Investigation of surface magnetic noise by shallow spins in diamond,  Phys. Rev. Lett. 112, 147602 (2014).
  * (46) A. Jarmola, V. M. Acosta, K. Jensen, S. Chemerisov, and D. Budker,  Temperature- and Magnetic-Field-Dependent Longitudinal Spin Relaxation in Nitrogen-Vacancy Ensembles in Diamond,  Phys. Rev. Lett. 108, 197601 (2012).