# Enhanced temperature sensing by multi-mode coupling in an on-chip
microcavity system

###### Abstract

The microcavity is a promising sensor platform, any perturbation would disturb
its linewidth, cause resonance shift or splitting. However, such sensing
resolution is limited by the cavity’s optical quality factor and mode volume.
Here we propose and demonstrate in an on-chip integrated microcavity system
that resolution of a self-referenced sensor could be enhanced with multi-mode
coupling. In experiments, inter-mode coupling strength is carefully optimized
with a pulley waveguide and observed a resolution improvement of nearly $3$
times in frequency domain. While experiencing small refractive index change
tuned by temperature, mode-coupled system shows a $7.2$ times sensitivity
enhancement that is than the uncoupled system on the same chip and a very
significant lineshape contrast ratio change as great reference for minor
frequency shifts. This approach will help design microcavity sensors to
improve detection sensitivity and resolution under limited manufacture
precision.

###### keywords:

Micro-optical device, Optical sensing and sensors, Mode coupling, Integrated
optics, Resonant modes

Xueyi Wang‡ Tingge Yuan‡ Jiangwei Wu Yuping Chen* Xianfeng Chen

‡These authors contributed equally to this work.

X. Wang, T. Yuan, J. Wu, Prof. Y. Chen, Prof. X. Chen
State Key Laboratory of Advanced Optical Communication Systems and Networks
School of Physics and Astronomy
Shanghai Jiao Tong University
Shanghai 200240, China
Email<EMAIL_ADDRESS>

Prof. X. Chen
Shanghai Research Center for Quantum Sciences
Shanghai 201315, China

Prof. X. Chen
Collaborative Innovation Center of Light Manipulations and Applications
Shandong Normal University
Jinan 250358, China

##  1 Introduction

Optical microcavity as one of the building blocks of photonic integrated
circuit has enabled a variety of applications including nonlinear optics[1-2],
low-threshold laser[3-4] and single molecule detection[5-12], its small mode
volume and high quality factor ($Q$ factor). Especially in unlabeled
sensing[13-19] or environmental monitoring[20-25] as a great supplement for
medical and environmental research. On the other hand, decrease in the
cavity’s mode volume would increase radiation losses that is no longer
neglectable[26] causing drop in its $Q$ factor. Overcoming such limitation by
introducing new principles to microcavity systems thus become urgent. There
have been works implementing microcavity lasers[27-28] to enhance light-matter
interaction or by utilizing opto-mechnical coupling[29] that boosts sensing
resolution by magnitudes.

Among all of the brand new solutions, one of the most approachable is by
introducing mode coupling into the system[30], for that it causes little extra
fabrication or experiment difficulties. When the two coupled modes are at
weak-coupling region[31], their coherent interaction can optimize the
spectrum’s lineshape for efficient sensing[32]. Within one single cavity, the
coupling condition can be satisfied by utilizing modes in different
polarization of a micro-toroidal cavity[33] or by applying UV curable adhesive
onto micro-bottle resonator to create a lossy mode[34], which achieved a 4.3
times refractive index change sensitivity amplification through its coupling
with another discrete mode. Meanwhile micro-ring resonator is the ideal
platform for on-chip integration, with as simple as a built-in Fabry–Pérot
(F-P) cavity on its coupled waveguide[35], multi-mode coupling between
cavities could be achieved. Such structure was first manufactured with polymer
platform[36] that increased the sensitivity of solution refractive index by
its sharp resonance slopes, and with a silicon-on-insulator chip as well[37]
which realized a tunable lineshape fitting a variety of applications.
Recently, mode coupling has also been controlled by scatterers to function at
its exceptional point showing possibility for unprecedented sensitivity[38].
Thus mode coupling could be a handy improvement to already widely studied
microcavity sensors.

In this work, we propose a design method to improve the resolution of
microcavity sensors through multi-mode coupling with a compact, on-chip
integrated micro-cavity system. Based on a waveguide to micro-racetrack
structure supporting three resonance modes simultaneously and a pulley coupler
with careful geometrical optimization, our design allows efficient and
distinct inter-mode coupling at $1520$ nm to $1555$ nm band for both racetrack
quasi-TE and TM modes leading to frequency shifts and sharp lineshape, which
helps to distinguish two modes during self-referenced sensing and breaks the
sensitivity’s dependence on the $Q$ factor that microcavity sensors always
suffer from. In frequency domain we achieved 3 times enhancement in resolution
and a sensitivity of 44 $\rm pm^{\circ}C^{-1}$ that is 7.2 times higher than
the uncoupled structure on the same chip, with lineshape contrast ratio (LCR)
of 24.1 times enhancement in the same time which function as great reference
for minor turbulence. Our proposed approach will benefit applications in
optical sensors that require integration and high sensitivity probing weak
signals under a compromised fabrication efficiency.

##  2 Theory

Conventionally when two modes are weakly coupled, for instance one discrete
mode and one continuous mode, the discrete mode would experience a frequency
shift and linewidth sharpening determined by their detuned wavelength and
coupling strength[31]. While the coupling includes two discrete modes
simultaneously, it is very likely that they will experience different shifts
for that they possess distinct coupling strength and eigenfrequencies. Thus by
controlling the composition of three modes we could manipulate the relative
frequency difference of them after coupling happened, that in certain
scenarios would help us to distinguish two discrete modes with higher
resolution. Here we first introduce the theory and how it compose with our on-
chip system.

In the waveguide micro-ring resonator (WGMRR) three modes co-exist in the
system as in Figure 1a, one waveguide (WG) mode reflected by built-in
gratings, two micro-ring resonance (MRR) modes with quasi-TE and TM
polarization. They possess different coupling efficiency $\kappa_{j}$,
internal loss $\gamma_{j}$ and resonant frequency separately $\omega_{j}$
($j=0$ for WG mode and $j=1,2$ for MRR quasi-TE, TM modes respectively as
shown in Figure 1a). The system Hamiltonian should be,

$H_{SYS}=\sum_{j,j=0,1,2}\hbar\omega_{i}a_{i}^{{\dagger}}a_{i}+i\hbar\kappa_{1}(a_{0}^{{\dagger}}a_{1}-a_{1}^{{\dagger}}a_{0})+i\hbar\kappa_{2}(a_{0}^{{\dagger}}a_{2}-a_{2}^{{\dagger}}a_{0})-ig\hbar(a_{1}^{{\dagger}}a_{2}+a_{2}^{{\dagger}}a_{1}).$
(1)

Figure 1: a) Scheme of the WGMRR system with WG mode in blue and MRR quasi-TE
and quasi-TM in red and black respectively. The insets are calculated with
conformal transformation[39] indicating the polarization of $\rm TE_{00}$ and
$\rm TM_{00}$ in straight and curved waveguides (see also Supporting
Information section 1.1). Frquency shift b) and linewidth sharpening c) for
different MRR modes caused by coupling, insets are transmission spectrum of
quasi-TE, TM modes after coupling.

In which $a_{j}^{{\dagger}}$ $(a_{j})$ are photon creation (annihilation)
operators. For the oblique nature of the side walls ridge waveguide and
birefringence of X-cut LN (Lithium Niobate), the quasi-TE (TM) modes are not
perfectly parallel(vertical) to the substrate plane (see inset of Figure 1a
and Supporting Information section 1.1) which causes them to couple with a
coefficient $g$[40]. At the same time enabling WG $\rm TE_{00}$ photon $a_{0}$
to generate MRR TM mode $a_{2}$ with different polarization. Similarly, at the
time $a_{2}$ couples back into the waveguide it is projected into TE polarized
beams that lead to its coherent interaction with $a_{0}$ lights.

Under the first Markov approximation
$\kappa_{0}^{2}(\omega)=\kappa_{0}/2\pi$[41], we approach the Langevin
equations of motion as follow,

$\left(\begin{array}[]{c}\dot{a_{0}}\\\ \dot{a_{1}}\\\
\dot{a_{2}}\end{array}\right)=\left(\begin{array}[]{ccc}-i\omega_{0}-\frac{\kappa_{0}+\gamma_{0}}{2}&\kappa_{1}&\kappa_{2}\\\
-\kappa_{1}&-i\omega_{1}-\frac{\kappa_{1}+\gamma_{1}}{2}&-g\\\
-\kappa_{2}&-g&-i\omega_{2}-\frac{\kappa_{2}+\gamma_{2}}{2}\end{array}\right)\left(\begin{array}[]{c}a_{0}\\\
a_{1}\\\
a_{2}\end{array}\right)+\left(\begin{array}[]{c}\sqrt{\kappa_{0}}a_{IN}\\\
0\\\ 0\end{array}\right),$ (2)

$a_{IN}=\sqrt{2P_{IN}\kappa_{0}/\hbar\omega_{0}}$ is the input amplitude in TE
polarization and so as $a_{0}$. Then we solve Equation (2) in frequency
domain,

$(\omega-\widetilde{\omega}_{0})a_{0}=-i\kappa_{1}a_{1}-i\kappa_{2}a_{2}-i\sqrt{\kappa_{0}}a_{IN}$
(3) $(\omega-\widetilde{\omega}_{1})a_{1}=iga_{2}+i\kappa_{1}a_{0},$ (4)
$(\omega-\widetilde{\omega}_{2})a_{2}=iga_{1}+i\kappa_{2}a_{0},$ (5)

in which
$\widetilde{\omega}_{j}=\omega_{j}-i\frac{\kappa_{j}+\gamma_{j}}{2}(j=0,1,2)$
is the complex eigenfrequency. Equation (3) implies that the output field is
made up by the coherent superposition of $a_{1},a_{2}$ to $a_{0}$ mode with
coupling constant $\kappa_{1}$ and $\kappa_{2}$ respectively. For that inter-
modal coupling coefficient $g$ within the MRR is relatively small, output
amplitude $a_{OUT}=\sqrt{\kappa_{0}}a_{0}-a_{IN}$ approximates to,

$a_{OUT}=\kappa_{0}\xi a_{IN}-a_{IN},$ (6)
$\xi=\frac{i(\omega-\widetilde{\omega}_{1})(\omega-\widetilde{\omega}_{2})}{(\omega-\widetilde{\omega}_{0})(\omega-\widetilde{\omega}_{1})(\omega-\widetilde{\omega}_{2})-\kappa_{1}^{2}(\omega-\widetilde{\omega}_{2})-\kappa_{2}^{2}(\omega-\widetilde{\omega}_{1})},$
(7)

the transmission spectrum should be,

$T=\left|\frac{a_{OUT}}{a_{IN}}\right|^{2}\approx 1-2\kappa_{0}\xi,$ (8)

where $|\kappa_{0}\xi|^{2}$ is ignored in above equation. Solving the
denominator of $\xi$ there are (details in Supporting Information 2),

$\widetilde{\omega}_{j\pm}=(\widetilde{\omega}_{0}+\widetilde{\omega}_{j}\pm\delta_{j})/2,\delta_{j}^{2}=(\widetilde{\omega}_{0}-\widetilde{\omega}_{j})^{2}+4\kappa_{j}^{2},$
(9)

the eigenfrequencies of MRR are shifted by complex frequencies
$\widetilde{\Delta}_{\pm
j}=\widetilde{\omega}_{j\pm}-\widetilde{\omega}_{j}=(\widetilde{\omega}_{0}-\widetilde{\omega}_{1}\pm\delta_{j})/2$,
$(+:\omega_{0}-\omega_{j}>0,-:\omega_{0}-\omega_{j}<0)$ respectively, in which
$\Delta_{\omega}=Re(\widetilde{\Delta}_{\pm j})$ stands for shifts in
frequency and $\Delta_{\kappa+\gamma}=-Im(\widetilde{\Delta}_{\pm j})$ for
changes in linewidth. Consequently, MRR modes experience a red (blue) shift if
they are blue (red) detuned to $\omega_{0}$ and a linewidth reduction either
way as shown in Figure 1b and c. For that quasi-TE mode has greater $\kappa$
that leads to bigger $|\Delta_{\omega}|$ and $|\Delta_{\kappa+\gamma}|$ for
$a_{1}$.

Consequently, if two MRR modes locate across the continuum mode eigenfrequency
they would be ”pulled closer” as their frequencies are shifted towards each
other
$|\Delta^{\prime}_{12}|=|\Delta_{12}|-\left[|Re(\delta_{1})|+|Re(\delta_{2})|\right]/2$
as case I in Figure 2a. Or they could be ”pushed apart” if coupled to two
different continuum modes as the cases in Figure 2a II. that
$|\Delta^{\prime}_{12}|=|\Delta_{12}|+\left[|Re(\delta_{1})|+|Re(\delta_{2})|\right]/2$.
In the case when their frequency differs even less, within one side of the
background’s FSR (free spectrum range), then they are ”pushed apart” when mode
with larger $\kappa$ locates closer to $\omega_{0}$ only then
$|\Delta^{\prime}_{12}|=|\Delta_{12}|+\left[|Re(\delta_{1})|-|Re(\delta_{2})|\right]/2$
could be enlarged (as the case in Figure 2a III, the width of modes’ stripes
indicates their relative coupling strength), and ”pulled closer” if in the
opposite case. Under the circumstances when multi-mode coupling leads to
”pushed apart”, it enhances the observation resolution of two adjacent MRR
modes, which is ideal for sensing applications or mode measurements. From our
following experiments, clearly in Figure 2b, c and d the transmission spectrum
of the WGMRR and WG (in black and blue) has strongly asymmetric dips and peaks
indicate multi-mode coupling happening in the above (I to III) scenarios.

Figure 2: a) Schematic of the proposed mode coupling involving three modes
arranged in 3 fashions leading to different mode resolution: when MRR modes
locate I across background eigenmode then they are ”pulled closer”, II over
different background eignmodes and III within one side of background mode
while mode with larger $\kappa$ is closer to $\omega_{0}$ when they are
”pushed apart”, wider stripe indicates stronger coupling strength of the
eigenmode to the background mode. Transmission spectrum of the WGMRR and WG in
black and blue which satisfies the above scheme I b), II c), and III d)
respectively.

##  3 Results

###  3.1 Exprimental setup

Figure 3: a) Experiment setup. TL: Tunable Laser, PC: polarization controller,
PD: photon detector, OSC: oscillator, TEC: thermoelectric cooler, electrical
wire and optical fiber in black and yellow respectively. Microscopic pictures
of the single waveguide b), gratings c), coupling region d), racetrack e) and
ring f) resonators.

The experimental setup is in Figure 3a, it involves a tunable infrared laser
source (New Focus TLB-6728) which is adjusted by a PC (polarization
controller) and then coupled into the chip through built-in gratings on the
waveguide. After the output port a PD (photodetector) collects the
transmission information that is then shown on the OSC (oscilloscope).
Meanwhile the chip is loaded on a TEC (thermoelectric cooler) stage with a
precision up to $0.01$ degree Celsius. From this setup any multi-mode coupling
effect is acquired from the the lineshape of transmission spectrum.

To achieve multi-mode coupling and compare its effects, we integrate a
waveguide, a waveguide to micro-ring and a waveguide to micro-racetrack all
manufactured on one single X-cut LNOI (lithium niobate on insulator) chip with
standard electron beam lithography and plasma-reactive etching (see
fabrication details in Supporting Information section 3). As marked in insets
of Figure 1a, the chip thickness is $h=0.6$ $\mu m$, while top width of
waveguide is $w=1$ $\mu m$, thickness $t=0.38$ $\mu m$ and side wall angle
$\theta=60^{\circ}$, which are ideal to support only fundamental modes. The
radius of both micro-racetrack and micro-ring is the same $R=129.03$ $\mu m$
while the racetrack contains an extra straight waveguide of $82.54$ $\mu m$.
At the coupling region Figure 3c, the waveguide width shrinks down to $0.8$
$\mu m$ with a gap of $G=0.6$ $\rm\mu m$ to achieve sufficient evanescent
field coupling, next we calculate such coupling strength and analyse its
effect in detail.

###  3.2 Multi-mode coupling in frequency domain

From the system Hamiltonian we could tell that the inter-modal coupling
strength is determined by the mode coupling efficient $\kappa$ at the pulley
waveguide, and it can be calculated with temporal perturbation theory as[42]
(details in Supporting Information section 1.2),

$\kappa=\int_{-\psi_{0}}^{\psi_{0}}\left[\frac{i\omega}{4}\int_{0}^{R+G}\int_{0}^{t}\left(\epsilon-\epsilon_{0}\right)\mathbf{E}_{WG}\cdot\mathbf{E}_{MRR}rdrdz\right]e^{i\varphi}d\psi,$
(10)

in which $\mathbf{E}_{WG}\cdot\mathbf{E}_{MRR}$ corresponds to the mode
overlap of normalized WG and MRR fields, permittivity $\epsilon$ can be
obtained from the electric field and is shown in the mode’s effective
refractive index ($n_{eff}$) in Figure 4a  and b. And
$\varphi=k_{0}n_{WG}R_{WG}\psi-m\psi$, $m$ is the radial mode number of
resonance mode inside cavity, $\varphi$ reflects the phase mismatch between
the waveguide and cavity mode that has major impact on $\kappa$. Meanwhile
X-cut LN is a anisotropic material as $n_{eff}$ shifts across azimuth angle
$\psi$ infecting phase matching condition along its way, thus it introduces
change in $\kappa$ proportional to $sinc(\varphi)$, see Figure 4c, where
$\kappa_{TE}$ increases and $\kappa_{TM}$ decrease as $\psi$ grows. It offers
another degree of freedom to manipulate coupling strength between different
polarization with the angle and length of pulley coupling scheme. Also in
frequency domain (Figure 4d) $\kappa$ degenerates significantly at longer
input wavelength due to non-ideal phase matching and mode overlaping, in the
following experiments we clearly observed weaker mode coupling at longer
wavelength as shown in Figure 4g. Reminding that the pulley waveguide needs to
be carefully designed to achieve sufficient coupling at target work waveband.
In our design, $\kappa$ is set to generate sufficient coupling for both
polarization across C-band at the same time differs with nearly one degree of
magnitude so they appear with diverse mode shifts during coupling to be
distinguished.

Figure 4: The simulated effective refractive index ($n_{eff}$) at different
angles a) on the X-cut LN chip and across resonance wavelength b). Coupling
coefficient $lg(\kappa)$ and $sinc(\varphi)$ at different angles c) on the
X-cut LN chip and across resonance wavelength d). e)Transmission spectrum of
the WG-micro-racetrack and the sole waveguide (in blue) when there is multi-
mode coupling. g)Transmission spectrum of the WG-micro-ring and the sole
waveguide (in blue) when there is no multi-mode coupling. Measured wavelength
difference $|\Delta_{12}|$, FWHM $2/(w_{1}+w_{2})$ and their product
$2|\Delta_{12}|/(w_{1}+w_{2})$ declaring the degree of mode separation of TM
mode and its closest TE mode in frequency domain for coupled f) and uncoupled
h) system, the red circle marks where inter-modal separation is significantly
enlarged by coupling.

With above system we first obtain the transmission spectrum of the micro-
racetrack (Figure 4e) the cavity modes preserved a strongly asymmetric
lineshape due to the coupling with WG modes (the transmission spectrum of
waveguide is in blue). The effective $Q$ factor calculated by applying Lorentz
fitting of the dip reaches $Q_{E}=11.07k$ and $Q_{M}=14.09k$ for TE and TM
respectively. With calculated coupling $Q_{CE}=17.10k$, $Q_{CM}=62.31k$ which
ends up in coupling coefficients $\kappa_{E}=2.67\times 10^{5}$,
$\kappa_{M}=1.40\times 10^{5}$ slightly larger than the simulation in Figure
4d supposedly caused by coupling region $\psi>15^{\circ}$ which has a larger
gap but still allows evanescent field coupling. Compared to their intrinsic
$Q$ factors $Q_{0E}=27.52k$, $Q_{0M}=18.21k$, the TE modes are over-coupled
leading to strong interference by WG modes, theoretically it causes mode shift
up to $72.52$ kHz as shown in Figure 1b. Experimentally, due to the different
FSRs of TE and TM modes, the wavelength difference between two closest modes
$|\Delta_{12}|$ would first increase and then decrease to zero during a
certain wavelength range as the insets of both Figure 4f and h. In a coupled
system, $|\Delta_{12}|$ and modes’ FWHM (full width half maximum, $w_{1}$,
$w_{2}$ for TE and TM respectively) were tuned by coupling strength and
relative background phase as in Figure 4f, determined by
$\widetilde{\omega}_{0}-\widetilde{\omega}_{1}$ according Equation (9). The
measured degree of mode separation $\frac{|\Delta_{12}|}{(w_{1}+w_{2})/2}$ is
enhanced for TM modes with radial mode number $1153$ to $1158$ (marked in red
in Figure 4f), satisfying either scenario II or III from Figure 2a, nearly 3
times larger compared to their siblings (mode number $1153$ to $1163$).

Meanwhile in the uncoupled system with micro-ring resonator (Figure 3e), whose
loaded $Q$ factors are $Q_{E}^{\prime}=22.61k$, $Q_{M}^{\prime}=10.76k$ and
coupling $Q_{CE}^{\prime}=137.37k$, $Q_{CM}^{\prime}=203.82k$ that is
significantly under-coupled. Calculated from its spectrum Figure 4g the mode
separation $\frac{|\Delta_{12}|}{(w_{1}+w_{2})/2}$ is proportional to
$|\Delta_{12}|$ despite of the fluctuation in FWHM caused by coupling depth as
in Figure 4h showing no signs of resolution enhancement nor inter-mode
coupling. Next we run tests over their sufficiency of detecting minor
temperature shifts for both systems.

###  3.3 Temperature sensing enhanced by multi-mode coupling

Figure 5: a) Transmission spectrum the WGMRR with multi-mode coupling leading
to ”pushing apart”. b) Transmission spectrum the WGMRR without mode coupling
by adjusting input polarization and wavelength. c) Measured wavelength
separation $|\Delta_{12}|$ of a) and b)(dashed) at different temperature. d)
Calculated lineshape contrast ratio
$\frac{(I_{1}-I_{0})+(I_{2}-I_{0})}{I_{1}+I_{2}}$ of a) and b)(dashed), under
multi-mode coupling their lineshapes experience a contrast shift $24.1$ times
greater.

In previous works utilizing double mode coupling such as fano resonance for
sensing purpose, apart from the sharp asymmetric lineshape that declines the
effective mode linewidth, those two modes are self-referenced so the system
accuracy is no longer limited by experiment equipment[42]. According to above
analysis, MRR quasi-TE and TM modes would sufficiently form a self-referenced
sensing system while coupling with WG mode further improves its capacity. Here
to compare the sensitivity for weak turbulence of a mode-coupled and an
uncoupled system, we keep the laser source sweeping in wavelength and slightly
adjust the stage’s temperature control which tunes $n_{eff}$ of the MRR
consequently. At around $1550$ nm wavelength the micro-racetrack forms an EIT-
like spectrum in Figure 5a, the three groups of modes experienced different
shifts due to their coupled background phase, the gap between two closest
modes $|\Delta_{12}|$ shifts as large as $44$ $\rm pm^{\circ}C^{-1}$ which is
$7.2$ times larger than the uncoupled system in Figure 5b that has very
similar mode gaps in the first place, as shown in Figure 5c.

On the other hand, a mode-coupled system possess a lineshape sensitive to its
background phase[43]. In our setup the background modes have a linewidth of
$0.23$ nm and FSR of $0.45$ nm (and it is also insensitive to thermal changes,
see Supporting Information section 4), thus even a minor shift at $\sim\rm pm$
would lead to observable changes to the spectrum lineshape. That experiences a
change in its LCR
$\left(\frac{(I_{1}-I_{0})+(I_{2}-I_{0})}{I_{1}+I_{2}}\right)$ as large as
$6.46\times 10^{-3}$ $\rm pm^{-1}$, while the according uncoupled modes only
changes $0.27\times 10^{-3}$ $\rm pm^{-1}$ as in Figure 5d. So there is a
maximum $24.1$ times enhancement compared with the uncoupled spectrum at
similar temperature working as another crucial criterion for minor turbulence
in the system’s $n_{eff}$.

It is still worth noting that our design did not reach the best performance
among all of the microcavity thermal sensors. Apart from the fact that our on-
chip microcavities were built large in the first place to compensate the WG
mode’s FSR, the pulley waveguide was designed to reach over-coupling condition
which lead to decrease in MRR modes’ $Q$ factor and coupling depth. Taking
those information under consideration, a further improvement should utilize
the coupling of multiple high-$Q$ factor modes either by connecting several
micro-resonators[30] or by optimizing cavities that support a number of
resonance modes.

##  4 Conclusion

In this paper, we have in depth analysed dimensions of multi-mode coupling,
revealing the intrinsic connection between frequency, $Q$ factor enhancement
and modes’ composition. Based on theoretical discussions and experiments
demonstrated with an on-chip integrated WGMRR, we have achieved higher
resolution in frequency domain and improved sensitivity as a self-referenced
sensor with micro-racetrack’s quasi-TE and TM modes coupled to waveguide
resonance mode at the same time. Here we have confined our discussions of
mode-coupled sensors within the fundamental case, practiced between discrete
and continuum modes, and it is expected that a natural extension to more
sophisticated configurations with several high-$Q$ modes or even multiple
exceptional points would reach much higher magnitudes of enhancement for
microcavity sensors without over-investing into manufacture techniques. Though
our design was practiced with a LNOI micro-racetrack, we believe that this
method can be applied to any cavity-based sensors or other material platforms
with great integration capability.

Supporting Information

Supporting Information is available from the Wiley Online Library or from the
author.

Acknowledgements

The authors would like to acknowledge support from National Natural Science
Foundation of China (Grant Nos. 12134009) and SJTU (No. 21X010200828).

Conflict of Interest

The authors declare no conflict of interest.

Data Availability Statement

The data that support the findings of this study are available from the cor-
responding author upon reasonable request.

Keywords

Micro-optical device, Optical sensing and sensors, Mode coupling, Integrated
optics, Resonant modes

References

[1] T. J. Kippenberg, S. M. Spillane, K. J. Vahala, Phys. Rev. Lett. 2004, 93
083904

[2] X. Ye, S. Liu, Y. Chen, Y. Zheng, X. Chen, Opt. Lett. 2020, 45, 2 523.

[3] J. W. B. Z. Y. C. . X. C. YiAn Liu, XiongShuo Yan, Sci. China Phys. Mech.
Astron. 2021, 64, 64234262.

[4] X. Liu, X. Yan, Y. Liu, H. Li, Y. Chen, X. Chen, Opt. Lett. 2021, 46, 21
5505.

[5] F. Vollmer, L. Yang, Nanophotonics 2012, 1, 3-4 267. [6] Y. Zhi, X.-C. Yu,
Q. Gong, L. Yang, Y.-F. Xiao, Advanced Materials 2017, 29, 12 1604920.

[7] X. Jiang, A. Qavi, S. Huang, L. Yang, Matter 2020, 3, 2 371.

[8] M. R. Foreman, J. D. Swaim, F. Vollmer, Adv. Opt. Photon. 2015, 7, 2 168.

[9] S.-J. Tang, M. Zhang, J. Sun, J.-W. Meng, X. Xiong, Q. Gong, D. Jin, Q.-F.
Yang, Y.-F. Xiao, Nature Photonics 2023, 1–6.

[10] F. Vollmer, S. Arnold, D. Keng, Proceedings of the National Academy of
Sciences 2008, 105, 5220701.

[11] X. Y.-F. L. L. H. L. C. D.-R. Zhu Jiangang, Ozdemir Sahin Kaya, Y. Lan,
Nature Photonics 2010, 4, 1 46.

[12] T. Lu, H. Lee, T. Chen, S. Herchak, J.-H. Kim, S. E. Fraser, R. C.
Flagan, K. Vahala, Proceedings of the National Academy of Sciences 2011, 108,
15 5976.

[13] S. Frustaci, F. Vollmer, Current Opinion in Chemical Biology 2019, 51 66,
chemical Genetics and Epigenetics • Molecular Imaging.

[14] F. Vollmer, S. Arnold, D. Keng, Proceedings of the National Academy of
Sciences 2008, 105, 5220701.

[15] F. Vollmer, S. Arnold, Nature Methods 2008, 5, 7 591.

[16] R. W. Boyd, J. E. Heebner, Appl. Opt. 2001, 40, 31 5742.

[17] F. Vollmer, D. Braun, A. Libchaber, M. Khoshsima, I. Teraoka, S. Arnold,
Applied Physics Letters 2002, 80, 21 4057.

[18] W. Kim, S. K. Ozdemir, J. Zhu, F. Monifi, C. Coban, L. Yang, Opt. Express
2012, 20, 28 29426.

[19] O. Gaathon, J. Culic-Viskota, M. Mihnev, I. Teraoka, S. Arnold, Applied
Physics Letters 2006, 89, 22 223901.

[20] A. M. Armani, K. J. Vahala, Opt. Lett. 2006, 31, 12 1896.

[21] W. Kim,K.$\rm\ddot{O}$zdemir, J. Zhu, L. He, L. Yang, Applied Physics
Letters 2010, 97, 7 071111.

[22] Q. Lu, X. Chen, L. Fu, S. Xie, X. Wu, Nanomaterials 2019, 9, 3.

[23] X. Xu, W. Chen, G. Zhao, Y. Li, C. Lu, L. Yang, Light: Science and
Applications 2018, 7, 162.

[24] C.-H. Dong, L. He, Y.-F. Xiao, V. R. Gaddam, S. K. Ozdemir, Z.-F. Han,
G.-C. Guo, L. Yang, Applied Physics Letters 2009, 94, 23 231119.

[25] J. Liao, L. Yang, Light: Science and Applications 2021, 10, 1 32.

[26] C. Ciminelli, F. Dell’Olio, D. Conteduca, C. Campanella, M. Armenise,
Optics Laser Technology 2014, 59 60.

[27] L. He, K.$\rm\ddot{O}$zdemir, J. Zhu, W. Kim, L. Yang, Nature
Nanotechnology 2011, 6, 7 428.

[28] K.$\rm\ddot{O}$zdemir, J. Zhu, X. Yang, B. Peng, H. Yilmaz, L. He, F.
Monifi, S. H. Huang, G. L. Long, L. Yang, Proceedings of the National Academy
of Sciences 2014, 111, 37 E3836.

[29] W. Yu, W. C. Jiang, Q. Lin, T. Lu, Nature Communications 2016, 7, 1
12311.

[30] Y.-F. Xiao, V. Gaddam, L. Yang, Opt. Express 2008, 16, 17 12538.

[31] B. Peng, K.$\rm\ddot{O}$zdemir, W. Chen, F. Nori, L. Yang, Nature
Communications 2014, 5, 1 5082.

[32] C.-M. Chang, O. Solgaard, Opt. Express 2013, 21, 22 27209.

[33] B.-B. Li, Y.-F. Xiao, C.-L. Zou, Y.-C. Liu, X.-F. Jiang, Y.-L. Chen, Y.
Li, Q. Gong, Applied Physics Letters 2011, 98, 2, 021116.

[34] J. Liao, X. Wu, L. Liu, L. Xu, Opt. Express 2016, 24, 8 8574.

[35] S. Fan, Applied Physics Letters 2002, 80, 6 908.

[36] C.-Y. Chao, L. J. Guo, Applied Physics Letters 2003, 83, 8 1527.

[37] L. Gu, H. Fang, J. Li, L. Fang, S. J. Chua, J. Zhao, X. Gan,
Nanophotonics 2019, 8, 5 841.

[38] W. Chen, K.$\rm\ddot{O}$zdemir, G. Zhao, J. Wiersig, L. Yang, Nature
2017, 548, 7666 192.

[39] M. Heiblum, J. Harris, IEEE Journal of Quantum Electronics 1975, 11, 2
75.

[40] L. Cortes-Herrera, X. He, J. Cardenas, G. P. Agrawal, Phys. Rev. A 2021,
103 063517.

[41] C. Gardiner, P. Zoller, Quantum noise: a handbook of Markovian and non-
Markovian quantum stochastic methods with applications to quantum optics,
Springer Science and Business Media, 2004.

[42] S.-L. Chuang, Journal of Lightwave Technology 1987, 5, 1 5.

[43] U. Fano, Phys. Rev. 1961, 124 1866.


# A discrete formulation for three-dimensional winding number

Ken Shiozaki  Center for Gravitational Physics and Quantum Information, Yukawa
Institute for Theoretical Physics, Kyoto University, Kyoto 606-8502, Japan

###### Abstract

For a smooth map $g\colon X\to U(N)$, where $X$ is a three-dimensional,
oriented, and closed manifold, the winding number or the map’s degree is
defined by
$W_{3}=\frac{1}{24\pi^{2}}\int_{X}\mathrm{Tr}\left[(g^{-1}dg)^{3}\right]$. We
introduce a method to compute $W_{3}$ using a discrete approximation of $X$ so
that the result is manifestly quantized.

††preprint: YITP-24-14

Introduction— Consider a three-dimensional closed and oriented manifold $X$.
For a smooth map $g:X\to U(N)$, with $U(N)$ representing the group of $N\times
N$ unitary matrices, the winding number, an integer value, is defined by the
following expression:

$\displaystyle W_{3}[g]$ $\displaystyle=\frac{1}{2\pi}\int_{X}H\in\mathbb{Z},$
(1) $\displaystyle H$
$\displaystyle=\frac{1}{12\pi}\mathrm{Tr}\left[(g^{-1}\mathrm{d}g)^{3}\right].$
(2)

The winding number $W_{3}[g]$ is pivotal in various branches of physics,
including topological band theory, where it acts as the topological invariant
for three-dimensional superconductors with time-reversal symmetry [1, 2], and
in non-Abelian (lattice) gauge theory, appears in instanton number
calculations [3, 4]. Often in these applications, the function $g$ is defined
only on a finite set of lattice points for numerical analysis. Therefore, an
efficient numerical formulation with lattice approximation of manifolds is an
important issue.

For the first Chern number $ch_{1}$ of a line bundle with connection, a sort
of two-dimensional counterpart of the winding number, a well-established
discrete formulation with evident quantization exists [5, 6]. Furthermore,
discrete line bundles over finite simplicial complexes have been explored,
especially concerning applications in computer graphics [7]. This paper
develops a method for evaluating $W_{3}[g]$ via a discretized approximation of
$X$, ensuring the result remains manifestly quantized, provided the
approximation of $X$ is sufficiently refined relative to the spatial
variation’s length scale.

Figure 1: (a) Illustration of a $\theta$-gap. (b) Projection onto the
eigenspace between two $\theta$-gaps. (c) Smearing eigenvalues over a finite
set of vertices. In each panel, the unit circle in the complex plane is
depicted, with blue arcs representing the spectrum of the $U(N)$-valued matrix
$g(x)$ within a local region. Red lines mark the locations of $\theta$-gaps.

Formulation— Given that $H$ is a closed three-form, it is locally exact,
meaning that for a local patch, there exists a two-form $B$ such that $H=dB$.
To construct $B$ explicitly, we introduce a gap condition for elements of
$U(N)$. A matrix $g\in U(N)$ exhibits a $\theta$-gap if none of its
eigenvalues are $e^{i\theta}$ for a given real number $\theta\in[0,2\pi)$ [8],
as illustrated in Fig. 1 (a). Furthermore, we define $\log_{\theta}z$ for
nonzero complex number $z\in\mathbb{C}^{\times}$ as

$\displaystyle\log_{\theta}z=\log|z|+i\arg z,\quad\theta\leq\arg
z<\theta+2\pi.$ (3)

For two distinct $\theta$-gaps, $\theta_{1}$ and $\theta_{2}$, the following
relation holds:

$\displaystyle\log_{\theta_{1}}z-\log_{\theta_{2}}z$
$\displaystyle=\begin{cases}2\pi
i\times\mathrm{sgn}(\theta_{1}-\theta_{2})&(\min(\theta_{1},\theta_{2})<\arg
z<\max(\theta_{1},\theta_{2})),\\\ 0&(\mathrm{otherwise}).\end{cases}$ (4)

Consider $U\subset X$ as a three-dimensional subspace where $g(x)$ maintains a
$\theta$-gap for $x\in U$. Let $\gamma(x)=(u_{1}(x),\dots,u_{N}(x))\in U(N)$
be a unitary matrix diagonalizing $g(x)$, i.e.,
$g(x)=\gamma(x)\Lambda(x)\gamma(x)^{\dagger}$ with
$\Lambda(x)=\mathrm{diag}(\lambda_{1},\dots,\lambda_{N})$, where
$\lambda_{n}\in U(1)$ for $n=1,\dots,N$. The exact form $B_{\theta}$ is given
by [9]

$\displaystyle B_{\theta}$ $\displaystyle=Q+R_{\theta},$ (5) $\displaystyle Q$
$\displaystyle=\frac{1}{4\pi}\mathrm{Tr}[\gamma^{-1}d\gamma\Lambda\gamma^{-1}d\gamma\Lambda^{-1}],$
(6) $\displaystyle R_{\theta}$
$\displaystyle=\frac{1}{2\pi}\mathrm{Tr}[\log_{\theta}\Lambda(\gamma^{-1}d\gamma)^{2}].$
(7)

Note that while $Q$ is independent of $\theta$, but $R_{\theta}$ does. It is
evident that $X$ can be covered by $\\{U_{i}\\}_{i}$ such that in each patch
$U$, $g(x)$ exhibits a specific $\theta$-gap $\theta_{i}$.

The unitary matrix $\gamma$ is not unique due to the transformation
$\gamma\mapsto\gamma W$, where $W\in U(N)$ commutes with $\Lambda$, satisfying
$W\Lambda W^{-1}=\Lambda$. This ambiguity, however, does not affect $Q$ and
$R_{\theta}$. To illustrate, consider the $N$ eigenvalues divided into groups
of $|I|$ degenerate ones, each with eigenvalue $\lambda_{I}$ so that
$\Lambda=\bigoplus_{I}\lambda_{I}{\bf 1}_{|I|}$. Introduce a block matrix
notation $A_{IJ}=(u_{i}^{\dagger}du_{j})_{i\in I,j\in J}$. The transformation
matrix $W$ is expressed as $W=\bigoplus_{I}W_{I}$ as well, with $W_{I}\in
U(|I|)$, modifying $A_{IJ}$ to
$A_{IJ}=W_{I}^{\dagger}A_{IJ}W_{J}+\delta_{IJ}W_{I}^{-1}dW_{I}$. Consequently,
$Q$ and $R_{\theta}$ can be represented as:

$\displaystyle Q$ $\displaystyle=\frac{1}{4\pi}\sum_{I,J;I\neq
J}\mathrm{Tr}_{I}[A_{IJ}A_{JI}]\lambda_{J}\lambda_{I}^{-1},$ (8)
$\displaystyle R_{\theta}$ $\displaystyle=\frac{1}{2\pi}\sum_{I,J;I\neq
J}\mathrm{Tr}_{I}[A_{IJ}A_{JI}]\log_{\theta}\lambda_{I},$ (9)

where $\mathrm{Tr}_{I}$ denotes the trace over indices $i\in I$. In the
summation $\sum_{I,J}$, terms with $I=J$ can be excluded due to
$\mathrm{Tr}_{I}[(A_{II})^{2}]=0$. This demonstrates the invariance of $Q$ and
$R_{\theta}$ under the transformation $\gamma\mapsto\gamma W$.

Another noteworthy aspect is that the difference in $B_{\theta}$ between two
$\theta$-gaps is a total derivative. For $0\leq\theta_{1},\theta_{2}<2\pi$,
and using (4), it follows that:

$\displaystyle B_{\theta_{1}}-B_{\theta_{2}}=d\alpha_{\theta_{1},\theta_{2}},$
(10)

where

$\displaystyle\alpha_{\theta_{1},\theta_{2}}$ $\displaystyle=-i\
\mathrm{sgn}(\theta_{1}-\theta_{2})\mathrm{Tr}[P_{\theta_{1},\theta_{2}}\gamma^{-1}d\gamma]$
$\displaystyle=-i\ \mathrm{sgn}(\theta_{1}-\theta_{2})$
$\displaystyle\quad\times\sum_{n;\min(\theta_{1},\theta_{2})<\arg\lambda_{n}<\max(\theta_{1},\theta_{2})}u_{n}^{\dagger}du_{n},$
(11)

with

$\displaystyle
P_{\theta_{1},\theta_{2}}=\sum_{n;\min(\theta_{1},\theta_{2})<\arg\lambda_{n}<\max(\theta_{1},\theta_{2})}u_{n}u_{n}^{\dagger},$
(12)

the orthogonal projection onto the eigenspace for eigenvalues that fulfill
$\min(\theta_{1},\theta_{2})<\arg\lambda_{n}<\max(\theta_{1},\theta_{2})$. In
cases where $\theta_{1}=\theta_{2}$ or no eigenvectors meet the condition
$\min(\theta_{1},\theta_{2})<\arg\lambda_{n}<\max(\theta_{1},\theta_{2})$,
$\alpha_{\theta_{1},\theta_{2}}$ is simply null.

Now, we express the winding number $W_{3}[g]$ as a sum of line integrals,
utilizing a cubic decomposition $L$ of the manifold $X$. (Any simplicial
decomposition is equally valid.) Within each cube $c$ of the lattice $L$, we
select $\theta_{c}\in[0,2\pi)$ such that $g_{x}$ for $x\in c$ exhibits a
$\theta$-gap of $\theta_{c}$. Thus, $W_{3}[g]$ can be reformulated as a sum of
integrals over all plaquettes:

$\displaystyle W_{3}[g]$
$\displaystyle=\frac{1}{2\pi}\sum_{c}\int_{c}dB_{\theta_{c}}$
$\displaystyle=\frac{1}{2\pi}\sum_{c}\int_{\partial c}B_{\theta_{c}}$
$\displaystyle=\frac{1}{2\pi}\sum_{p}\int_{p}(B_{\theta_{p}^{-}}-B_{\theta_{p}^{+}}).$
(13)

Here, $\sum_{p}$ runs over all plaquettes $p$ in the lattice $L$, with each
$p$ being oriented. The gap parameters $\theta_{p}^{+}$ and $\theta_{p}^{-}$
correspond to the cubes adjacent to plaquette $p$, in directions parallel and
antiparallel to $p$’s normal vector, respectively, as depicted in Fig. 2 (a).
This formulation further simplifies to a sum of line integrals:

$\displaystyle W_{3}[g]$
$\displaystyle=\frac{1}{2\pi}\sum_{p}\int_{p}d\alpha_{\theta_{p}^{-},\theta_{p}^{+}}$
$\displaystyle=\frac{1}{2\pi}\sum_{p}\oint_{\partial
p}\alpha_{\theta_{p}^{-},\theta_{p}^{+}}.$ (14)

The transition to the last expression is contingent upon
$\alpha_{\theta_{p}^{-},\theta_{p}^{+}}$ being smoothly defined across
plaquette $p$. When $\alpha_{\theta_{p}^{-},\theta_{p}^{+}}$ is solely
determined along the loop $\partial p$ bordering plaquette $p$, a $2\pi$
ambiguity may arise from large gauge transformations $u_{n}\to
u_{n}e^{i\chi_{n}}$, where $\oint_{\partial p}d\chi_{n}=2\pi$, potentially
altering $W_{3}[g]$ by an integer. However, if the cubic lattice $L$ is
sufficiently fine relative to $g$’s spatial variations, the integral
$\oint_{\partial p}\alpha_{\theta_{p}^{-},\theta_{p}^{+}}$ approximates 0 mod
$2\pi$, permitting its interpretation as an $\mathbb{R}$-valued quantity
devoid of the $2\pi$ ambiguity.

Figure 2:  (a) A plaquette $p$ within the cubic lattice, showing
$\theta_{p}^{+}$ and $\theta_{p}^{-}$ as the $\theta$-gaps of cubes adjacent
to $p$, aligned parallel and anti-parallel to $p$’s normal vector,
respectively. The vertices $v_{0},v_{1},v_{2},$ and $v_{3}$ are sequentially
labeled around the perimeter of plaquette $p$. (b) An edge $v_{a}v_{b}$ of the
cubic lattice, illustrating the $\theta$-gaps of cubes adjacent to the edge
$v_{a}v_{b}$.

We claim that the winding number as expressed in (14) can be calculated solely
using the diagonalizing matrices $\gamma(v)$ at the vertices $v$ of lattice
$L$. Diagonalizing $g(v)$ at vertices $v\in L$ yields the eigenvector and
eigenvalue pairs $\\{u_{n}(v),\lambda_{n}(v)\\}_{n=1,\dots,N}$ for each vertex
$v$. The ordering of eigenvectors $u_{n}(v)$ is such that the angles of
eigenvalues ascend, satisfying
$0\leq\lambda_{1}(v)\leq\cdots\leq\lambda_{N}(v)<2\pi$. (Note that eigenvalues
$\lambda$ near $0$ may reorder significantly under minor perturbations, yet
this does not contribute to the discrete formula below.) The gap parameter
$\theta_{c}$ for each cube $c$ is determined as follows: From the eight
vertices of cube $c$, denoted as $v\in c$, we derive $8N$ eigenvalues
$\\{\lambda_{n}(v)\\}_{v\in c,n=1,\dots,N}$. By smearing all eigenvalues
$\lambda_{n}(v)$, we have a set of intervals:

$\displaystyle I_{c}$ $\displaystyle=\bigcup_{v\in
c,n=1,\dots,N}\Bigg{\\{}\arg(\lambda_{n}(v)e^{i\delta\phi})\in[0,2\pi)\Bigg{|}$
$\displaystyle\hskip
30.0pt-\frac{\beta}{2N}<\delta\phi<\frac{\beta}{2N}\Bigg{\\}}.$ (15)

Here, $0<\beta<1$ is a constant smearing parameter ensuring that adjacent
eigenvalues fall within the same smeared interval. For example. we can set as
$\beta=1/2$. We select a $\theta$ from the set $[0,2\pi)\backslash I_{c}$ to
serve as $\theta_{c}$ for cube $c$. (Refer to Fig. 1 (c) for visualization.)
With $\theta$-gaps for all cubes in lattice $L$ thus defined, the gap
parameters $\theta_{p}^{+}$ and $\theta_{p}^{-}$ for each plaquette $p$ are
specified. For each corner vertex $v_{0},v_{1},v_{2},v_{3}$ of plaquette $p$,
we define $\gamma_{\\{\theta_{p}^{+},\theta_{p}^{-}\\}}(v_{j})$ as a $N\times
N_{q}$ matrix:

$\displaystyle\gamma_{\\{\theta_{p}^{+},\theta_{p}^{-}\\}}(v_{a})=\left(u_{n_{1}}(v_{a}),\dots,u_{n_{N_{q}}}(v_{a})\right),$
(16)

comprising $N_{q}$ eigenvectors, with eigenvalue angles satisfying
$\min(\theta_{p}^{+},\theta_{p}^{-})<\arg\lambda_{n_{1}}(v_{a})\leq\cdots\leq\arg\lambda_{n_{N_{q}}}(v_{a})<\max(\theta_{p}^{+},\theta_{p}^{-})$.
The nonnegative integer $N_{q}$, indicating the count of eigenvalues between
$e^{i\theta_{p}^{+}}$ and $e^{i\theta_{p}^{-}}$, should be in common for the
four vertices $v_{0},v_{1},v_{2},v_{3}$ of plaquette $p$, assuming the lattice
$L$ is fine enough. Then, the line integral $\oint_{\partial
p}\alpha_{\theta_{p}^{-},\theta_{p}^{+}}$ can be approximated as

$\displaystyle\oint_{\partial p}\alpha_{\theta_{p}^{-},\theta_{p}^{+}}$
$\displaystyle\cong\mathrm{sgn}(\theta_{p}^{+}-\theta_{p}^{-})\times\mathrm{Arg}\,\det\left[\gamma_{\\{\theta_{p}^{+},\theta_{p}^{-}\\}}(v_{0})^{\dagger}\gamma_{\\{\theta_{p}^{+},\theta_{p}^{-}\\}}(v_{3})\gamma_{\\{\theta_{p}^{+},\theta_{p}^{-}\\}}(v_{3})^{\dagger}\gamma_{\\{\theta_{p}^{+},\theta_{p}^{-}\\}}(v_{2})\right.$
$\displaystyle\hskip
50.0pt\left.\times\gamma_{\\{\theta_{p}^{+},\theta_{p}^{-}\\}}(v_{2})^{\dagger}\gamma_{\\{\theta_{p}^{+},\theta_{p}^{-}\\}}(v_{1})\gamma_{\\{\theta_{p}^{+},\theta_{p}^{-}\\}}(v_{1})^{\dagger}\gamma_{\\{\theta_{p}^{+},\theta_{p}^{-}\\}}(v_{0})\right]=:\Phi_{p},$
(17)

where $\mathrm{Arg}$ denotes the principal value ${-\pi<\mathrm{Arg}\ z<\pi}$.
Consequently, we obtain a discrete formula for the winding number:

$\displaystyle W^{\rm
dis}_{3}[g]=\frac{1}{2\pi}\sum_{p}\Phi_{p}\in\mathbb{Z},$ (18)

which relies solely on the diagonalization of matrix $g(v)$ at vertices $v\in
L$ of the cubic lattice $L$ approximating the manifold $X$.

The discrete formula (18) is inherently quantized. To demonstrate this, we
express $e^{2\pi iW^{\mathrm{dis}}_{3}[g]}$ as the product of edge
contributions:

$\displaystyle e^{2\pi iW^{\mathrm{dis}}_{3}[g]}=\prod_{p}e^{i\Phi_{p}}$
$\displaystyle=\prod_{v_{a}v_{b}\in\\{\mathrm{edges}\\}}\exp\Bigg{[}i\arg\Big{\\{}\det\left[\gamma_{\\{\theta_{1},\theta_{0}\\}}(v_{a})^{\dagger}\gamma_{\\{\theta_{1},\theta_{0}\\}}(v_{b})\right]^{\mathrm{sgn}(\theta_{1}-\theta_{0})}\det\left[\gamma_{\\{\theta_{2},\theta_{1}\\}}(v_{a})^{\dagger}\gamma_{\\{\theta_{2},\theta_{1}\\}}(v_{b})\right]^{\mathrm{sgn}(\theta_{2}-\theta_{1})}$
$\displaystyle\quad\times\det\left[\gamma_{\\{\theta_{3},\theta_{2}\\}}(v_{a})^{\dagger}\gamma_{\\{\theta_{3},\theta_{2}\\}}(v_{b})\right]^{\mathrm{sgn}(\theta_{3}-\theta_{2})}\det\left[\gamma_{\\{\theta_{0},\theta_{3}\\}}(v_{a})^{\dagger}\gamma_{\\{\theta_{0},\theta_{3}\\}}(v_{b})\right]^{\mathrm{sgn}(\theta_{0}-\theta_{3})}\Big{\\}}\Bigg{]}.$
(19)

Here, $v_{a}v_{b}$ denotes an individual edge, and
$\theta_{0},\dots,\theta_{3}$ are the gap parameters for cubes adjacent to the
edge $v_{a}v_{b}$, ordered counterclockwise from the vector
$\overrightarrow{v_{b}v_{a}}$ pointing out of the page. See Fig. 2 (b).
Regardless of the relative magnitudes of $\theta_{0},\theta_{1},\theta_{2}$,
and $\theta_{3}$, each edge’s contribution in (19) cancels out exactly due to
the property that for $\theta<\theta^{\prime}<\theta^{\prime\prime}$,

$\displaystyle\gamma_{\\{\theta,\theta^{\prime\prime}\\}}(v)=\left(\gamma_{\\{\theta,\theta^{\prime}\\}}(v_{a}),\gamma_{\\{\theta^{\prime},\theta^{\prime\prime}\\}}(v_{a})\right).$
(20)

Thus, $e^{2\pi iW^{\mathrm{dis}}_{3}[g]}=1$ is valid for any sufficiently fine
discrete approximation of $X$.

Model calculation— Our formulation extends to computing the winding number for
maps $g:X\to GL_{N}(\mathbb{C})$, where the target space consists of
invertible matrices. Invertible matrices that cannot be diagonalized, known as
”exceptional points,” constitute a ring in the three-dimensional parameter
space and are stable under minor perturbations. To circumvent these
exceptional points, one can employ the singular value decomposition $g=U\Sigma
V^{\dagger}$ to derive the unitary matrix $UV^{\dagger}$ at each vertex within
the discretized parameter space. We verified our formula (18) with the model
$g(k_{x},k_{y},k_{z})=t(\sin k_{x}\sigma_{x}+\sin k_{y}\sigma_{y}+\sin
k_{z}\sigma_{z})-i(m+\cos k_{x}+\cos k_{y}+\cos k_{z})\mathbf{1}_{2}$ on the
three-torus $(k_{x},k_{y},k_{z})\in[-\pi,\pi]^{\times 3}$, employing a cubic
lattice of $20\times 20\times 20$ mesh. Here,
$\sigma_{\mu}\in\\{\sigma_{x},\sigma_{y},\sigma_{z}\\}$ denotes the Pauli
matrices. We have checked that the winding number $W^{\rm dis}_{3}[g]$ equals
$-2\mathrm{sgn}(t)$ for $|m|<1$, $\mathrm{sgn}(t)$ for $1<|m|<3$, and $0$ for
$|m|>3$, which is consistent with the direct calculation of the analytic form
$W_{3}[g]$.

Summary— In this work, we presented a formulation for calculating the three-
dimensional winding number $W_{3}[g]$ for smooth maps $g:X\to U(N)$, utilizing
a discrete approximation of the manifold $X$. Our approach allows for the
computation of $W_{3}[g]$ exclusively through the diagonalization of matrices
$g(v)$ at a finite number of vertices, ensuring the result is explicitly
quantized to integer values. Discrete formulations that explicitly quantify
topological invariants are currently limited in scope. Examples such as
instanton numbers represented by the second Chern numbers, higher-dimensional
winding numbers, and degrees of maps to more general symmetric spaces have yet
to be explored. We look forward to future studies shedding more light on these
topics.

###### Acknowledgements.

We were supported by JST CREST Grant No. JPMJCR19T2, and JSPS KAKENHI Grant
No. 22H05118 and 23H01097.

## References

  * Volovik [2003] G. E. Volovik,  _The universe in a helium droplet_ , Vol. 117 (OUP Oxford, 2003).
  * Schnyder  _et al._ [2008] A. P. Schnyder, S. Ryu, A. Furusaki, and A. W. W. Ludwig, Classification of topological insulators and superconductors in three spatial dimensions, Phys. Rev. B 78, 195125 (2008).
  * Nakahara [2018] M. Nakahara,  _Geometry, topology and physics_ (CRC press, 2018).
  * Lüscher [1982] M. Lüscher, Topology of lattice gauge fields, Communications in Mathematical Physics 85, 39 (1982).
  * Fujiwara  _et al._ [2001] T. Fujiwara, H. Suzuki, and K. Wu, Topological Charge of Lattice Abelian Gauge Theory, Progress of Theoretical Physics 105, 789 (2001), https://academic.oup.com/ptp/article-pdf/105/5/789/5186102/105-5-789.pdf .
  * Fukui  _et al._ [2005] T. Fukui, Y. Hatsugai, and H. Suzuki, Chern numbers in discretized brillouin zone: Efficient method of computing (spin) hall conductances, Journal of the Physical Society of Japan 74, 1674 (2005), https://doi.org/10.1143/JPSJ.74.1674 .
  * Knöppel and Pinkall [2016] F. Knöppel and U. Pinkall, Complex line bundles over simplicial complexes and their applications, in  _Advances in Discrete Differential Geometry_, edited by A. I. Bobenko (Springer Berlin Heidelberg, Berlin, Heidelberg, 2016) pp. 197–239.
  * Carpentier  _et al._ [2015] D. Carpentier, P. Delplace, M. Fruchart, K. Gawedzki, and C. Tauber, Construction and properties of a topological index for periodically driven time-reversal invariant 2d crystals, Nuclear Physics B 896, 779 (2015).
  * Gawedzki and Reis [2002] K. Gawedzki and N. Reis, WZW branes and gerbes, Reviews in Mathematical Physics 14, 1281 (2002), https://doi.org/10.1142/S0129055X02001557 .


# SafeRL-Kit: Evaluating Efficient Reinforcement Learning Methods
for Safe Autonomous Driving

Linrui Zhang  Qin Zhang  Li Shen  Bo Yuan  Xueqian Wang

###### Abstract

Safe reinforcement learning (RL) has achieved significant success on risk-
sensitive tasks and shown promise in autonomous driving (AD) as well.
Considering the distinctiveness of this community, efficient and reproducible
baselines are still lacking for safe AD. In this paper, we release SafeRL-Kit
to benchmark safe RL methods for AD-oriented tasks. Concretely, SafeRL-Kit
contains several latest algorithms specific to zero-constraint-violation
tasks, including Safety Layer, Recovery RL, off-policy Lagrangian method, and
Feasible Actor-Critic. In addition to existing approaches, we propose a novel
first-order method named Exact Penalty Optimization (EPO) and sufficiently
demonstrate its capability in safe AD. All algorithms in SafeRL-Kit are
implemented (i) under the off-policy setting, which improves sample efficiency
and can better leverage past logs; (ii) with a unified learning framework,
providing off-the-shelf interfaces for researchers to incorporate their
domain-specific knowledge into fundamental safe RL methods. Conclusively, we
conduct a comparative evaluation of the above algorithms in SafeRL-Kit and
shed light on their efficacy for safe autonomous driving. The source code is
available at this https URL.

Machine Learning, ICML

##  1 Introduction

Figure 1: The overall framework of SafeRL-Kit. The trajectory collector
interacts with specified AD environments (e.g., MetaDrive (Li et al., 2021))
and stores transitions in the memory. SafeRL-Kit contains several safe RL
agents that efficiently learn from past experiences, including Safety Layer,
Recovery RL, Off-policy Lagrangian, Feasible Actor Critic, and newly proposed
Exact Penalty Optimization.

Reinforcement Learning (RL) has achieved superhuman performance in many
decision-making problems (Mnih et al., 2015; Vinyals et al., 2019). Typically,
the agent learns from trial and error and requires minimal prior knowledge of
the environment. Such a paradigm has natural advantages in mastering complex
skills for highly nonlinear systems like autonomous vehicles (Kiran et al.,
2021).

Nevertheless, concerns about the systematic safety limit the widespread use of
standard RL in real-world applications (Amodei et al., 2016). As an
alternative, safe RL takes safety requirements as hard constraints and
optimizes policies in the feasible domain. In recent years, it has been deemed
as a practical solution to resource allocation (Liu et al., 2021), robotic
locomotion (Yang et al., 2022), etc.

There have also been studies introducing safe RL into autonomous driving (AD)
(Isele et al., 2018; Chen et al., 2021; Li et al., 2022). Despite those
ongoing efforts, a unified benchmark is of great relevance to facilitate
further research on safe AD. We notice some risk-sensitive simulated
environments (Li et al., 2021; Herman et al., 2021) have been proposed, but an
efficient safe RL toolkit is still absent for this community. Considering the
distinctiveness of AD-oriented tasks, common code-bases (Ray et al., 2019;
Yuan et al., 2021) lack the following pivotal characteristics:

(1) Being safety-critical.

The agent must maintain zero cost-return as much as possible since any
inadmissible behavior in autopilot leads to catastrophic failures. Instead,
the previous code-base is built for a general-purpose with trajectory-based
constraints and non-zero thresholds.

(2) Being sample-efficient.

Off-policy algorithms can better leverage past logs and human demonstrations,
which is crucial for AD. By contrast, the previous code-base requires tens of
millions of interactions due to its on-policy algorithms, like CPO and PPO-L
(Ray et al., 2019).

(3) Being up-to-date.

There has been a fast-growing body of RL-based safe control. Nevertheless, the
previous code-base merely contains elder baselines (Achiam et al., 2017; Chow
et al., 2017) and lacks the latest advances.

(4) Being easy-to-use.

Most work on learning-based safe AD tends to incorporate domain-specific
knowledge into fundamental safe RL. Thus the toolkit is supposed to provide
off-the-shelf interfaces for extended studies. However, the modules of the
previous code-base are highly coupled and are implemented with the deprecated
TensorFlow version.

To provide a such as toolkit for safe RL algorithms and understand which of
them are best suited for AD-oriented tasks, our contributions in this work are
summarized as the following three-folds:

  * •

We release SafeRL-Kit, which contains the latest advances in safe RL (Dalal et
al., 2018; Ha et al., 2020; Thananjeyan et al., 2021; Ma et al., 2021). All
algorithms are implemented efficiently under off-policy settings and with a
unified training framework.

  * •

We propose a novel first-order method coined Exact Penalty Optimization (EPO)
and incorporate it into SafeRL-Kit. EPO utilizes a single penalty factor and a
ReLU operator to construct an equivalent unconstrained objective. Empirical
results show the simple technique is surprisingly effective and robust for AD-
oriented tasks.

  * •

We benchmark SafeRL-Kit in a representative toy environment and a simulated
platform with realistic vehicle dynamics. To the best of our knowledge, this
paper is the first to provide unified off-policy safe RL baselines and a fair
comparison of them specific to AD.

##  2 Related Work

###  2.1 Safe RL Algorithms

A number of works tackle RL-based safe control for autonomous agents, and we
divide them into three genres. The first type of method, coined as safe policy
optimization, incorporates safety constraints into the standard RL objective
and yields a constrained sequential optimization problem (Chow et al., 2017;
Achiam et al., 2017; Zhang et al., 2020; Ma et al., 2021; Zhang et al., 2022).
The second type of method, coined as safety correction, projects initial
unsafe behaviors to the feasible region (Dalal et al., 2018; Zhao et al.,
2021). The third type of method, coined as safety recovery, learns an
additional pair of safe actor-critic to take over control when encountering
potential risks (Thananjeyan et al., 2021; Yang et al., 2022).

There have also been studies on safe RL specific to AD-oriented tasks. Isele
et al. (2018) utilizes a prediction module to generate masks on dangerous
behaviors, which merely works in discrete action spaces. Wen et al. (2020)
extend Constrained Policy Optimization (CPO) (Achiam et al., 2017) to AD and
employ synchronized parallel actors to accelerate the convergence speed for
on-policy CPO. Chen et al. (2021) take the ego-camera view as input and train
an additional recovery policy via a heuristic objective based on Hamilton-
Jacobi reachability. Li et al. (2022) propose a human-in-loop approach to
learn safe driving efficiently.

###  2.2 Safe RL Benchmarks

For general scenarios, a set of benchmarks are commonly used to evaluate the
efficacy of safe RL algorithms. The classic
environments111https://github.com/SvenGronauer/Bullet-Safety-Gym include Robot
with Limit Speed (Zhang et al., 2020), Circle and Gather (Achiam et al.,
2017), etc. Safety-gym222https://github.com/openai/safety-gym (Ray et al.,
2019) contains several tasks (goal, button, push) and agents (point, car,
doggo) that are representative in robot control problems. Meanwhile, the
authors provide popular baselines333https://github.com/openai/safety-starter-
agents, including CPO and some on-policy Lagrangian methods. Safe-control-
gym444https://github.com/utiasDSL/safe-control-gym (Yuan et al., 2021) bridges
the gap between control and RL communities. The authors also developed an
open-sourced toolkit supporting both model-based and data-driven control
techniques.

For AD-oriented tasks, there have been some existing environments for safe
driving. Li et al. (2021) release
Metadrive555https://github.com/metadriverse/metadrive that benchmarks
reinforcement learning algorithms for vehicle autonomy, including safe
exploitation and exploration. Herman et al. (2021) propose Learn-to-
Race666https://github.com/learn-to-race/l2r that focuses on safe control in
high speed. Nevertheless, it still lacks a set of strong baselines specific to
the AD community considering its distinctiveness depicted above in Section 1.
To our best knowledge, this paper is the first to provide unified off-policy
safe RL baselines and a fair comparison of them for the purpose of autonomous
driving.

##  3 Preliminaries

(a) Cost Signal = 0 (b) Cost Signal = 1

Figure 2: SpeedLimit Benchmark. The vehicle is rewarded for driving along the
avenue, but receives a cost signal if $vel>1.5m/s$.

A Markov Decision Process (MDP) (Sutton & Barto, 1998) is defined by a tuple
$(\mathcal{S},\mathcal{A},\mathcal{P},\mathcal{R},\mu,\gamma)$.
${\mathcal{S}}$ and ${\mathcal{A}}$ denote the state space and the action
space respectively.
${\mathcal{P}}:{\mathcal{S}}\times\mathcal{A}\times\mathcal{S}\mapsto[0,1]$ is
the transition probability function to describe the dynamics of the system.
$\mathcal{R}:\mathcal{S}\times\mathcal{A}\mapsto\mathbb{R}$ is the reward
function. $\mu:\mathcal{S}\mapsto[0,1]$ is the initial state distribution.
$\gamma$ is the discount factor for future reward. A stationary policy
$\pi:S\mapsto P(A)$ maps the given states to probability distributions over
action space. The goal of standard RL is to find the optimal policy $\pi^{*}$
that maximizes the expected discounted return
$J_{R}(\pi)=\mathop{\mathbb{E}}_{\tau\sim\pi}\big{[}\sum^{\infty}_{t=0}\gamma^{t}R(s_{t},a_{t})\big{]},$
where $\tau=\\{(s_{t},a_{t})\\}_{t\geq 0}$ is a sample trajectory and
$\tau\sim\pi$ accounts for the distribution over trajectories depending on
$s_{0}\sim\mu,a_{t}\sim\pi(\cdot|s_{t}),s_{t+1}\sim P(\cdot|s_{t},a_{t})$.

A Constrained Markov Decision Process (CMDP) (Altman, 1999) extends MDP to
$(\mathcal{S},\mathcal{A},\mathcal{P},\mathcal{R},\mathcal{C},\mu,\gamma)$.
The cost function $\mathcal{C}:\mathcal{S}\times\mathcal{A}\mapsto[0,+\infty]$
reflects the violation of systematic safety. The goal of safe RL is to find

$\pi^{*}={\arg\max}_{\pi}J_{R}(\pi)\quad\mathrm{s.t.}\ \ \\{a_{t}\\}_{t\geq
0}\text{ is feasible}.$

In a CMDP, the cost function is typically constrained in the following two
ways. The first is _Instantaneous Constrained MDP_. This type of Safe RL
formualtion requires the selected actions enforce the constraint at every
decision-making step, namely $C(s_{t},a_{t})\leq\epsilon$. The second is
_Cumulative Constrained MDP_. This type of Safe RL formualtion requires the
discounted sum of cost signals in the whole trajectory within a certain
threshold, namely
$J_{C}(\pi)=\mathop{\mathbb{E}}_{\tau\sim\pi}\big{[}\sum^{\infty}_{t=0}\gamma^{t}C(s_{t},a_{t})\big{]}\leq
d.$

##  4 Problem Setup

In this paper, we develop SafeRL-Kit to evaluate efficient RL algorithms for
safe autonomous driving on existing benchmarks. We simplify the cost function
as the following risk-indicator:

$C(s,a)=\begin{cases}1,&\text{if the transition is unsafe}\\\
0,&\text{otherwise}\end{cases}.$ (1)

This formulation is generalizable to different AD-oriented tasks without
cumbersome reward and cost shaping. The goal of the autonomous vehicle is to
reach the destination as fast as possible while adhering to zero cost signals
at every time steps. Specifically, we conduct comparative evaluations on a
representative toy environment and a simulated platform with realistic vehicle
dynamics respectively.

(a) Cost Signal = 0 (b) Cost Signal = 1

Figure 3: MetaDrive Benchmark. The vehicle aims to reach virtual markers, but
receives a cost signal if it collides with obstacles and other vehicles or it
is out of the road. Table 1: Comparison of different safe reinforcement
learning algorithms for AD-oriented tasks.

Algorithm | Constraint Type | Policy Type
---|---|---
Cumulative/Instantaneous | State-wise/Trajectory-wise | Deterministic | Stochastic
CPO (Ray et al., 2019) | Cumulative | Trajectory-wise | $\times$ | $\surd$
PPO-L (Ray et al., 2019) | Cumulative | Trajectory-wise | $\times$ | $\surd$
TRPO-L (Ray et al., 2019) | Cumulative | Trajectory-wise | $\times$ | $\surd$
Safety Layer | Instantaneous | State-wise | $\surd$ | $\times$
Recovery RL | Cumulative | State-wise | $\surd$ | $\surd$
Off-policy Lagrangian | Cumulative | Trajectory-wise | $\surd$ | $\surd$
Feasible Actor-Critic | Cumulative | State-wise | $\surd$ | $\surd$
Exact Penalty Optimization | Cumulative | Both | $\surd$ | $\surd$

###  4.1 SpeedLimit Benchmark

The task is inspired by Zhang et al. (2020), as illustrated in Figure 2. In
SpeedLimit task, the agent is a four-wheeled race-car whose observation is ego
position, velocity and yaw. The selected action controls the Revolution Per
Minute (RPM) and steering of wheels. The agent is rewarded for approaching
$x_{dest}=+\infty$ and the cost function is

$C(s,a)=\begin{cases}1,&\text{if vehicle's velocity}>1.5m/s\\\
0,&\text{otherwise}\end{cases}.$ (2)

The toy environment is simple yet representative since speed control is a
classic problem in vehicle autonomy. Besides, the speed limit is easy to reach
and thus undesirable algorithms may violate the safety constraint at almost
every time step. That is, the toy environment enables us to see which
algorithms can effectively degrade the dense cost return and are best suited
for safe AD tasks.

###  4.2 MetaDrive Benchmark

This task is inspired by Li et al. (2021), as illustrated in Figure 3.
Metadrive is a compositional, lightweight and realistic platform for vehicle
autonomy. Most importantly, it provides pre-defined environments for safe
policy learning in autopilots. Concretely, the observation is encoded by a
vector containing ego-state, navigation information and surrounding
information detected by the Lidar. We control the speed and steering of the
car to hit virtual land markers for rewards, and the cost function is defined
as

$C(s,a)=\begin{cases}1,&\text{if collides or out of the road}\\\
0,&\text{otherwise}\end{cases}$ (3)

It worth mentioning that we set the traffic density twice than the original
paper to construct a more challenging scenario.

##  5 Efficient Safe RL Algorithms

###  5.1 Overall Implementation

The current version of SafeRL-Kit contains some latest RL-based methods,
including _Safety Layer_ (Dalal et al., 2018), _Recovery RL_ (Thananjeyan et
al., 2021), _Off-policy Lagrangian_ (Ha et al., 2020), _Feasible Actor-Critic_
(Ma et al., 2021) and newly proposed _Exact Penalty Optimization_. We compare
above methods along with some existing on-policy baselines (Ray et al., 2019)
in Table 1.

Before diving into algorithmic details, we first explain the overall
implementation of SafeRL-Kit and its benefits:

(1)

The adopted algorithms address safe policy learning from different
perspectives (Safety Layer for safety correction; Recovery RL for safety
recovery; Lagrangian, FAC, and EPO for constrained optimization). Thus, users
can combine AD-specific knowledge with the proper type of safe RL baselines in
their studies.

(2)

All the algorithms are implemented under the off-policy Actor-Critic
architecture. Thus, they enjoy better sample efficiency and can leverage human
demonstration if needed.

(3)

All the algorithms are implemented with a unified training framework. By
default, all networks are MLPs with (256,256) hidden layers activated via the
ReLU function. The essential updates of backbone networks follow TD3 (Fujimoto
et al., 2018) without pre-training processes. Thus, we can conduct a fair
comparison to see which of them are best suited for AD-oriented tasks.

###  5.2 Safety Layer

Safety Layer, added on top of the original policy network, conducts a
quadratic-programming-based constrained optimization to find the ”nearest”
action to the feasible region.

Specifically, Safety Layer utilizes a parametric linear model

$C(s_{t},a_{t})\approx g(s_{t};\omega)^{\top}a_{t}+c_{t-1}$ (4)

to approximate the single-step cost function with supervised training and
yields the following QP problem

$\displaystyle a_{t}^{*}=$ $\displaystyle\ \ {\arg\min}_{a}\
\frac{1}{2}||a-\mu_{\theta}(s)||^{2}$ (5) $\displaystyle\mathrm{s.t.}\quad
g(s_{t};\omega)^{\top}a_{t}+c_{t-1}\leq\epsilon,$

which projects the unsafe action back to the feasible region.

Since there is only one compositional cost signal in our problem, the closed-
form solution of problem (5) is

$a_{t}^{*}=\mu_{\theta}(s_{t})-\bigg{[}\frac{g(s_{t};\omega)^{\top}\mu_{\theta}(s)+c_{t-1}-\epsilon}{g(s_{t};\omega)^{\top}g(s_{t};\omega)}\bigg{]}^{+}g(s_{t};\omega)$
(6)

Thus, Safety Layer is the type of method that addresses state-wise,
instantaneous constraints.

By the way, the $g_{\omega}$ is trained from offline data in Dalal et al.
(2018). SafeRL-Kit instead learns the linear model with the policy network
synchronously, considering the side-effect of distribution shift. We employ a
warm-up in the training process to avoid meaningless, inaccurate corrections.

Table 2: Hyper-parameters of different safety-aware algorithms in SafeRL-Kit.

Hyper-parameter |  Safety Layer |  Recovery RL |  Lagrangian |  FAC |  EPO
---|---|---|---|---|---
Cost Limit |  0.02 |  0.1 |  0.1 |  0.1 |  0.1
Reward Discount |  0.99 |  0.99 |  0.99 |  0.99 |  0.99
Cost Discount |  0.99 |  0.99 |  0.99 |  0.99 |  0.99
Warm-up Ratio |  0.2 |  0.2 |  N/A |  N/A |  N/A
Batch Size |  256 |  256 |  256 |  256 |  256
Critic LR |  3E-4 |  3E-4 |  3E-4 |  3E-4 |  3E-4
Actor LR |  3E-4 |  3E-4 |  3E-4 |  3E-4 |  3E-4
Safe Critic LR |  3E-4 |  3E-4 |  3E-4 |  3E-4 |  3E-4
Safe Actor LR |  N/A |  3E-4 |  N/A |  N/A |  N/A
Multiplier LR |  N/A |  N/A |  1E-5 |  1E-5 |  N/A
Multiplier Init |  N/A |  N/A |  0.0 |  N/A |  N/A
Policy Delay |  2 |  2 |  2 |  2 |  2
Multiplier Delay |  N/A |  N/A |  N/A |  12 |  N/A
Penalty Factor |  N/A |  N/A |  N/A |  N/A |  5

###  5.3 Recovery RL

The critical insight behind Recovery RL is to introduce an additional policy
that recovers potential unsafe states. Consequently, it trains two independent
RL agents instead of solving a cumbersome constrained optimization problem.

Specifically, Recovery RL learns a safe critic to estimate the future
probability of constraint violation as

$Q^{\pi}_{\text{risk}}(s_{t},a_{t})=c_{t}+(1-c_{t})\gamma\mathbb{E}_{\pi}Q^{\pi}_{\text{risk}}(s_{t+1},a_{t+1}).$
(7)

This formulation is slightly different from the standard Bellman equation
since it assumes the episode terminates when the agent receives a cost signal.
We found in experiments that such an early stopping makes it intractable for
agents to master desirable skills in AD. Thus, we remove the early-stopping
condition but still preserve the original formulation of
$Q^{\pi}_{\text{risk}}$ in (7) since it limits the upper bound of the safe
critic and eliminates the over-estimation in Q-learning.

In the phase of policy execution, the recovery policy takes over the control
when the predicted value of the safe critic exceeds the given threshold, as

$a_{t}=\begin{cases}\pi_{\text{task}}(s_{t}),&\text{if
}Q^{\pi}_{\text{risk}}\big{(}s_{t},\pi_{\text{task}}(s_{t})\big{)}\leq\epsilon\\\
\pi_{\text{risk}}(s_{t}),&\text{otherwise}\end{cases}$ (8)

The optimization objective of $\pi_{\text{task}}$ is to maximize the
cumulative rewards, whereas the goal of $\pi_{\text{risk}}$ is to minimize
$Q^{\pi}_{\text{risk}}$, namely to degrade the potential risk of the agent.

It is important to store $a_{\text{task}}$ and $a_{\text{risk}}$
simultaneously in the replay buffer, and utilize them to train
$\pi_{\text{task}}$ and $\pi_{\text{risk}}$ respectively in Recovery RL. This
technique ensures that $\pi_{\text{task}}$ can learn from the new MDP, instead
of proposing same unsafe actions continuously.

Similar to Safety Layer, Recovery RL in SafeRL-Kit also has a warm-up stage
where $Q^{\pi}_{\text{risk}}$ is trained but is not utilized.

###  5.4 Off-policy Lagrangian

Lagrangian Relaxation is commonly-used to address constrained optimization
problem. Safe RL as well can be formulated as a constrained sequential
optimization problem

$\displaystyle\mathop{\max}_{\pi}$
$\displaystyle\mathop{\mathbb{E}}_{s\sim\mu}V_{0}^{\pi}(s)$ (9)
$\displaystyle\mathrm{s.t.}$
$\displaystyle\mathop{\mathbb{E}}_{s\sim\mu}U^{\pi}_{0}(s)\leq\epsilon,$

where
$V^{\pi}_{0}(s)=\mathop{\mathbb{E}}_{\tau\sim\pi}\big{[}\sum^{\infty}_{t=0}\gamma^{t}r_{t}\big{|}s_{0}=s]$
and
$U^{\pi}_{0}(s)=\mathop{\mathbb{E}}_{\tau\sim\pi}\big{[}\sum^{\infty}_{t=0}\gamma^{t}c_{t}\big{|}s_{0}=s]$.

Strong duality holds for primal problem (9) (Paternain et al., 2022), thus it
can be tackled via the dual problem

$\mathop{\max}_{\lambda\geq
0}\mathop{\min}_{\pi}\mathop{\mathbb{E}}_{s\sim\mu}-V_{0}^{\pi}(s)+\lambda\big{(}U^{\pi}_{0}(s)-\epsilon\big{)}.$
(10)

The off-policy objective of problem (10) in the parametric space (Ha et al.,
2020) can be formulated as

$\mathop{\max}_{\lambda\geq
0}\mathop{\min}_{\theta}\mathbb{E}_{\mathcal{D}}-Q^{\pi}(s,\pi_{\theta}(s))+\lambda\big{(}Q^{\pi}_{c}(s,\pi_{\theta}(s))-\epsilon\big{)}.$
(11)

Stochastic primal-dual optimization (Luenberger et al., 1984) is applied here
to update primal and dual variables alternatively, which follows as

$\begin{cases}\theta\leftarrow\theta-\eta_{\theta}\nabla_{\theta}\mathbb{E}_{\mathcal{D}}\big{(}-Q^{\pi}(s,\pi_{\theta}(s))+\lambda
Q^{\pi}_{c}(s,\pi_{\theta}(s))\big{)}\\\
\lambda\leftarrow\big{[}\lambda+\eta_{\lambda}\mathbb{E}_{\mathcal{D}}\big{(}Q^{\pi}_{c}(s,\pi_{\theta}(s))-\epsilon\big{)}\big{]}^{+}\end{cases}$
(12)

Notably, the timescale of primal variable updates is required to be faster
than the timescale of Lagrange multipliers. Thus, we set
$\eta_{\theta}\gg\eta_{\lambda}$ in SafeRL-Kit.

###  5.5 Feasible Actor-Critic

The constraint of Off-policy Lagrangian in Section 5.4 is based on the
expectation of whole trajectories, which inevitably allows some unsafe roll-
outs. Ma et al. (2021) introduce a new concept, namely state-wise constraints
for cumulative cost-return which follows as

$\displaystyle\mathop{\max}_{\pi}$
$\displaystyle\mathop{\mathbb{E}}_{s\sim\mu}V_{0}^{\pi}(s)$ (13)
$\displaystyle\mathrm{s.t.}$ $\displaystyle U^{\pi}_{0}(s)\leq\epsilon,\forall
s\in\mathcal{I_{F}}.$

Here $s\in\mathcal{I_{F}}$ stands for all ”feasible” initial states. Also,
their theoretical results show that problem (13) is a stricter version than
problem (9), which provides strong safety guarantees for state-wise safe
control.

The dual problem of (13) is derived by rescaling the state-wise constraints
and follows as

$\mathop{\max}_{\lambda\geq
0}\mathop{\min}_{\pi}\mathop{\mathbb{E}}_{s\sim\mu}-V_{0}^{\pi}(s)+\lambda(s)\big{(}U^{\pi}_{0}(s)-\epsilon\big{)}.$
(14)

The distinctiveness of problem (14) is there are infinitely many Lagrangian
multipliers that are state-dependent. In SafeRL-Kit, we employ an neural
network $\lambda(s;\xi)$ activated by _Softplus_ function to map the given
state $s$ to its corresponding Lagrangian multiplier $\lambda(s)$.

The primal-dual ascents of policy network is similar to (12); the updates of
multiplier network is given by

$\xi\leftarrow\xi+\eta_{\xi}\nabla_{\xi}\mathbb{E}_{\mathcal{D}}\lambda(s;\xi)\big{(}Q^{\pi}_{c}(s,\pi_{\theta}(s))-\epsilon\big{)}.$
(15)

Besides, SafeRL-Kit also sets a different interval schedule $m_{\pi}$ (for
$\pi_{\theta}$ delay steps) and $m_{\lambda}$ (for $\lambda_{\xi}$ delay
steps) to stabilize the training process (Ma et al., 2021).

###  5.6 Exact Penalty Optimization

Algorithm 1 State-wise Exact Penalty Optimization

0: deterministic policy network $\pi(s;\theta)$; critic networks
$\hat{Q}(s,a;\phi)$ and $\hat{Q}_{c}(s,a;\varphi)$

1: for t in $1,2,...$ do

2: $a_{t}=\pi(s_{t};\theta)+\epsilon,\ \ \epsilon\sim\mathcal{N}(0,\sigma)$.

3: Apply $a_{t}$ to the environment.

4: Store the transition $(s_{t},a_{t},s_{t+1},r_{t},c_{t},d_{t})$ in
$\mathcal{B}$.

5: Sample a mini-batch of $N$ transitions from $\mathcal{B}$.

6:
$\varphi\leftarrow{\arg\min}_{\varphi}\mathop{\mathbb{E}}_{\mathcal{B}}\big{[}\hat{Q}_{c}(s,a;\varphi)-\big{(}c+\gamma_{c}(1-d)\hat{Q}_{C}(s^{\prime},\pi(s^{\prime};\theta);\varphi)\big{)}\big{]}^{2}$.

7:
$\phi\leftarrow{\arg\min}_{\phi}\mathop{\mathbb{E}}_{\mathcal{B}}\big{[}\hat{Q}(s,a;\phi)-\big{(}r+\gamma(1-d)\hat{Q}(s^{\prime},\pi(s^{\prime};\theta);\phi)\big{)}\big{]}^{2}$.

8:
$\theta\leftarrow{\arg\min}_{\theta}\mathop{\mathbb{E}}_{\mathcal{B}}\big{[}-\hat{Q}(s,\pi(s;\theta);\phi)+\kappa\cdot\max\\{0,\hat{Q}_{c}(s,\pi(s;\theta);\varphi)-\delta\\}\big{]}$.

9: end for

In this paper, we propose a simple-yet-effective approach motivated by the
exact penalty method.

######  Theorem 5.1.

Considering the following two problems

$\displaystyle\min f(x)\ \ \mathrm{s.t.}\ g_{i}(x)\leq 0,i=1,2,...$ (P)
$\displaystyle\min f(x)+\kappa\cdot\sum_{i}\max\\{0,g_{i}(x)\\}$ (Q)

Suppose $\lambda^{*}$ is the optimal Lagrange multiplier vector of problem
(P). If the penalty factor $\kappa\geq||\lambda^{*}||_{\infty}$, problem (P)
and problem (Q) share the same optimal solution set.

###### Proof.

See our recent work (Zhang et al., 2022). ∎

The above theorem enables us to construct an equivalent function whose
unconstrained minimizing points also solve the previous constrained problem.
Meanwhile, the unconstrained problem can tackle multiple constraints with
exactly one consistent penalty factor.

Thus, we simplify Lagrangian-based methods (i.e., Off-policy Lagrangian and
FAC) with this technique, considering that the single-constrained optimization
problem (9) and the multi-constrained optimization problem (13) are suited for
exact penalty method in Theorem 5.1. In this way, we can employ a single
minimization on primal variables with fixed penalty terms instead of
cumbersome min-max optimization over both primal and dual variables.

Below we merely summarize the state-wise Exact Penalty Optimization (EPO) in
Algorithm 1 as an alternative to FAC, since FAC provides stricter safety
guarantees but suffers from the oscillation and instability of the multiplier
network. The off-policy surrogate objective of state-wise EPO follows as

$\ell(\theta)=\mathbb{E}_{\mathcal{D}}-Q^{\pi}(s,\pi_{\theta}(s))+\kappa\big{[}Q^{\pi}_{c}(s,\pi_{\theta}(s))-\epsilon\big{]}^{+},$
(16)

where $\kappa$ is a fixed, sufficiently large hyper-parameter.

##  6 Empirical Analysis

We benchmark RL-based algorithms on SpeedLimit task (Zhang et al., 2020) and
MetaDrive platform (Li et al., 2021). Below, we give a comparative evaluation
according to the empirical results.

Table 3: Mean performance with normal 95% confidence for safety-aware
algorithms on benchmarks.

Environments | Safety Layer | Recovery RL | Lagrangian | FAC | EPO
---|---|---|---|---|---
SpeedLimit | Ep-Reward | $651.59\pm 10.70$ | $623.67\pm 99.58$ | $565.50\pm 69.29$ | $631.55\pm 34.92$ | $\bm{684.86\pm 3.19}$
Ep-Cost | $76.30\pm 9.07$ | $187.14\pm 96.50$ | $7.28\pm 3.11$ | $7.83\pm 5.23$ | $5.44\pm 0.53$
CostRate | $0.33\pm 0.01$ | $0.43\pm 0.06$ | $0.06\pm 0.01$ | $0.07\pm 0.01$ | $0.02\pm 0.01$
MetaDrive | SuccessRate | $0.73\pm 0.05$ | $0.78\pm 0.06$ | $0.74\pm 0.05$ | $0.68\pm 0.04$ | $0.73\pm 0.05$
Ep-Cost | $12.91\pm 1.10$ | $14.18\pm 1.92$ | $9.23\pm 4.88$ | $3.29\pm 0.50$ | $4.29\pm 0.71$
CostRate | $0.04\pm 0.001$ | $0.05\pm 0.001$ | $0.02\pm 0.01$ | $0.01\pm 0.01$ | $0.01\pm 0.01$

#### Unconstrained Reference.

We utilize TD3 (Fujimoto et al., 2018) as the unconstrained reference for
upper bounds of reward performance and constraint violations. For the
SpeedLimit task (500 max_episode_horizon), TD3 exceeds the velocity threshold
at almost every step with a near 100% cost rate. For the MetaDrive environment
(1000 max_episode_horizon), the agent receives sparse cost signals when it
collides with obstacles or is out of the road. Besides, the cost signals are
encoded into the reward function; otherwise, it would be too hard to learn
desirable behaviors (Li et al., 2021). Consequently, TD3 with reward-shaping
(TD3-RS) would not have that high cumulative costs as it does in the toy
environment.

#### Overall Performance.

The mean performances are summarized in Table 2 and the learning curves over
five seeds are shown in Figure 4 and 5. We conclude that Safety Layer and
Recovery RL are less effective in degrading cost return. They still have
around 10% safety violations in SpeedLimit, and the safety improvement in
MetaDrive is also limited. As for Safety Layer, the main reasons are that the
linear approximation to the cost function brings about non-negligible errors,
and the single-step correction is myopic for future risks. As for Recovery RL,
the estimation error of $Q_{\text{risk}}$ is probably the biggest problem
affecting the recovery effects. By contrast, Off-policy Lagrangian and FAC
have significantly lower cumulative costs. However, the Lagrangian-based
methods have the inherent problems from primal-dual ascents. For one thing,
the Lagrangian multiplier tuning causes oscillations of learning curves. For
another thing, those algorithms are susceptible to Lagrangian multipliers’
initialization and learning rate. We conclude that constrained optimization
still outperforms safety correction and recovery if the hyper-parameters are
appropriately settled. At last, we find that the newly proposed EPO is
surprisingly effective for learning safe AD. In SpeedLimit, it converges to a
high plateau quickly while adhering to an almost zero cost return. In
MetaDrive, it is still competitive with SOTA baselines. We regard the
underlying reason as that EPO is an equivalent form to FAC but reduces state-
dependent Lagrangian multipliers to one fixed hyper-parameter. The consistent
loss function stabilizes the training process compared with primal-dual
optimization.

#### Sensitivity Analysis.

In this paper, we study the sensitivity to hyper-parameters of Lagrangian-
based methods and EPO in Figure 6 and Figure 7 respectively. We found that
Lagrangian-based methods are susceptible to the learning rate of the
Lagrangian multiplier(s) in stochastic primal-dual optimization. First, the
oscillating $\lambda$ causes non-negligible deviations in the learning curves.
Besides, the increasing $\eta_{\lambda}$ may degrade the performance
dramatically. The phenomenon is especially pronounced in FAC, which has a
multiplier network to predict the state-dependent $\lambda(s;\xi)$. Thus, we
suggest $\eta_{\lambda}\ll\eta_{\theta}$ in practice. As for EPO, we find if
the penalty factor $\kappa$ is too small, the cost return may fail to
converge. Nevertheless, if $\kappa$ is sufficiently large, the learning curves
are robust and almost identical. Thus, we suggest $\kappa>5$ in experiments
and a grid search for better performance.

#### Sample Complexity.

Considering the difficulty of the above two tasks, we run $5\times 10^{5}$ and
$1\times 10^{6}$ interactive steps respectively to obtain admissible results.
Notably, previous on-policy codebases require significantly more samples for
convergence; for example, Ray et al. (2019) run $1\times 10^{7}$ interactive
steps even for toy environments. Thus, SafeRL-Kit with off-policy
implementations is much more sample-efficient compared to theirs, emphasizing
the applicability of SafeRL-Kit to data-expensive AD-oriented tasks.

##  7 Further Discussion

The released SafeRL-kit contains several SOTA off-policy safe RL methods that
are suited for safety-critical autonomous driving. We conduct the comparative
evaluation of those baselines over one representative toy environment and one
simulated AD platform, respectively. The proposed Exact Penalty Optimization
in this paper is easy-to-implement and surprisingly effective on AD-oriented
tasks. We think future work on SafeRL-kit from two aspects:

  * •

The off-policy implementation of SafeRL-Kit can naturally leverage offline
data, including past logs and human demonstrations, which are commonly used
and highly effective for AD-oriented tasks.

  * •

We only benchmark safe RL methods with vector input (ego-state, navigation
information, Lidar signals, etc.) in this paper. Nevertheless, vision-based AD
is still less studied in the current version of SafeRL-Kit.

(a) Eval Episode Reward (b) Eval Episode Cost

(c) Training Cost rate

Figure 4: Learning curves on the SpeedLimit benchmark. The x-axis is the
number of interactions with the simulator (500,000 total).

(a) Eval Success Rate (b) Eval Episode Cost

(c) Training Cost rate

Figure 5: Learning curves on the MetaDrive Benchmark. The x-axis is the number
of interactions with the simulator (1,000,000 total).

(a) Reward-Lag-SpeedLimit (b) Cost-Lag-SpeedLimit

(c) Reward-FAC-MetaDrive (d) Cost-FAC-MetaDrive

Figure 6: Sensitivity studies of Lagrangian-based methods. The first two
figures are reward and cost plots of Off-policy Lagrangian on SpeedLimit task
with different $\lambda$ learning rates. The last two figures are success rate
and cost plots of Feasible Actor-Critic on MetaDrive benchmark with different
$\lambda(s;\xi)$ learning rates.

(a) Eval Episode Reward (b) Eval Episode Cost

(c) Training Cost rate (d) Training Cost rate

Figure 7: Sensitivity studies of Exact Penalty Optimization. The first two
figures are reward and cost plots of EPO on the SpeedLimit task with different
penalty factors $\kappa$. The last two figures are the success rate and cost
plots of EPO on the MetaDrive benchmark with different penalty factors
$\kappa$.

## References

  * Achiam et al. (2017) Achiam, J., Held, D., Tamar, A., and Abbeel, P.  Constrained policy optimization.  In _International Conference on Machine Learning_ , pp. 22–31. PMLR, 2017.
  * Altman (1999) Altman, E.  _Constrained Markov decision processes_ , volume 7.  CRC Press, 1999.
  * Amodei et al. (2016) Amodei, D., Olah, C., Steinhardt, J., Christiano, P., Schulman, J., and Mané, D.  Concrete problems in ai safety.  _arXiv preprint arXiv:1606.06565_ , 2016.
  * Chen et al. (2021) Chen, B., Francis, J., Nyberg, J. O. E., and Herbert, S. L.  Safe autonomous racing via approximate reachability on ego-vision.  _arXiv preprint arXiv:2110.07699_ , 2021.
  * Chow et al. (2017) Chow, Y., Ghavamzadeh, M., Janson, L., and Pavone, M.  Risk-constrained reinforcement learning with percentile risk criteria.  _The Journal of Machine Learning Research_ , 18(1):6070–6120, 2017.
  * Dalal et al. (2018) Dalal, G., Dvijotham, K., Vecerik, M., Hester, T., Paduraru, C., and Tassa, Y.  Safe exploration in continuous action spaces.  _arXiv preprint arXiv:1801.08757_ , 2018.
  * Fujimoto et al. (2018) Fujimoto, S., Hoof, H., and Meger, D.  Addressing function approximation error in actor-critic methods.  In _International conference on machine learning_ , pp. 1587–1596. PMLR, 2018.
  * Ha et al. (2020) Ha, S., Xu, P., Tan, Z., Levine, S., and Tan, J.  Learning to walk in the real world with minimal human effort.  _arXiv preprint arXiv:2002.08550_ , 2020.
  * Herman et al. (2021) Herman, J., Francis, J., Ganju, S., Chen, B., Koul, A., Gupta, A., Skabelkin, A., Zhukov, I., Kumskoy, M., and Nyberg, E.  Learn-to-race: A multimodal control environment for autonomous racing.  In _Proceedings of the IEEE/CVF International Conference on Computer Vision_ , pp. 9793–9802, 2021.
  * Isele et al. (2018) Isele, D., Nakhaei, A., and Fujimura, K.  Safe reinforcement learning on autonomous vehicles.  In _2018 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS)_ , pp. 1–6. IEEE, 2018.
  * Kiran et al. (2021) Kiran, B. R., Sobh, I., Talpaert, V., Mannion, P., Al Sallab, A. A., Yogamani, S., and Pérez, P.  Deep reinforcement learning for autonomous driving: A survey.  _IEEE Transactions on Intelligent Transportation Systems_ , 2021.
  * Li et al. (2021) Li, Q., Peng, Z., Xue, Z., Zhang, Q., and Zhou, B.  Metadrive: Composing diverse driving scenarios for generalizable reinforcement learning.  _arXiv preprint arXiv:2109.12674_ , 2021.
  * Li et al. (2022) Li, Q., Peng, Z., and Zhou, B.  Efficient learning of safe driving policy via human-ai copilot optimization.  _arXiv preprint arXiv:2202.10341_ , 2022.
  * Liu et al. (2021) Liu, Y., Ding, J., and Liu, X.  Resource allocation method for network slicing using constrained reinforcement learning.  In _2021 IFIP Networking Conference (IFIP Networking)_ , pp. 1–3. IEEE, 2021.
  * Luenberger et al. (1984) Luenberger, D. G., Ye, Y., et al.  _Linear and nonlinear programming_ , volume 2.  Springer, 1984.
  * Ma et al. (2021) Ma, H., Guan, Y., Li, S. E., Zhang, X., Zheng, S., and Chen, J.  Feasible actor-critic: Constrained reinforcement learning for ensuring statewise safety.  _arXiv preprint arXiv:2105.10682_ , 2021.
  * Mnih et al. (2015) Mnih, V., Kavukcuoglu, K., Silver, D., Rusu, A. A., Veness, J., Bellemare, M. G., Graves, A., Riedmiller, M., Fidjeland, A. K., Ostrovski, G., et al.  Human-level control through deep reinforcement learning.  _nature_ , 518(7540):529–533, 2015.
  * Paternain et al. (2022) Paternain, S., Calvo-Fullana, M., Chamon, L. F., and Ribeiro, A.  Safe policies for reinforcement learning via primal-dual methods.  _IEEE Transactions on Automatic Control_ , 2022.
  * Ray et al. (2019) Ray, A., Achiam, J., and Amodei, D.  Benchmarking safe exploration in deep reinforcement learning.  _arXiv preprint arXiv:1910.01708_ , 7:1, 2019.
  * Sutton & Barto (1998) Sutton, R. S. and Barto, A. G.  _Reinforcement learning: An introduction_.  MIT press, 1998.
  * Thananjeyan et al. (2021) Thananjeyan, B., Balakrishna, A., Nair, S., Luo, M., Srinivasan, K., Hwang, M., Gonzalez, J. E., Ibarz, J., Finn, C., and Goldberg, K.  Recovery rl: Safe reinforcement learning with learned recovery zones.  _IEEE Robotics and Automation Letters_ , 6(3):4915–4922, 2021.
  * Vinyals et al. (2019) Vinyals, O., Babuschkin, I., Czarnecki, W. M., Mathieu, M., Dudzik, A., Chung, J., Choi, D. H., Powell, R., Ewalds, T., Georgiev, P., et al.  Grandmaster level in starcraft ii using multi-agent reinforcement learning.  _Nature_ , 575(7782):350–354, 2019.
  * Wen et al. (2020) Wen, L., Duan, J., Li, S. E., Xu, S., and Peng, H.  Safe reinforcement learning for autonomous vehicles through parallel constrained policy optimization.  In _2020 IEEE 23rd International Conference on Intelligent Transportation Systems (ITSC)_ , pp. 1–7. IEEE, 2020.
  * Yang et al. (2022) Yang, T.-Y., Zhang, T., Luu, L., Ha, S., Tan, J., and Yu, W.  Safe reinforcement learning for legged locomotion.  _arXiv preprint arXiv:2203.02638_ , 2022.
  * Yuan et al. (2021) Yuan, Z., Hall, A. W., Zhou, S., Brunke, L., Greeff, M., Panerati, J., and Schoellig, A. P.  safe-control-gym: a unified benchmark suite for safe learning-based control and reinforcement learning.  _arXiv preprint arXiv:2109.06325_ , 2021.
  * Zhang et al. (2022) Zhang, L., Shen, L., Yang, L., Chen, S.-Y., Yuan, B., Wang, X., and Tao, D.  Penalized proximal policy optimization for safe reinforcement learning.  _arXiv preprint arXiv:2205.11814_ , 2022.
  * Zhang et al. (2020) Zhang, Y., Vuong, Q., and Ross, K. W.  First order constrained optimization in policy space.  _arXiv preprint arXiv:2002.06506_ , 2020.
  * Zhao et al. (2021) Zhao, W., He, T., and Liu, C.  Model-free safe control for zero-violation reinforcement learning.  In _5th Annual Conference on Robot Learning_ , 2021.