[a]James C. Osborn

# Tuning HMC parameters with gradients

###### Abstract

We investigate the effectiveness of tuning HMC parameters using information
from the gradients of the HMC acceptance probability with respect to the
parameters. In particular, the optimization of the trajectory length and
parameters for higher order integrators will be studied in the context of pure
gauge and dynamical fermion actions.

Gauge field generation in lattice field theory (LFT) generally relies on the
Hybrid (or Hamiltonian) Monte Carlo (HMC) algorithm [1] which integrates the
Hamiltonian equations of motion to generate a new trial configuration. Higher-
order integration schemes have been shown to improve the efficiency of HMC
simulations for LFT applications [8]. These higher-order schemes can have
several parameters that, in general, should be tuned for the specific
application. Furthermore, the action can have preconditioners, such as that by
Hasenbusch [2], which introduce more free parameters to tune. While the tuning
can be, and often is, done by hand by monitoring the acceptance rate and/or
norms of the force terms, this can become tedious as the number of parameters
increases. Here we investigate automating the tuning process using techniques
commonly employed in machine learning applications.

##  1 HMC

The Hybrid Monte Carlo method proceeds by first choosing a random Gaussian
momentum field, $p$, which is associated with the gauge field degrees of
freedom, $U$. Then the equations of motion for the Hamiltonian
$H(p,U)=\frac{1}{2}p^{2}+S(U)$, with $S(U)$ the lattice action to be sampled
from, are approximately integrated for some integration time $\tau$. This
produces a new set of fields

$\displaystyle\left(\begin{array}[]{c}p^{\prime}\\\
U^{\prime}\end{array}\right)=F_{\alpha}\left(\begin{array}[]{c}p\\\
U\end{array}\right)$ (5)

where $F_{\alpha}$ is the operator implementing the integration scheme and
$\alpha$ represents the integration time, $\tau$, and any other tunable
integration or action parameters. The new configuration is conditionally
accepted with the probability

$\displaystyle P_{a}(p^{\prime},U^{\prime},p,U)=\min[1,\exp(-\Delta H)]$ (6)

with

$\displaystyle\Delta H=H(p^{\prime},U^{\prime})-H(p,U)~{}.$ (7)

###  1.1 Improved integrators

For simplicity we’ll consider integrators with a single force type which could
be either a gauge force or a combined gauge and fermion force. The integrator
$F_{\alpha}$ is an approximation of $\exp(\tau[A+B])$ where the operator $A$
is the gauge field update and $B$ is the momentum update using the force term.
This update can be divided into $n$ copies of a base integrator as
$\exp(\tau[A+B])=[\exp(\epsilon[A+B])]^{n}$ with $\epsilon=\tau/n$.

To start, we consider a simple two-step base integrator given by the
approximation

$\displaystyle e^{\epsilon(A+B)}\approx e^{\epsilon\lambda
A}e^{(\epsilon/2)B}e^{\epsilon(1-2\lambda)A}e^{(\epsilon/2)B}e^{\epsilon\lambda
A}$ (8)

where $\lambda$ is a free parameter to tune. Omelyan, et al. [5] gave an
optimal value for this parameter of $\lambda\approx 0.193$, under the
assumption that the two operators that appear in the leading order error term
have equal magnitudes. While this assumption produces a reasonably good
integrator, for LFT applications, the magnitudes may not be very close, and
one can typically reduce the error by tuning the parameter for specific
actions [8].

One can also consider higher order integrators, such as force-gradient
variants. The force-gradient integrator we will consider here is

$\displaystyle e^{\epsilon(A+B)}\approx e^{\epsilon\theta A}e^{\epsilon\lambda
B}e^{(\epsilon/2)(1-2\theta)A}e^{\epsilon(1-2\lambda)B+\epsilon^{3}\chi
C}e^{(\epsilon/2)(1-2\theta)A}e^{\epsilon\lambda B}e^{\epsilon\theta A}$ (9)

where $C=[B,[A,B]]$ is the force-gradient term. For simplicity we implement
this term using the method introduced in [9]. This force-gradient update has
three parameters, $\theta$, $\lambda$, and $\chi$, that can be tuned, however
the cancellation of the order $\epsilon^{3}$ errors places a constraint on two
of the parameters, leaving only one free parameter. Here, for testing
purposes, we will leave all three unconstrained when tuning in order to
evaluate the tuning process, and to allow for a compromise in reducing the
order $\epsilon^{5}$ error at the expense of not completely canceling the
order $\epsilon^{3}$ errors.

In full LFT simulations with dynamical quarks, there can also be different
integration scales used for the gauge force and the fermion forces, using a
recursive integrator [7], which can introduce many more parameters, although
we will not consider those here. As the number of free parameters in the
integration increases, the prospect of automated tuning becomes a particularly
attractive alternative to the conventional methods of tuning by hand.

##  2 Loss function

To measure the efficiency of the integrator, we consider a cost function based
on the effective integration time. For a set of $N$ HMC update steps, each of
integration time $\tau$, the total effective integration time is

$\displaystyle\tau\sqrt{\langle P_{a}\rangle N}$ (10)

where $\langle P_{a}\rangle$ is the average acceptance probability and the
square root comes from the random walk due to choosing a new momentum at the
beginning of each HMC update. The number of update steps needed to go $T$
effective integration time units is then

$\displaystyle N_{T}=\frac{T^{2}}{\langle P_{a}\rangle\tau^{2}}~{}.$ (11)

We can then define the cost of updates for an effective integration time $T$
as the cost per update of length $\tau$, times $N_{T}$. For simplicity we set
the cost per update based on the number of force evaluations, $N_{force}$, and
take $T=1$. This gives the final cost metric

$\displaystyle Cost=N_{T}N_{force}=\frac{N_{force}}{\langle
P_{a}\rangle\tau^{2}}~{}.$ (12)

We want to minimize the cost function (12), however this quantity has an
average in the denominator which makes it inconvenient to work with. Instead
we minimize the loss function

$\displaystyle Loss=-\langle P_{a}\rangle\tau^{2}$ (13)

which has a similar effect and is much easier to work with. For the
minimization we use machine learning (ML) methods and use the gradient of the
loss function to update the HMC parameters. We use the Adam optimizer [4] and
accumulate the gradient during the update using a single update stream (batch
size of one). Since the Adam optimizer already accumulates the gradient
information during the optimization process, we do not need to explicitly
include the ensemble average of $P_{a}$ in the loss function, and instead just
calculate the gradient of $P_{a}$ for the current HMC update. We then update
the tunable parameters based on the accumulated gradient at the end of every
HMC update.

##  3 Calculating gradients

Optimizing the parameters $\alpha$ requires calculating the gradient of the
loss function in Eq. (13) with respect to each of the parameters $\alpha$.
This is typically done through back propagation, which starts from the
definition of the loss function and evaluates the chain rule back through all
the intermediate fields.

The main task is evaluating the gradient of $P_{a}$, which is a function of
$H(p^{\prime},U^{\prime})\equiv H^{\prime}$. Note that if $\Delta H\leq 0$
then the gradient of $P_{a}$ is zero. In the other case, we consider the
integrator as a series of individual update steps. For example, for $n$ copies
of the two-step integrator, written in operator notation, we have

$\displaystyle\left(\begin{array}[]{c}p^{\prime}\\\
U^{\prime}\end{array}\right)=\left[e^{\epsilon\lambda
A}e^{(\epsilon/2)B}e^{\epsilon(1-2\lambda)A}e^{(\epsilon/2)B}e^{\epsilon\lambda
A}\right]^{n}\left(\begin{array}[]{c}p\\\ U\end{array}\right)~{}.$ (18)

We can calculate the gradients from chain rule

$\displaystyle\frac{\partial
H^{\prime}}{\partial\alpha}=\sum_{k}\frac{\partial H^{\prime}}{\partial
q_{k}}\left.\frac{q_{k}}{\partial\alpha}\right|_{q_{k-1}}$ (19)

where

$\displaystyle q_{k}=\left(\begin{array}[]{c}p_{k}\\\ U_{k}\end{array}\right)$
(22)

is the field state after the $k$th integrator step (either applying $A$ or
$B$) starting with $q_{0}=(p,U)$. The first term in the chain rule, Eq. (19),
is accumulated from the end of the update working backwards using

$\displaystyle\frac{\partial H^{\prime}}{\partial q_{k}}=\frac{\partial
H^{\prime}}{\partial q_{k+1}}\frac{\partial q_{k+1}}{\partial q_{k}}$ (23)

and the second term is the gradient of the individual update step with respect
to the parameters keeping the input fields, $q_{k-1}$, fixed. One then needs
to calculate the gradients for the individual update steps which we do not
include here. One issue with the calculation of the gradients is the need to
save all the intermediate fields during an update step. This requires a large
increase in the memory usage, however since the HMC is typically scaled out to
many nodes to reduce the time per update as much as possible, memory use is
not likely to be a limiting factor.

##  4 Results

The code was implemented using the QEX LFT framework [6, 3]. We are using a
simple “tape” implementation which saves a list of operations forming the HMC
update. The tape is then run in the forwards direction to perform the HMC
update, and can be run in the backwards direction to produce the gradient of
the full HMC update step. All runs were done on a relatively small
$12^{3}\times 24$ lattice running on a single desktop.

###  4.1 Pure gauge results

The initial tests were performed on a pure gauge theory with a plaquette
action at $\beta=5.6$. For all runs we started from a thermalized
configuration and then ran 200 tuning updates and 400 measurement updates. The
number of tuning updates was taken as a constant for all runs for simplicity
and was chosen to ensure that all runs were sufficiently tuned. In many cases
the number of tuning updates could have been lowered without impacting
results.

Figure 1:  Pure gauge simulation with single copy ($n=1$, $\epsilon=\tau$) of
a two-step integrator with and without tuning on the trajectory length,
$\tau$, and the integrator parameter $\lambda$. Left: RMS integration error
versus trajectory length. Right: Cost versus trajectory length.

We first examined tuning HMC using an integrator with a single copy of a two-
step update with the field update first. Here there are two parameters we can
tune, the trajectory length, $\tau$, (which is the same as the step size,
$\epsilon$, for a single copy) and the integrator parameter $\lambda$.

On the left panel of Figure 1 we show results for the root mean squared (RMS)
integrator error, $\sqrt{\langle\Delta H^{2}\rangle}$, versus $\tau$ for tuned
and untuned runs. The green x’s are for no tuning, using the Omelyan, et al.
value of $\lambda$. The error scales as $\tau^{3}$ for small $\tau$ as
expected for a second order integrator, and is demonstrated by the purple
line. If we tune $\tau$, then we end up with the points marked with blue
stars, which are at larger values of $\tau$ where the asymptotic $\tau^{3}$
scaling starts to break down. Each star represents an independent tuning run.
The black triangles are for tuning $\lambda$ at each value of $\tau$. Again
this scales as $\tau^{3}$ for small $\tau$, but now with a smaller coefficient
as shown by the red line. The blue circles are the values obtained from runs
tuning both $\lambda$ and $\tau$. Again the tuned value of $\tau$ ends up
right around the value where the asymptotic scaling with $\tau$ starts to
break down.

On the right panel of Figure 1 we show results for the cost function, Eq.
(12), versus $\tau$ for tuned and untuned runs. For the case where $\lambda$
wasn’t tuned, we can see that tuning $\tau$ was able to find values of $\tau$
that were near the minimum of the cost function. For the case of tuned
$\lambda$, the minimum of the cost is very shallow, so the spread of the runs
with tuned $\tau$ is larger, but the values found were still close to optimal.

Figure 2: . Pure gauge simulation with single copy ($n=1$, $\epsilon=\tau$) of
a two-step integrator. Left: Tuned parameters. Right: Parameter paths during
tuning.

In the left panel of Figure 2 we show the final values of the tuned parameters
for several tuning runs with and without tuning $\lambda$. The tuned values
are fairly consistent within each case, but there is a large difference
between the two cases.

In the right panel of Figure 2 we show the trajectories taken in the parameter
space during tuning when both parameters are tuned, all starting from the same
point. While the trajectories may look very different, the final values are
fairly consistent. We also saw a similar convergence when using different
starting parameters.

Figure 3:  Cost versus integrator step size for pure gauge simulation with
multiple copies ($n=$ # ABABA) of a two-step integrator with all parameters
tuned.

In Figure 3 we show the cost versus step size for the fully tuned HMC with an
integrator containing different numbers of copies of the two-step integrator.
We can see that increasing the number of copies decreases the cost up to the
16 copies used here. Also, the tuned value of $\epsilon$ varies with the
number of copies.

Figure 4:  Pure gauge simulation with single copy ($n=1$, $\epsilon=\tau$) of
a force-gradient integrator. Left: RMS integration error versus trajectory
length. Right: Cost versus trajectory length.

In the left panel of Figure 4 we show results for the RMS integration error
versus $\tau$ for the pure gauge action using the force-gradient integrator.
We found that the tuning rate used in the Adam optimizer had a large effect on
the final result when $\tau$ was held fixed at a small value. We had to
decrease the tuning rate from the default of $10^{-3}$ to $10^{-5}$ in order
to see the expected $\tau^{5}$ scaling of the error when tuning the
parameters. However the error with the untuned parameters, did achieve the
expected scaling and performed at least as well as the tuned case. So in this
case, tuning the force-gradient integrator did not seem to be necessary.

In the right panel of Figure 4 we show results for the cost versus $\tau$ for
the pure gauge action using the force-gradient integrator. The runs with
tuning $\tau$ were consistently able to find the minimal cost value as seen
from comparing to the runs with different fixed values of $\tau$.

###  4.2 Staggered quark results

Next we performed tests using a single copy (4 tastes) of unimproved staggered
quarks with the plaquette gauge action. The lattice volume, $12^{3}\times 24$,
and $\beta=5.6$ were the same as in the pure gauge runs. The staggered mass
was set to $m=0.04$ which gives a pion mass of about $0.5$ in lattice units.
Since the fermion force is much more expensive than the gauge force, we only
count the number of fermion force computations in calculating the cost
function. We also put the gauge and fermion forces on the same integrator
timescale for simplicity, instead of using a nested integrator.

Figure 5:  Cost versus trajectory length for staggered quark simulations with
a single integrator copy ($n=1$, $\epsilon=\tau$). Left: Two-step integrator.
Right: Force-gradient integrator.

Figure 5 shows the cost versus $\tau$ for staggered fermion HMC with the two-
step (left panel) and force-gradient (right panel) integrators. Note that the
force-gradient integrator can be generalized such that the integrator
parameters $\lambda$ and $\chi$ have different values for the gauge and
fermion forces, and we allowed them to be tuned independently. We find that
tuning $\tau$ is able to reliably find the minimal cost and that tuning the
other HMC parameters gives a noticeable improvement. Note that for the case of
staggered quarks, as opposed to pure gauge, we see a significant improvement
when tuning the force-gradient integrator.

Figure 6:  Cost versus integrator step size for staggered quark simulations
with multiple copies ($n=$ # ABABA) of a two-step integrator with all
parameters tuned.

In Figure 6 we show the cost versus step size for staggered HMC with multiple
copies of the two-step integrator. As with the pure gauge case, we see that
the cost decreases as the number of copies increases up to the 4 copies
tested. Note that for two copies of the two-step integrator, the cost is
similar to that of the tuned force-gradient integrator (which has the same
number of fermion force evaluations), with the force-gradient cost function
being slightly lower. Without tuning however, the force-gradient integrator is
much worse than two copies of the two-step integrator. As the volume increases
we expect that the force-gradient integrator will perform even better, but
only if properly tuned.

##  5 Summary

We explored implementing and optimizing HMC simulations using gradient
information with methods borrowed from ML applications. Overall we found that
tuning HMC parameters using gradient information works well, at least for the
cases tested here. In most cases tuning the HMC parameters made the HMC more
efficient with the only exception being the case of pure gauge action with the
force-gradient integrator. Using the gradients is a very convenient way to
tune HMC, once the initial investment to develop the implementation has been
done.

We also found that the force-gradient integrators may work much better when
tuned, at least in the case of staggered fermions, and could be competitive
with other integrators even on small volumes. We plan to continue testing this
approach with other actions (improved gauge, improved fermions and with
rooting) and with tuning the mass parameters in Hasenbusch preconditioners.

This research was supported by the Exascale Computing Project (17-SC-20-SC), a
collaborative effort of the U.S. Department of Energy Office of Science and
the National Nuclear Security Administration. This work was performed at
Argonne National Laboratory which is supported under Contract DE-
AC02-06CH11357.

## References

  * [1] Simon Duane, A.D. Kennedy, Brian J. Pendleton, and Duncan Roweth.  Hybrid monte carlo.  Physics Letters B, 195(2):216–222, 1987.
  * [2] Martin Hasenbusch.  Speeding up the hybrid monte carlo algorithm for dynamical fermions.  Physics Letters B, 519(1):177–182, 2001.
  * [3] Xiao-Yong Jin and James C. Osborn.  QEX: a framework for lattice field theories.  PoS, ICHEP2016:187, 2017.
  * [4] Diederik P. Kingma and Jimmy Ba.  Adam: A method for stochastic optimization.  Proc. 3rd International Conference for Learning Representations, 2015.
  * [5] I. P. Omelyan, I. M. Mryglod, and R. Folk.  Optimized verlet-like algorithms for molecular dynamics simulations.  Physical Review E, 65(5), May 2002.
  * [6] James C. Osborn et al.  Quantum EXpressions lattice field theory framework.  https://github.com/jcosborn/qex.
  * [7] J. C. Sexton and D. H. Weingarten.  Hamiltonian evolution for the hybrid monte carlo algorithm.  Nuclear Physics B, 380(3):665–677, August 1992.
  * [8] Tetsuya Takaishi and Philippe de Forcrand.  Testing and tuning symplectic integrators for the hybrid monte carlo algorithm in lattice qcd.  Physical Review E, 73(3), March 2006.
  * [9] Hantao Yin and Robert Mawhinney.  Improving DWF Simulations: Force Gradient Integrator and the Mobius Accelerated DWF Solver.  PoS, Lattice 2011:051, 2012.


††thanks: These authors contributed equally.††thanks: These authors
contributed equally.††thanks: These authors contributed equally.

# Chiral switching of many-body steady states in a dissipative Rydberg gas

Chongwu Xie  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  Konghao Sun  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  Kang-Da Wu  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  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  Hefei National Laboratory,
University of Science and Technology of China, Hefei 230088, 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  Hefei National Laboratory,
University of Science and Technology of China, Hefei 230088, People’s Republic
of China  Wei Yi<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  Hefei National Laboratory,
University of Science and Technology of China, Hefei 230088, People’s Republic
of China  Guo-Yong Xiang<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  Hefei National Laboratory,
University of Science and Technology of China, Hefei 230088, People’s Republic
of China

###### Abstract

Dissipative Rydberg gases are an outstanding platform for the investigation of
many-body quantum open systems. Despite the wealth of existing studies, the
non-equilibrium dynamics of dissipative Rydberg gases are rarely examined or
harnessed from the perspective of non-Hermitian physics, which is but
intrinsic to open systems. Here we report the experimental observation of a
chiral switching between many-body steady states in a dissipative thermal
Rydberg vapor, where the interplay of many-body effects and non-Hermiticity
plays a key role. Specifically, as the parameters are adiabatically varied
around a closed contour, depending on the chirality of the parameter
modulation, the Rydberg vapor can change between two collective steady states
with distinct Rydberg excitations and optical transmissions. Adopting a mean-
field description, we reveal that both the existence of the bistable steady
states and chiral dynamics derive from an exceptional structure in the
parameter space, where multiple steady states of the many-body Liouvillian
superoperator coalesce. We demonstrate that both the exceptional structure and
the resulting state-switching dynamics are tunable through microwave dressing
and temperature variations, confirming their reliance on the many-body
dissipative nature of the Rydberg vapor.

##  I Introduction

Many-body quantum open systems host rich collective behaviors that are absent
in isolated settings. Understanding and manipulating these often complicated
phenomena have become an active field of research in recent years, thanks to
the rapid progress in engineering quantum simulators across a range of
physical systems including ultracold atoms, superconducting qubits, trapped
ions, and Rydberg gases [1]. A quantum open system is intrinsically non-
Hermitian [2], and subject to effective non-Hermitian descriptions, in
contrast to their Hermitian counterparts of isolated (or closed) systems. As
such, many intriguing properties of non-Hermitian Hamiltonians also manifest
in the dynamics of quantum open systems, offering unique perspectives and
fresh insights.

An outstanding example here is the branching singularity in non-Hermitian
models known as the exceptional point (EP), where eigenvalues and eigenvectors
of the non-Hermitian Hamiltonian simultaneously coalesce [3, 4, 5, 6, 7, 10,
8, 9, 11]. Intriguingly, these Hamiltonian EPs give rise to critical behaviors
and topological phenomena, which have stimulated extensive studies on EP-
enhanced sensing [12, 13, 14], laser-mode selection [15, 16, 17], directional
wave-packet transport [18], and chiral state transfer [19, 20, 21, 24, 31, 22,
23, 25, 26, 27, 30, 28, 29]. In quantum open systems where the density-matrix
dynamics is driven by the Liouvillian superoperator, EPs can also emerge in
the Liouvillian eigenspectrum (dubbed Liouvillian EPs) [32, 33, 34], and are
responsible for the recently observed chiral-state transfer in transient
dynamics [26, 27] and efficient quantum heat engines [35]. However, the
experimental studies of the Liouvillian-EP-related dynamics are limited to the
single-particle case, as an experimentally accessible quantum many-body system
that features EP-induced dynamics remains elusive. Further, since the
Liouvillian EPs typically occur in the excited eigenstates (away from the low-
lying steady states with vanishing Liouvillian eigenvalues), their impact is
often limited to transient dynamics, which presents a practical challenge for
experimental observation.

Figure 1:  Exceptional structures in a thermal Rydberg gas. (a) Schematics of
the coupling scheme, with a microwave dressing. Labeling of the energy levels:
$|1\rangle\equiv|5S_{1/2}\rangle$, $|2\rangle\equiv|5P_{3/2}\rangle$,
$|3\rangle\equiv|50D_{5/2}\rangle$, $|4\rangle\equiv|51P_{3/2}\rangle$. The
parameters $\Omega_{c}$, $\Omega_{p}$, and $\Omega_{\text{mw}}$ represent the
Rabi frequencies of the corresponding couplings, respectively. (b)
Experimental setup, with detailed parameters given in the Methods. The
parameters of lasers and microwave in this experiment can be found in
Appendix. Meaning of the abbreviations: PD, Photoelectric detector; DM,
dichroic mirror; MW, microwave. (c) Numerically calculated total Rydberg-state
population $n=\rho_{33}+\rho_{44}$ of the steady states as a function of
$\Omega_{p}$ and $\Delta_{c}$, without the microwave field. (d) Numerically
calculated total Rydberg-state population $n$ as a function of $\Omega_{mw}$
and $\Delta_{c}$ with microwave field. (e)(f) Experimentally measured phase
diagrams under the parameters of (c)(d), respectively. The color bar
represents the normalized transmission difference, $\delta T$, between the
forward and backward scans of $\Delta_{c}$. The green lines represent the
numerically fitted second-order ELs, and the yellow star is the third order
EP. The parameters are: (c)(e) $\Omega_{c}/2\pi=$12.21\text{\,}\mathrm{MHz}$$.
(d)(f) $\Omega_{p}/2\pi=$13.07\text{\,}\mathrm{MHz}$$,
$\Omega_{c}/2\pi=$12.21\text{\,}\mathrm{MHz}$$. See Methods for fitting
parameters used in the theoretical model. (g)(h) Sample data for the
extraction of $\Delta T$, showing hysteresis behavior in the presence of
bistable steady states. The red and blue arrows indicate the scanning
direction, and $\Delta T$ is the transmission difference between two opposite
scanning directions. We take $\Omega_{p}/2\pi=11.55~{}\text{MHz}$ and
$\Omega_{mw}/2\pi=0$ in (g), and $\Omega_{p}/2\pi=13.07~{}\text{MHz}$ and
$\Omega_{mw}/2\pi=15.67~{}\text{MHz}$ in (h). The small peak in (h) at
$\Delta_{c}/2\pi\approx-100$ MHz corresponds to the EIT of the state
$|50D_{3/2}\rangle$.  Figure 2:  Chiral mode switching without the microwave
field. (a1)(c1)/(b1)(d1) The purple dots, connected by the purple lines to
guide the eye, represent the measured transmission along the purple trajectory
in Fig. 1(e). We show samples of the measured hysteresis of the transmission
in the background, to better visualize the steady states. Here, the black dots
represent the initial states, and the arrows indicate the encircling
direction. (a2)-(d2) The corresponding numerically calculated encircling
trajectories on the landscape of the steady-state Rydberg population. The
parameters are the same as those in Fig. 1(c).

In this work, we experimentally demonstrate adiabatic chiral state switching
in a dissipative Rydberg vapor, where, under the many-body effects therein,
Liouvillian EPs emerge in the degenerate steady-state subspace and impact the
long-time dynamics. Characterized by high-lying electronic states and strong
interactions, Rydberg atoms are an ideal candidate for simulating many-body
phenomena such as quantum magnetism [36, 37, 38, 39], topological order [40,
41, 42], as well as many-body dynamics and thermalization [43, 44, 45, 46,
47]. Dissipative Rydberg gases, particularly thermal Rydberg vapors, further
provide direct access to non-equilibrium dynamics in many-body quantum open
systems [48, 49, 50, 51, 52].

Our experiment is based on the bistable transition in a thermal Rydberg vapor,
where two collective steady states feature different Rydberg excitations, and
are characterized by distinct light transmissions in the standard
electromagnetically induced transparency (EIT). While the density-matrix
dynamics of the system can be captured by a many-body Lindblad master equation
under the mean-field approximation, an exceptional structure emerges in the
parameter space where degenerate steady states of the Liouvillian
superoperator coalesce. Ramification of the steady states gives rise to
second-order exceptional lines (ELs) in the parameter space that terminate at
a third-order EP. We experimentally demonstrate that the dynamics near the
exceptional structure exhibits chiral features, exemplified by a chiral state
switching between different steady states—the light transmission can be
switched to a different value when the parameters are adiabatically tuned
around a closed contour to come back to the starting point, but the switching
only occurs one-way around. Such a chiral state switching is the many-body
counterpart of the EP encircling observed in single-particle non-Hermitian
models [20, 21]. It sensitively depends on the many-body parameters of the
Rydberg vapor, and is tunable through microwave dressing and temperature
modulation. Further, since both the chiral dynamics and the exceptional
structure vanish on the single-atom level, our observation is fundamentally
different from the previously studied chiral dynamics in classical non-linear
non-Hermitian systems [54, 55] where the emergence of exceptional structures
in the non-linear case hinges upon the presence of Hamiltonian EPs without
non-linearity. Our findings introduce a fresh non-Hermitian perspective to the
dissipative dynamics in Rydberg atoms, with interesting prospects for future
applications.

Figure 3:  Chiral mode switching with the microwave-field dressing.
(a1)(c1)/(b1)(d1) The cyan (magenta) dots, connected by the cyan (magenta)
lines to guide the eye, represent the measured transmission along the cyan
(magenta) trajectory in Fig. 1(f). We show samples of the measured
transmission hysteresis in the background. The black dots represent the
initial states, and the arrows indicate the encircling direction,
respectively. (a2)-(d2) The corresponding numerically calculated encircling
trajectories on the landscape of the steady-state Rydberg population. The
parameters are the same as those in Fig. 1(d).  Figure 4:  Tuning Chiral mode
switching using temperature. (a) The transmission of the probe laser at the
end of the trajectory, as a function of the temperature. The red (cyan) line
represents measurements from the clockwise (counter-clockwise) encircling
direction. The error bars are the standard deviation for 50 repeated
measurements. (b1)-(b4) The relative location of the exceptional structure to
the fixed trajectory (black loop), as the temperature changes: (b1) 69 oC,
(b2) 70 oC, (b3) 72.5 oC and (b4) 74 oC.

##  II Experimental setup and theoretical model

For our experiment, we use a 87Rb vapor cell at a temperature of $\sim
70^{\circ}$C. As illustrated in Fig. 1(b), the cell is illuminated by a pair
of counter propagating laser beams: a weak probe light with a central
wavelength of 780 nm (resonant with the transition
$5S_{1/2}\Longleftrightarrow 5P_{3/2}$) and a typical power of several micro
Watts; and a strong coupling light with a central wavelength of 480 nm (blue-
detuned to the transition $5P_{3/2}\Longleftrightarrow 50D_{5/2}$). An
additional microwave field, with a central frequency of 17.047 GHz and the
same polarization as that of the incident light, is applied to the atoms to
resonantly couple the Rydberg states $50D_{5/2}$ and $51P_{3/2}$. The coupling
scheme is illustrated in Fig. 1(a), where the relevant internal states are
labeled $|1\rangle\equiv|5S_{1/2}\rangle$, $|2\rangle\equiv|5P_{3/2}\rangle$,
$|3\rangle\equiv|50D_{5/2}\rangle$, $|4\rangle\equiv|51P_{3/2}\rangle$
throughout the work.

In the absence of the microwave field, the atoms undergo the standard EIT in a
ladder-type configuration. Under the many-body effects of the Rydberg gas, the
system is known to host two collective steady states with distinct occupations
in the Rydberg state $|3\rangle$. The addition of the microwave field further
dresses the Rydberg states, affecting the population of the Rydberg
excitations, as well as the properties of the steady states. Invoking a mean-
field treatment of the many-body effects [56], the dynamics of the atomic
ensemble is described by the Lindblad equation

$\displaystyle\dot{\rho}=\mathcal{L}\rho:=-i[H,\rho]+\sum_{i=1}^{6}\left(L_{i}\rho
L_{i}^{\dagger}-\frac{1}{2}\\{L_{i}^{\dagger}L_{i},\rho\\}\right),$ (1)

where $\mathcal{L}$ is the Liouvillian superoperator, and

$\displaystyle H=\begin{pmatrix}0&\Omega_{p}/2&0&0\\\
\Omega^{*}_{p}/2&0&\Omega_{c}/2&0\\\
0&\Omega_{c}^{*}/2&-\tilde{\Delta}_{1}&\Omega_{mw}/2\\\
0&0&\Omega_{mw}^{*}/2&-\tilde{\Delta}_{2}\end{pmatrix},$ (2)

in the basis $\\{|1\rangle,|2\rangle,|3\rangle,|4\rangle\\}$. The dissipative
processes are described by the quantum jump operators $L_{i}$, with
$L_{1}=\gamma_{21}|1\rangle\langle 2|$, $L_{2}=\gamma_{32}|2\rangle\langle
3|$, $L_{3}=\gamma_{31}|1\rangle\langle 3|$, and
$L_{5}=\gamma_{41}|1\rangle\langle 4|$ accounting for the corresponding
spontaneous emissions, respectively. The dephasing of the Rydberg states are
given by $L_{4}=\gamma_{3}|3\rangle\langle 3|$ and
$L_{6}=\gamma_{4}|4\rangle\langle 4|$. The effective detuning
$\tilde{\Delta}_{1}=\Delta_{c}-\alpha\rho_{33}-\beta\rho_{44}$ and
$\tilde{\Delta}_{2}=-\gamma\rho_{33}-\epsilon\rho_{44}$ incorporate the many-
body effects in the Rydberg gas on the mean-field level, which are
proportional to the population of the Rydberg states. The many-body effects
can derive from the long-range Rydberg interactions [48, 50, 56] or ion-
induced Stark shifts [57, 58]. Here
$\gamma_{21},\gamma_{31},\gamma_{32},\gamma_{41}$ and $\gamma_{3},\gamma_{4}$
are the corresponding spontaneous decay rate and dephasing rate, respectively,
and $\\{\alpha,\beta,\gamma,\epsilon\\}$ are the many-body parameters. Note
that in the absence of the microwave field, the Rydberg state $|4\rangle$ is
decoupled, and one should neglect jump operators $L_{5,6}$, focusing only on
the dynamics within the subspace spanned by
$\\{|1\rangle,|2\rangle,|3\rangle\\}$.

Setting $\dot{\rho}=0$ in Eq. (1) allows one to solve for the steady states of
the system, which correspond to the eigenstates of the Liouvillian
superoperator with vanishing eigenvalues. Under the influence of the many-body
effects, the steady state can be degenerate, characterized by multiple
solutions of the Rydberg-state population under the same parameters. In Fig.
1(c)(d), we show the steady-state populations of the Rydberg states
$n=\rho_{33}+\rho_{44}$ (shaded in red), both without and with the microwave
field. In either case, the steady state is three-fold degenerate within a
parameter window (shaded in blue). Further analysis shows that only the
highest- and lowest-lying solutions are stable, indicating bistability.
Nevertheless, at the boundaries of the bistable region (green curves), a
stable steady state coalesces with the unstable one. The boundaries are hence
identified as the second-order Liouvillian ELs, which further terminate at a
third-order EP (yellow pentagrams), where all three steady-state solutions
merge.

Note that the steady-state degeneracy and the Liouvillian exceptional
structures derive from the many-body effects of the Rydberg gas, as they
disappear for vanishing many-body parameters
$\\{\alpha,\beta,\gamma,\epsilon\\}$. The exceptional structure here is
therefore fundamentally different from that in classical non-linear non-
Hermitian models [54, 55, 59, 60], where EPs already emerge in the
eigenspectrum of the linear Hamiltonian and become modified under non-
linearity. Furthermore, despite the emergence of the exceptional structure in
the steady-state space, all steady states possess vanishing Liouvillian
eigenvalues, in contrast to the branching eigenvalues at the Hamiltonian EPs.
The dynamic impact of the exceptional structures therefore manifests at long
times, when the system is in the vicinity of the steady state. This contrasts
previous studies of Liouvillian EPs in the excited eigenstates, which focuses
on the dynamics at short times.

Experimentally, the bistable region outlined by the exceptional structure is
observable through the transmission of the probe laser in the EIT.
Specifically, with other parameters fixed, we scan the detuning $\Delta_{c}$
forth and back, and record the difference in the transmission $\Delta T$. In
Fig. 1(e)(f), we show the experimental data for the normalized transmission
difference $\delta T=\Delta T/\Delta T_{max}$ between the forward and backward
scans, in the planes of $\Omega_{p}-\Delta_{c}$ and $\Omega_{mw}-\Delta_{c}$,
respectively. Here $\Delta T_{max}$ is the maximum value of $\Delta T$
throughout the scanned parameter space. The bistable region is clearly visible
in Fig. 1(e)(f), with non-vanishing $\delta T$. Since $\Delta T$ is positively
correlated to the difference of the Rydberg-state population $n$ between the
two steady states, the boundary of the bistable region corresponds to that of
the Liouvillian exceptional structure, which agrees well with the numerical
fit (green lines) using Eq. (1). In Fig. 1(g)(h), we plot the measured
transmission data under fixed $\Omega_{p}$ and $\Omega_{mw}$. The observed
hysteresis is nothing but the well-known collective bistable transition,
between a steady state with a large Rydberg excitation and transmission, and
another with limited excitation and transmission [48]. Importantly, now that
we identify the bistable transition points with the Liouvillian ELs and EPs, a
natural question arises as to their dynamic consequences.

##  III Chiral mode switching across the EL

A notable dynamic consequence of EPs in non-Hermitian models is the chiral
state transfer. As the system parameters slowly vary along a closed loop in
the vicinity of an EP, the final state of the encircling dynamics depends on
the chirality of the parameter change. The phenomenon is a direct result of
the branching eigenspectral landscape of the non-Hermitian Hamiltonian near
the EP, and has been observed in classical systems with gain and loss [22,
23], and in quantum open settings where the dynamics is dominated by single
qubits [19, 25, 26, 27, 29, 30, 31]. Here we show that chiral mode switching
persists in the many-body setting of Rydberg vapor, thanks to the similarity
of the collective steady-state landscape near the Liouvillian exceptional
structure to that of the eigenspectrum near a Hamiltonian EP.

For our experiment, we initialize the system in one of the steady states
within the bistable region, and slowly vary $\Omega_{p}$ and $\Delta_{c}$
along the purple loop in Fig. 1(e). Each encirclement takes the time $T=40$ s,
and for each data point, we stay the system for $100$ ms before taking data,
to ensure its relaxation to the local steady state. Both timescales are much
larger than the relaxation time $T_{r}\approx 0.2$ ms near the bistable
boundary (see Supplementary Material for detailed protocol). Figure 2(a1)(c1)
show the measured transmission when the initial state is the high-lying steady
state, and in the absence of microwave dressing. Under a clockwise (CW)
parameter change, the measured transmission changes to a lower value on
returning to the initial point, indicating a state switch. By contrast, under
a counterclockwise (CCW) rotation, the measured transmission recovers its
initial value, as the system returns to its initial steady state. The
chirality of the state switching is reversed when the system is initialized in
the low-lying steady state, as illustrated in Fig. 2(b1)(d1). For each case,
we show the corresponding Liouvillian steady-state landscape to the right [see
Fig. 2(a2)(b2)(c2)(d2)], where the purple trajectory corresponds to the
numerically simulated population using Eq. (1) with the encircling period
$T=5$ ms. The experimental observations are in good agreement with numerical
simulations. Note that the trajectory of the encirclement intersects the
second-order EL on either side of the third order EP at least once, which, as
we show later, is important for the occurrence of chiral state switching.

##  IV Tuning the many-body effects

Since the observed chiral state switching is many-body in nature, it is also
subject to controls that impact the many-body effects in the Rydberg vapor.
Indeed, according to Eq. (1) and Fig. 1(f), through a resonant coupling of two
Rydberg states, the microwave dressing significantly changes the bistable
region and the exceptional structure. Further, since the power of the
microwave field can be tuned over a wider range than that of the probe field,
the chiral dynamics can now be examined over a broader range of parameters.

With the microwave dressing, we experimentally demonstrate the chiral state
switching in the parameter space of $\Omega_{mw}$ and $\Delta_{c}$. We vary
the parameters along the closed trajectories shown in Fig. 1(d), all with an
encircling period of $T=20$ s. Facilitated by the improved tunability of
$\Omega_{mw}$, we initialize the system in a high-lying steady state within
the bistable region, but choose different encircling trajectories. As
illustrated in Fig. 3(a1)(c1), the chiral state switching persists when the
trajectory cuts through the bistable region [magenta in Fig. 1(d)], thus
intersecting the EL on either side of the third-order EP once. This is similar
to the trajectory in Fig. 2. However, when the trajectory only intersects with
the EL on one side of the EP [green in Fig. 1(d)], the measured transmission
always recovers the initial value, as the system returns to the same steady
state. This is illustrated in Fig. 3(b1)(d1). Based on these observations, we
summarize a set of sufficient conditions for the occurrence of chiral state
switching: i) the initial state is prepared in the bistable region as one of
the steady states; ii) with respect to the initial position, the two nearest
intersection points between the trajectory and the ELs should be on either
side of the third order EP.

Alternatively, the atomic density of the vapor is tunable through temperature,
which offers a novel control over the chiral state switching. To demonstrate
such a possibility, we fix the encircling contour in the parameter space of
$\Omega_{mw}$ and $\Delta_{c}$, as well as the initial state which is chosen
as the high-lying steady state. We then perform CW and CCW encirclings over a
fixed trajectory at different temperatures. Since the exceptional structure
shifts toward regions with larger detuning $|\Delta_{c}|$ [see Fig. 4(b1-b4)],
varying the temperature adjusts the relative position of the fixed trajectory
to the exceptional structure, enabling the tuning of the chiral state
switching. We now speed up the encircling process and vary the parameters
around the loop within $T=33$ ms. Since this time scale is still much larger
than the local relaxation time, we expect the system to remain close to the
local steady states along the trajectory. We then measure the transmission of
the final state of the CW and CCW rotations. As illustrated in Fig. 4(a), the
chiral state switching emerges in the intermediate temperature range, where
the two conditions for the chiral state switching hold.

##  V Summary

We experimentally demonstrate the chiral state switching between many-body
steady states in a dissipative Rydberg vapor. Through theoretical analysis, we
reveal that the phenomenon, along with the well-known bistable transition in
Rydberg gases, is a direct consequence of the Liouvillian exceptional
structure, consisting of ELs and higher-order EPs. We establish the condition
for the occurrence of the chiral state switching through experiments with
various parameters and encircling trajectories. Importantly, both the
exceptional structure and the chiral dynamics studied here hinge upon the
many-body nature of the Rydberg gas, and are therefore fundamentally different
from previous studies on two- or three-level systems [26, 27]. Such a many-
body origin not only provides novel control schemes in the form of the
microwave dressing and temperature modulation, but also offers a practical
route toward the rich non-Hermitian physics in the quantum many-body regime.
For instance, various extensively researched properties of EPs in non-
Hermitian settings, the EP-enhanced sensing in particular, may find
applications in the context of many-body quantum open systems that are
experimentally accessible.

## References

  * [1] I. M. Georgescu, S. Ashhab, and Franco Nori, Rev. Mod. Phys. 86, 153 (2014).
  * [2] Y. Ashida, Z. Gong, and M. Ueda, Adv. Phys. 69, 249 (2020).
  * [3] T. Kato, $Perturbation\,Theory\,of\,Linear\,Operators$ (Springer, 1966).
  * [4] M. V. Berry, and D. H. J. O’Dell, J. Phys. Math. Gen. 31, 2093 (1998).
  * [5] W. D. Heiss, Eur. Phys. J. D 7, 1 (1999).
  * [6] M. Berry, Czech. J. Phys. 54, 1039 (2004).
  * [7] W. D. Heiss, J. Phys. A 37, 2455 (2004).
  * [8] M.-A. Miri, and A. Alù, Science 363, 42 (2019).
  * [9] Ş. K. Özdemir, S. Rotter, F. Nori and L. Yang, Nat. Mat. 18, 783 (2019).
  * [10] K. Ding, C. Fang, and G.-C. Ma, Nat. Rev. Phys. 4, 745 (2022).
  * [11] A.-D. Li, H. Wei, M. Cotrufo, W.J. Chen, S. Mann, X. Ni, B.-C. Xu, J.-F. Chen, J. Wang, S.-H. Fan, C.-W. Qiu, A. Alù, and L. Chen, Nat. Nanotech. 18, 706 (2023).
  * [12] W.-J. Chen, Ş. K. Özdemir, G.-M. Zhao, J. Wiersig, and L. Yang, Nature 548, 192 (2017).
  * [13] Y.-H. Lai, Y.-K. Lu, M.-G. Suh, Z.-Q. Yuan, and K. Vahala, Nature 576, 65 (2019).
  * [14] M.-Z. Zhang, W. Sweeney, C. W. Hsu, L. Yang, A. D. Stone, and L. Jiang, Phys. Rev. Lett. 123, 180501 (2019).
  * [15] B. Peng, Ş. K. Özdemir, S. Rotter, H. Yilmaz, M. Liettzer, F. Monifi, C. M. Bender, F. Nori, and L. Yang, Science 346 6207 (2014).
  * [16] L. Feng, Z.-J. Wong, R.-M. Ma, Y. Wang, and X. Zhang, Science 346 6207 (2014).
  * [17] H. Hodaei, M. Miri, M. Heinrich, D. N. Christodoulides, and M. Khajavikhan, Science 346 6212 (2014).
  * [18] B. Peng, Ş. K. Özdemir, M. Liertzer, W.-J. Chen, J. Kramer, H. Yilmaz, J. Wiersig, S. Rotter, and L. Yang, PANS 113, 25 (2016).
  * [19] R. Lefebvre, O. Atabek, M. Šindelka, and N. Moiseyev, Phys. Rev. Lett. 103, 123003 (2009).
  * [20] R. Uzdin, A. Mailybaev, and N. Moiseyev, J. Phys. A: Math. Theor. 44, 435302 (2011).
  * [21] T. J. Milburn, J. Doppler, C. A. Holmes, S. Portolan, S. Rotter, and P. Rabl, Phys. Rev. A 92, 052124 (2015).
  * [22] J. Doppler, A. A. Mailybaev, J. Böhm, U. Kuhl, A. Girschik, F. Libisch, T. J. Milburn, P. Rabl, N. Moiseyev, and S. Rotter, Nature 537, 76 (2016).
  * [23] X.-L. Zhang, S.-B. Wang, B. Hou, and C. T. Chan, Phys. Rev. X 8, 021066 (2018).
  * [24] Ş. K. Özdemir, S. Rotter, F. Nori, and L. Yang, Nat. Mat. 18, 783 (2019).
  * [25] Q.-H. Song, M. Odeh, J. Zúñiga-Pérez, B. Kanté, and P. Genevet, Science 373, 6559 (2021).
  * [26] W.-J. Chen, M. Abbasi, Yogesh N. Joglekar, and K. W. Murch, Phys. Rev. Lett. 127, 140504 (2021).
  * [27] W.-J. Chen, M. Abbasi, B. Ha, S. Erdamar, Y. N. Joglekar, and K. W. Murch, Phys. Rev. Lett. 128, 110402 (2022).
  * [28] W.-Q. Liu, Y. Wu, C.-K. Duan, X. Rong, and J.-F. Du, Phys. Rev. Lett. 126, 170506 (2021).
  * [29] L.-Y. Ding, K.-Y. Shi, Y.-X. Wang, Q.-X. Zhang, C.-H. Zhu, L.-D. Zhang, J.-Q. Yi, S.-N. Zhang, X. Zhang, and W. Zhang, Phys. Rev. A. 105, L010204 (2022).
  * [30] Z.-J. Ren, D. Liu, E.-T. Zhao, C.-D. He, K. K. Pak, J. Li, and G.-B. Jo, Nat. Phys. 18, 385 (2022).
  * [31] H. Nasari, G. Lopez-Galmiche, H. E. Lopez-Aviles, A. Schumer, A. U. Hassan, Q. Zhong, S. Rotter, P. LiKamWa, D. N. Christodoulides, and M. Khajavikhan, Nature 605, 256 (2022).
  * [32] F. Minganti, A. Miranowicz, R. W. Chhajlany, and F. Nori, Phs. Rev. A 100, 062131 (2019).
  * [33] F. Minganti, A. Miranowicz, R. W. Chhajlany, I. I. Arkhipov, and F. Nori, Phys. Rev. A 101, 062112 (2020).
  * [34] P. Kumar, K. Snizhko, and Y. Gefen, Phys. Rev. A 104, L050405 (2021).
  * [35] J.-W. Zhang, J.-Q. Zhang, G.-Y. Ding, J.-C. Li, J.-T. Bu, B. Wang, L.-L. Yan, S.-L. Su, L. Chen, F. Nori, Ş. K. Özdemir, F. Zhou, H. Jing, and M. Feng, Nat.Commun 13 6225 (2022).
  * [36] H. Labuhn, D. Barredo, S. Ravets, S. d. Léséleuc, T. Macrì, T. Lahaye, and A. Browaeys, Nature 534, 667 (2016).
  * [37] H. Bernien, S. Schwartz, A. Keesling, H. Levine, A. Omran, H. Pichler, S. Choi, A. S. Zibrov, M. Endres, M. Greiner, V. Vuletić, and M. D. Lukin, Nature 551 579 (2017).
  * [38] V. Lienhard, S. d. Léséleuc, D. Barredo, T. Lahaye, A. Browaeys, M. Schuler, L.-P. Henry, and A. M. Läuchli, Phys. Rev. X 8, 021070 (2018).
  * [39] P. Scholl, M. Schuler, H. J. Williams, A. A. Eberharter, D. Barredo, K.-N. Schymik, V. Lienhard, L.-P. Henry, T. C. Lang, T. Lahaye, A. M. Läuchli, and A. Browaeys, Nature 595, 233 (2021).
  * [40] R. Verresen, M. D. Lukin, and A. Vishwanath, Phys. Rev. X 11, 031005 (2021).
  * [41] G. Semeghini, H. Levine, A. Keesling, S. Ebadi, T. T. Wang, D. Bluvstein, R. Verresen, H. Pichler, M. Kalinowski, R. Samajdar, A. Omran, S. Sachdev, A. Vishwanath, M. Greiner, V. Vuletić, and M. D. Lukin, Science 374, 6572 (2021).
  * [42] R. Samajdar, D. G. Joshi, Y.-T. Teng, and S. Sachdev, Phys. Rev. Lett. 130, 043601 (2023).
  * [43] E. Guardado-Sanchez, P. T. Brown, D. Mitra, T. Devakul, D. A. Huse, P. Schauß, and W. S. Bakr, Phys. Rev. X 8, 021069 (2018).
  * [44] Y. Park, K. Kim, H.-S. Sim, and J. Ahn, Phys. Rev. Lett. 120, 180502 (2018).
  * [45] A. A. Michailidis, C. J. Turner, Z. Papić, D. A. Abanin, and M. Serbyn, Phys. Rev. X 10, 011055 (2020).
  * [46] D. Bluvstein, A. Omran, H. Levine, A. Keesling, G. Semeghini, S. Ebadi, T. T. Wang, A. A. Michailidis, N. Maskara, W. W. Ho, S. Choi, M. Serbyn, M. Greiner, V. Vuletić, and M. D. Lukin, Science 371, 6536 (2021).
  * [47] C. J. Turner, J.-Y. Desaules, K. Bull, and Z. Papić, Phys. Rev. X 11, 021021 (2021).
  * [48] C. Carr, R. Ritter, C. G. Wade, C. S. Adams, and K. J. Weatherill, Phys. Rev. Lett. 111, 113901 (2013).
  * [49] N. Šibalić, C. G. Wade, C. S. Adams, K. J. Weatherill, and T. Pohl, Phys. Rev. A 94, 011401(R) (2016).
  * [50] N. R. de Melo, C. G. Wade, N. Šibalić, J. M. Kondo, C. S. Adams, and K. J. Weatherill, Phys. Rev. A 93, 063863 (2016).
  * [51] D.-S. Ding, H. Busche, B.-S. Shi, G.-C. Guo, and C. S. Adams, Phys. Rev. X 10, 021023 (2020).
  * [52] D.-S. Ding, Z.-K. Liu, B.-S. Shi, G.-C. Guo, K. Mølmer, and C. S. Adams, Nat. Phys. 18, 1447 (2022).
  * [53] K. Wadenpfuhl and C. S. Adams, Phys. Rev. Lett. 131, 143002 (2023).
  * [54] H.-W. Wang, S. Assawaworrarit, and S.-H. Fan, Opt. Lett. 44, 3 (2019).
  * [55] Z. Li, C. Li, G. Xu, W. Chen, Z. Xiong, H. Jing, J.-S. Ho, and C.-W. Qiu, Sci. Adv 9, 27 (2023).
  * [56] C. G. Wade, M. Marcuzzi, E. Levi, J. M. Kondo, I. Lesanovsky, C. S. Adams, and K. J. Weatherill, Nat. Commun. 9, 3567 (2018).
  * [57] D. Weller, A. Urvoy, A. Rico, R. Löw, and H. Kübler, Phys. Rev. A 94, 063820 (2016).
  * [58] D. Weller, J. P. Shaffer, T. Pfau, R. Löw, and H. Kübler, Phys. Rev. A 99, 043418 (2019).
  * [59] S. Assawaworrarit, X.-F. Yu, and S.-H. Fan, Nature 546, 387 (2017).
  * [60] K. Bai, L. Fang, T.-R. Liu, J.-Z. Li, D.-D. Wan, and M. Xiao, Natl. Sci. Rev 10, 7 (2022).

Acknowledgement. This work was supported by the National Natural Science
Foundation of China (Grants Nos. 12134014, 61905234, 11974335, 12374337, and
12374479), the Key Research Program of Frontier Sciences, CAS (Grant No.
QYZDYSSW-SLH003), USTC Research Funds of the Double First-Class Initiative
(Grant No. YD2030002007) and the Fundamental Research Funds for the Central
Universities (Grant No. WK2470000035).

Supplemental Material for "Chiral switching of many-body steady states in a
dissipative Rydberg gas"

##  Appendix A Theoretical fitting

For the boundary of the bistable region in the absence of the microwave field,
we fix $\Omega_{c}/2\pi=$12.21\text{\,}\mathrm{MHz}$$ according to the
experiment, and vary $\gamma_{21}$, $\gamma_{32}$, $\gamma_{31}$,
$\gamma_{3}$, and $\alpha$ to fit the calculated ELs with the measured
boundaries of the bistable region. The fitted values of the parameters are:
$\gamma_{21}/2\pi=$27.22\text{\,}\mathrm{MHz}$$,
$\gamma_{32}/2\pi=$0.29\text{\,}\mathrm{MHz}$$,
$\gamma_{31}/2\pi=$0.016\text{\,}\mathrm{MHz}$$,
$\gamma_{3}/2\pi=$20.31\text{\,}\mathrm{MHz}$$ and
$\alpha/2\pi=$-273.89\text{\,}\mathrm{MHz}$$.

In the presence of the microwave field, we fix
$\Omega_{c}/2\pi=$12.21\text{\,}\mathrm{MHz}$$,
$\Omega_{p}/2\pi=$13.07\text{\,}\mathrm{MHz}$$ according to the experiment,
and vary $\gamma_{21}$, $\gamma_{32}$, $\gamma_{31}$, $\gamma_{3}$, $\alpha$,
$\beta$, $\gamma$, and $\epsilon$ to fit the calculated ELs with the measured
boundaries. We further require $\gamma_{41}=\gamma_{31}$ and
$\gamma_{3}=\gamma_{4}$ for the fitting process. The fitted parameters are:
$\gamma_{21}/2\pi=$7.34\text{\,}\mathrm{MHz}$$,
$\gamma_{32}/2\pi=$0.70\text{\,}\mathrm{MHz}$$,
$\gamma_{31}/2\pi=\gamma_{41}/2\pi=$0.13\text{\,}\mathrm{MHz}$$,
$\gamma_{3}/2\pi=\gamma_{4}/2\pi=$11.78\text{\,}\mathrm{MHz}$$,
$\alpha/2\pi=$-174.32\text{\,}\mathrm{MHz}$$,
$\beta/2\pi=$-954.93\text{\,}\mathrm{MHz}$$,
$\epsilon/2\pi=$1114.08\text{\,}\mathrm{MHz}$$, and
$\gamma/2\pi=$6843.66\text{\,}\mathrm{MHz}$$.

Stability of the steady state To study the stability of the steady-state
solution, we consider a small deviation from the steady-state solution
$\rho^{\prime}$, with $\rho=\rho^{\prime}+\Delta\rho$. Here $\rho^{\prime}$
satisfies $\mathcal{L}[\rho^{\prime}]=0$. We then expand Eq. (1) to the first
order in $\Delta\rho$ around $\rho^{\prime}$. For instance, in the absence of
the microwave field, we have

$\displaystyle\Delta\dot{\rho}=J\left[\rho^{\prime}\right]\Delta\rho,$ (S1)

where

$\displaystyle\leavevmode\resizebox{433.62pt}{}{$J\left[\rho^{\prime}\right]=\begin{pmatrix}0&0&-\Omega_{p}&0&0&\gamma_{21}&0&0&\gamma_{31}\\\
0&-\frac{\gamma_{21}}{2}&0&0&-\frac{\Omega_{c}}{2}&0&0&0&0\\\
\frac{\Omega_{p}}{2}&0&-\frac{\gamma_{21}}{2}&\frac{\Omega_{c}}{2}&0&-\frac{\Omega_{p}}{2}&0&0&0\\\
0&0&-\frac{\Omega_{c}}{2}&\frac{1}{2}(-\gamma_{3}-\gamma_{32}-\gamma_{31})&\Delta_{c}-\alpha\rho^{\prime}_{33}&0&0&\frac{\Omega_{p}}{2}&-\alpha\text{Im}\rho^{\prime}_{13}\\\
0&\frac{\Omega_{c}}{2}&0&\alpha\rho^{\prime}_{33}-\Delta_{c}&\frac{1}{2}(-\gamma_{3}-\gamma_{32}-\gamma_{31})&0&-\frac{\Omega_{p}}{2}&0&\alpha\text{Re}\rho^{\prime}_{13}\\\
0&0&\Omega_{p}&0&0&-\gamma_{21}&0&-\Omega_{c}&\gamma_{32}\\\
0&0&0&0&\frac{\Omega_{p}}{2}&0&\frac{1}{2}(-\gamma_{21}-\gamma_{3}-\gamma_{32}-\gamma_{31})&\Delta_{c}-\alpha\rho^{\prime}_{33}&-\alpha\text{Im}\rho^{\prime}_{23}\\\
0&0&0&-\frac{\Omega_{p}}{2}&0&\frac{\Omega_{c}}{2}&\alpha\rho^{\prime}_{33}-\Delta_{c}&\frac{1}{2}(-\gamma_{21}-\gamma_{3}-\gamma_{32}-\gamma_{31})&\frac{1}{2}\left(2\alpha\text{Re}\rho^{\prime}_{23}-\Omega_{c}\right)\\\
0&0&0&0&0&0&0&\Omega_{c}&-\gamma_{32}-\gamma_{31}\\\ \end{pmatrix}$},$ (S2)
$\displaystyle\Delta\rho=\left[\Delta\rho_{11},\mathrm{Re}\Delta\rho_{12},\mathrm{Im}\Delta\rho_{12},\mathrm{Re}\Delta\rho_{13},\mathrm{Im}\Delta\rho_{13},\Delta\rho_{22},\mathrm{Re}\Delta\rho_{23},\mathrm{Im}\Delta\rho_{23},\Delta\rho_{33}\right]^{\text{T}}.$
(S3)

We now examine the eigenvalues of the non-Hermitian matrix: when the real
components of the eigenvalues are all negative, the steady state
$\rho^{\prime}$ is stable; otherwise, $\rho^{\prime}$ is unstable. The stable
(unstable) regimes in Figs. 1(b) and are thus determined. In the presence of
the microwave field, we follow a similar procedure to calculate the stable
(unstable) regimes in Fig. 1(d).

Figures |  $\Omega_{p}/2\pi$ (MHz) |  $\Omega_{mw}/2\pi$ (MHz) | Temperature$~{}(^{o}C)$
---|---|---|---
Figure 2 | 5.51-12.96 | 0 | 70
Figure 3 | 13.07 | 0-94.02 | 70
Figure 4 | 13.07 | 0-94.02 | 69.0-74.0
Table 1:  Experimental parameters. Both the microwave and probe laser are
resonant to the corresponding energy levels in the whole experiment. The power
of the coupling laser is fixed at about 127.3 mW which corresponds to a Rabbi
frequency $\Omega_{c}/2\pi=12.21$ MHz. Both the two lasers are focused to the
same location and have a same $1/e^{2}$ waist of around 100 $\mu$m. All the
lasers and microwave have the same polarization.

##  Appendix B Experimental parameters

Table 1 shows the general experimental parameters. In the demonstration of the
chiral dynamics in the Rydberg system, as shown in Figs. 2 and 3, we change
the parameters every time interval $\tau=0.2$ s. At each time interval, we
wait $0.1$ s for the system to evolve to a steady state and record the
transmission of the probe laser with an integral time of $0.1$ s. The discrete
temporal parameters along the purple loop in Fig. 1(e) are as follows

$\displaystyle\Omega_{p}(n)$
$\displaystyle=\Omega_{0}\sqrt{1+p\cos{\theta_{n}}},\
\Delta_{c}(n)=\Delta_{0}+\delta\sin{\theta_{n}},$ (S4)

where $\Omega_{0}/2\pi=10.58$ MHz, $p=0.31$, $\Delta_{0}/2\pi=-52.80$ MHz,
$\delta/2\pi=14.92$ MHz. $\Omega_{p}(n)$ and $\Delta_{c}(n)$ are the
parameters at the time $t_{n}=n\tau$. The value of $\theta_{n}$ depends on the
encirclement direction: $\theta_{n}=2\pi n/N$ for the clockwise and
$\theta_{n}=2\pi(1-n/N)$ for the counterclockwise directions, respectively.
$n$ is taken from 0 to $N$, and $n=0$ is the starting while $n=N$ is the
ending. Here, we set $N=200$, which decides the time cost in an encirclement.
On the other hand, the discrete temporal parameters along the three loops in
Fig. 1(f) are given by

$\displaystyle\Omega_{mw}(n)$
$\displaystyle=\Omega_{0}-\tilde{\Omega}\cos{\theta_{n}},\
\Delta_{c}(n)=\Delta_{0}-\delta\sin{\theta_{n}},$ (S5)

with different $\Omega_{0},\ \tilde{\Omega},\ \Delta_{0}$ and $\delta$, which
can be found in Tbl. 2. As before, we set $\theta_{n}=2\pi n/N$ for the
clockwise and $\theta_{n}=2\pi(1-n/N)$ for the counterclockwise directions,
respectively, with $N=100$.

loops |  $\Omega_{0}/2\pi$ (MHz) |  $\tilde{\Omega}/2\pi$ (MHz) |  $\Delta_{0}/2\pi$ (MHz) |  $\delta/2\pi$ (MHz)
---|---|---|---|---
black | 34.48 | 18.81 | -57.80 | 7.54
cyan | 34.48 | 18.81 | -57.80 | 7.54 $\Theta(\pi-\theta_{n})$
magenta | 53.27 | 37.60 | -57.80 | 29.86 $\Theta(\theta_{n}-\pi)$
Table 2:  Description of loops. Here $\Theta(x)=1$ when $x>0$ and
$\Theta(x)=0$ when $x\leq 0$. These parameters are used in Eq. (S5)

##  Appendix C Experimental details

Figure S1:  Detailed experimental setup. The whole experimental setup can
enable us to experimentally control several different parameters, e.g., laser
detuning, laser power, MW power. Meaning of the abbreviations: R, RedPitaya;
PD, Photoelectric detector; DM, dichroic mirror; MW, microwave; IF,
interference filter; RF, radio frequency source; EOM, electro-optical
modulator; AOM, acousto-optical modulator; PBS, polarizing beam splitter; BS;
beam splitter; SHG, second harmonic generation; ULE, ultra-low expansion
materials; SAS, saturation absorption spectrum.

Parameter control. The detailed information of our experimental setup for
investigating parameter encircling is detailed in Fig. S1. In our experiment,
$\Omega_{p},\ \Delta_{c},$ and $\Omega_{mw}$ are the three control parameters.
First, a photoelectric detector (PD, PD1 in Fig. S1) is used to monitor the
power of the laser. An acousto-optical modulator (AOM) and a RedPitaya-
designed servo system (R1 in Fig. S1) help lock the power of the probe laser
$\Omega_{p}$ to a fixed value. Redpitaya can function as a lockbox by
utilizing the pyrpl Python module. We can easily modify the set point of the
servo system to vary the value of $\Omega_{p}$ continually using Python codes.
Second, the frequency of the coupling laser is locked on a cavity. Actually,
the coupling laser is generated through a second harmonic generation (SHG)
process using a basic laser with a central wavelength of 960 nm. We split a
little basic laser (about 10 mW) to a frequency lock system before the SHG
process (see Fig. S1). The split basic laser is modulated through an electro-
optic modulator (EOM) before it enters the cavity. The EOM not only makes a
modulation to generate an error signal but can also make a large frequency
shift on the basic laser. That means we can easily modify the frequency of the
basic laser used in the SHG by changing the frequency of the radio frequency
(RF) signal on the EOM. Then the coupling laser detuning $\Delta_{c}$ can be
controlled. However, if we directly switch the frequency of the RF by setting
another value, the RF will output nothing in its switching time. And the
frequency of the coupling laser cannot change continually using such a method.
So we apply the frequency modulation (FM) function to the RF source to control
the frequency of the RF, and $\Delta_{c}$ continually. The FM function can
make a frequency shift proportional to the modulation voltage, and a
continually varying voltage can be easily obtained experimentally using an
arbitrary signal generator (ASG). RedPitaya can work as an ASG, which can be
easily controlled with Python. So we can use a computer to control the
RedPitaya (R2 in Fig. S1), and subsequently $\Delta_{c}$ with the help of the
FM. For the same reason as that of the microwave source, we can apply an
amplitude modulation (AM) to the microwave source to control $\Omega_{mw}$
continually.

Generation of the probe and coupling lasers. A continues laser (Toptica DLC
pro) with a central wavelength of 780 nm and a total power of 80 mW is split
into three branches: the first is used for monitoring the absolute frequency
by saturation absorption spectrum (SAS) technology, the second is sent to
frequency locking to $5S_{1/2}\leftrightarrow 5P_{3/2}$ transition and
linewidth reduction using a high-finesse reference cavity, and the third
branch is used for experiments, i.e., the probe light in the EIT system.
Moreover, the power of the probe light is monitored by a photoelectric
detector (PD1) and stabilized using an AOM.

The coupling light, which has a central wavelength of 480 nm (resonant with
$5P_{3/2}\leftrightarrow 50D_{5/2}$), is generated through a SHG step with
Toptica tapered amplifier (TA) pro. In particular, the TA unit can deliver up
to 1.4 W at 960 nm and is then fed to the frequency conversion module using a
periodically poled potassium titanyl phosphate (PPKTP) crystal, which is anti-
reflecting coated at both 480 nm and 960 nm, and the temperature is stabilized
to around 28∘C to achieve the quasi-phase match condition. As a result, for an
input power of 1.385 W, an output power of 200 mW of the coupling laser is
achieved. Note that we have sacrificed conversion efficiency to achieve a more
stable performance.

Linewidth reduction and frequency control of the probe and coupling lasers.
For reducing technical noise and conveniently controlling the frequency of
both probe and coupling lasers, we use a double-coated 100-mm-long high
finesse-reference cavity made of ultra-low expansion (ULE) material, which is
loaded on a Zerodur mounting block with four contacting viton balls for
vibration isolation. The cavity is placed in a vacuum house (VH, model
VH-6020-4), which is enclosed by loose thermal insulation, with the
temperature stabilized to around 32.06 ∘C. The typical vacuum degree is on the
order of 10-8 mbar after baking. The finesses of 780 nm and 960 nm light are
around 15000 and 20000, respectively. After locking the frequency to the
cavity mode with a fast PID (FALC 110) using PDH technology, the linewidth is
supposed to be around 1 kHz. The whole frequency locking system, except the
locking electronics, is provided by a stable laser system (SLS).

Since the cavity mode cannot be coincident with both the probe and coupling
laser mode, we let both 960 and 780 nm light be coupled into a fiber electro-
optical modulator for generating sidebands with a large frequency range (9
kHz-10 GHz) and locked one of the sidebands of both lasers to one cavity mode.
Thus, the frequency of both the probe and coupling lasers can be arbitrarily
tuned by controlling the RF source.

Vapor cell and temperature control. The vapor cell used in this experiment has
a tiny internal geometric size of $8\times 8\times 8~{}\text{mm}^{3}$. It is
filled with natural rubidium vapor ($27.8\%$ for 87Rb and $72.2\%$ for 85Rb),
but we only used the part of 87Rb. The vapor cell is contained in a specially
designed ceramic oven (see Fig. S2). Ceramics have good thermal conductivity
and do not reflect or absorb microwaves like metals do. The oven is driven by
a thermoelectric cooler (TEC), which is controlled by a TEC controller
(TED4015, Thorlabs).

Figure S2:  Vapor cell and ceramic oven. The vapor cell is placed in a ceramic
oven.

##  Appendix D Additional data

Additional encircling paths. We have only displayed a portion of the
encirclement’s results due to the main text’s space constraints. In fact, we
have achieved the encirclement along a variety of loops experimentally; they
will be shown below as supplemental material. Figure S3 shows the more
detailed encircling results along the loops in Fig. 1(f). Figure S3(a1)-(a4)
do not show a chiral phenomenon where the final states are independent on the
encircling direction or the initial states. As expected, Figs. S3(b1)-(c4)
show chiral switching where the final states are only dependent on the
encircling direction. Since they satisfy the two conditions for the occurrence
of chiral switching.

Figure S3:  Additional encircling. The encircling trajectories are along the
path in Fig. 1(f) with the same color. The meanings of the points, lines and
arrows are the same as those in Figs. 2 and 3.

Relaxation time near the bistable region. When we scan $\Delta_{c}$ across the
bistable zone, the relationship between the transmission and the scanning time
is depicted in Fig. S4. It took the system approximately $0.2$ ms to evolve
from the low-lying to the high-lying state; thus, we can estimate the
approximate relaxation time $T_{r}\approx 0.2$ ms.

Figure S4:  State evolution. The measured transmission of the probe laser when
we scan $\Delta_{c}$. Here the parameters are the same as those in Fig. 1(g).


randomly generating two transition laws $(\bm{P}_{1},\bm{P}_{2})$. However,
the formal proof that the set where (97) is satisfied is dense in
$\mathfrak{P}\times\mathfrak{P}$ is novel.

###  C.7 Proof of Proposition 4.14

Proposition 4.14 Let Assumption 3.2, 3.3, 4.1, 4.9 hold with
$\bm{\mu}^{\text{E}}\in\mathsf{RL}_{\mathcal{F}}(\bm{r}^{\text{E}})$ for some
$\bm{r}^{\text{E}}\in\mathcal{R}$ and

$\mathcal{R}\subseteq\left\\{{\bm{r}_{\bm{w}}=\bm{\Phi}\bm{w}:\bm{\Phi}\in\mathbb{R}^{mn\times
d},\bm{w}\in\mathbb{R}^{d}}\right\\}.$ (100)

Then, if for
$\bm{\Xi}\vcentcolon=\begin{bmatrix}\bm{E}-\gamma\bm{P},\bm{\Psi}\end{bmatrix}$
it holds that

$\operatorname{\mathbf{rank}}\begin{bmatrix}\bm{\Phi},\bm{\Xi}\end{bmatrix}-\left({\operatorname{\mathbf{rank}}\bm{\Phi}+\operatorname{\mathbf{rank}}\bm{\Xi}}\right)=0,$
(101)

then we have
$\mathsf{IRL}_{\mathcal{R},\mathcal{F}}(\bm{\mu}^{\text{E}})=\left\\{{\bm{r}^{\text{E}}}\right\\}$.

###### Proof.

First, we note that the rank condition is equivalent to the condition that the
subspace spanned by the reward features intersects the Minkowski sum of the
subspace of potential shaping transformations and the subspace spanned by the
safety constraints, only at zero:

$\displaystyle\operatorname{\mathbf{span}}\bm{\Phi}\cap\left({\mathcal{U}+\operatorname{\mathbf{span}}\bm{\Psi}}\right)=\bm{0}$
(102) $\displaystyle\iff$
$\displaystyle\dim\left({\operatorname{\mathbf{span}}\bm{\Phi}\cap\left({\mathcal{U}+\operatorname{\mathbf{span}}\bm{\Psi}}\right)}\right)=0$
$\displaystyle\iff$
$\displaystyle\dim\left({\operatorname{\mathbf{span}}\bm{\Phi}}\right)+\dim\left({\operatorname{\mathbf{span}}\bm{\Xi}}\right)-\dim\left({\operatorname{\mathbf{span}}\begin{bmatrix}\bm{\Phi},\bm{\Xi}\end{bmatrix}}\right)=0$
$\displaystyle\iff$
$\displaystyle\operatorname{\mathbf{rank}}\begin{bmatrix}\bm{\Phi},\bm{\Xi}\end{bmatrix}-\left({\operatorname{\mathbf{rank}}\bm{\Phi}+\operatorname{\mathbf{rank}}\bm{\Xi}}\right)=0.$

To ease notation we define
$\mathcal{V}\vcentcolon=\operatorname{\mathbf{span}}\bm{\Psi}$ and
$\mathcal{W}\vcentcolon=\operatorname{\mathbf{span}}\bm{\Phi}$. Observe that
$\mathcal{U},\mathcal{V},\mathcal{W}$ are linear subspaces of
$\mathbb{R}^{nm}$. By Proposition 4.2, we have
$\bm{r}^{\text{E}}\in\mathsf{IRL}_{\mathcal{R},\mathcal{F}}(\bm{\mu}^{\text{E}})$.
Furthermore, due to (101) it holds
$\left({\mathcal{U}+\mathcal{V}}\right)\cap\mathcal{W}=\bm{0}$. We therefore
get

$\displaystyle\mathsf{IRL}_{\mathcal{R},\mathcal{F}}(\bm{\mu}^{\text{E}})$
$\displaystyle\stackrel{{\scriptstyle(i)}}{{=}}\left({\nabla
f(\bm{\mu}^{\text{E}})+\mathcal{U}+\mathcal{C}}\right)\cap\mathcal{R}$ (103)
$\displaystyle\stackrel{{\scriptstyle(ii)}}{{\subseteq}}\left({\nabla
f(\bm{\mu}^{\text{E}})+\mathcal{U}+\mathcal{C}}\right)\cap\mathcal{W}$
$\displaystyle\stackrel{{\scriptstyle(iii)}}{{\subseteq}}\left({\nabla
f(\bm{\mu}^{\text{E}})+\mathcal{U}+\mathcal{V}}\right)\cap\mathcal{W}$
$\displaystyle\stackrel{{\scriptstyle(iv)}}{{=}}\left({\bm{r}^{\text{E}}+\mathcal{U}+\mathcal{V}}\right)\cap\mathcal{W}$
$\displaystyle\stackrel{{\scriptstyle(v)}}{{=}}\bm{r}^{\text{E}}+\left({\mathcal{U}+\mathcal{V}}\right)\cap\mathcal{W}$
$\displaystyle\stackrel{{\scriptstyle(vi)}}{{=}}\bm{r}^{\text{E}}$

Here, we used Theorem 4.5 in $(i)$, the reward class (100) in $(ii)$, and the
inclusion $(iii)$ holds since $\mathcal{C}\subset\mathcal{V}$. Furthermore,
$(iv)$ and $(v)$ follow from
$\bm{r}^{\text{E}}\in\mathsf{IRL}_{\mathcal{R},\mathcal{F}}(\bm{\mu}^{\text{E}})$,
and $(vi)$ from (101). This concludes the proof. ∎

##  Appendix D Proof of Theorem 5.1

Theorem 5.1 Let Assumption 4.1 and let
$\bm{\mu}^{\text{E}}\in\mathsf{RL}_{\mathcal{F}}(\bm{r}^{\text{E}})$ for some
$\bm{r}^{\text{E}}\in\mathcal{R}\vcentcolon=\mathcal{R}^{\left\lVert\cdot\right\rVert_{1}}$.
Let
$\hat{\bm{\mu}}\in\mathsf{RL}_{\mathcal{F}}\circ\mathsf{IRL}_{\mathcal{R},\mathcal{F}}(\hat{\bm{\mu}}^{\text{E}}_{\mathcal{D}})$
and $R\vcentcolon=\max_{s,a}\left\lVert\bm{\Phi}(s,a)\right\rVert_{\infty}$.
Choosing

$N=\left\lceil\dfrac{32R^{2}}{\varepsilon^{2}}\log\left({\dfrac{2d}{\delta}}\right)\right\rceil\text{
and
}T=\left\lceil\log\left({\dfrac{\varepsilon}{8R}}\right)/\log(\gamma)\right\rceil,$
(104)

it holds with probability at least $1-\delta$

$\displaystyle
J(\bm{\mu}^{\text{E}},\bm{r}^{\text{E}})-J(\hat{\bm{\mu}},\bm{r}^{\text{E}})\leq\varepsilon,$
(105) $\displaystyle
J(\hat{\bm{\mu}},\hat{\bm{r}})-J(\bm{\mu}^{\text{E}},\hat{\bm{r}})\leq\varepsilon,$
$\displaystyle\forall\hat{\bm{r}}\in\mathsf{IRL}_{\mathcal{R},\mathcal{F}}(\hat{\bm{\mu}}^{\text{E}}_{\mathcal{D}}),$

where $J(\bm{\mu},\bm{r})\vcentcolon=\bm{r}^{\top}\bm{\mu}-f(\bm{\mu})$.
Moreover, if

  1. (a)

$f$ is $L$-strongly convex with respect to the norm
$\left\lVert\cdot\right\rVert$, it holds with probability at least $1-\delta$

$\left\lVert\hat{\bm{\mu}}-\bm{\mu}^{\text{E}}\right\rVert\leq\sqrt{\dfrac{2\varepsilon}{L}}.$
(106)

  2. (b)

$f(\bm{\mu})=-\beta\,\mathbb{E}_{(s,a)\sim\bm{\mu}}\left[{H\left(\bm{\pi}^{\bm{\mu}}(\cdot|s)\right)}\right]$
with $\beta>0$, it holds with probability at least $1-\delta$

$\mathbb{E}_{(s,a)\sim\bm{\mu}^{\text{E}}}\left[{\left\lVert\bm{\pi}^{\hat{\bm{\mu}}}(\cdot|s)-\bm{\pi}^{\text{E}}(\cdot|s)\right\rVert_{1}}\right]\leq\sqrt{\dfrac{2\varepsilon}{\beta}}.$
(107)

###### Proof.

Consider the idealized and the empirical min-max objective
$L(\bm{\mu},\bm{w})\vcentcolon=\bm{r}_{\bm{w}}^{\top}\left({\bm{\mu}-\bm{\mu}^{\text{E}}}\right)-f(\bm{\mu})$
and
$\hat{L}(\bm{\mu},\bm{w})\vcentcolon=\bm{r}_{\bm{w}}^{\top}\left({\bm{\mu}-\hat{\bm{\mu}}^{\text{E}}_{\mathcal{D}}}\right)-f(\bm{\mu})$,
the main idea of the proof it to get a uniform bound on
$\left|{L-\hat{L}}\right|$. The statement in (105) is then following from the
saddle point property of $(\bm{\mu}^{\text{E}},\bm{r}^{\text{E}})$ and
$(\hat{\bm{\mu}},\hat{\bm{r}})$. For (106) (and (107)) we are then using
strict concavity (and the soft-suboptimality Lemma) to translate proximity of
the optimal value to proximity of the optimal occupancy (and the optimal
policy).

Step 1:
First, we define the true expert feature expectations
$\bm{\sigma}^{\text{E}}\vcentcolon=\bm{\Phi}^{\top}\bm{\mu}^{\text{E}}$, the
empirical expert feature expectation
$\hat{\bm{\sigma}}^{\text{E}}_{\mathcal{D}}\vcentcolon=\bm{\Phi}^{\top}\hat{\bm{\mu}}^{\text{E}}_{\mathcal{D}}$,
as well as
$\bm{\sigma}^{\text{E}}_{T}\vcentcolon=(1-\gamma)\mathbb{E}\left[{\sum_{t=0}^{T}\gamma^{t}\bm{\Phi}(s_{t},a_{t})\big{|}\pi^{\text{E}}}\right]$.
We can then decompose $\max\left|{\hat{L}-L}\right|$ as follows:

$\displaystyle\max_{\bm{\mu}\in\mathcal{F},\left\lVert\bm{w}\right\rVert_{1}\leq
1}\left|{L(\bm{\mu},\bm{w})-\hat{L}(\bm{\mu},\bm{w})}\right|=\max_{\left\lVert\bm{w}\right\rVert_{1}\leq
1}\left|{\bm{w}^{\top}\left({\hat{\bm{\sigma}}^{\text{E}}_{\mathcal{D}}-\bm{\sigma}^{\text{E}}}\right)}\right|$
(108)
$\displaystyle\stackrel{{\scriptstyle(i)}}{{\leq}}\max_{\left\lVert\bm{w}\right\rVert_{1}\leq
1}\left|{\bm{w}^{\top}\left({\bm{\sigma}^{\text{E}}-\bm{\sigma}^{\text{E}}_{T}}\right)}\right|+\max_{\left\lVert\bm{w}\right\rVert_{1}\leq
1}\left|{\bm{w}^{\top}\left({\bm{\sigma}^{\text{E}}_{T}-\hat{\bm{\sigma}}^{\text{E}}_{\mathcal{D}}}\right)}\right|$
$\displaystyle\stackrel{{\scriptstyle(ii)}}{{\leq}}\underbrace{\left\lVert\bm{\sigma}^{\text{E}}-\bm{\sigma}^{\text{E}}_{T}\right\rVert_{\infty}}_{I_{1}}+\underbrace{\left\lVert\bm{\sigma}^{\text{E}}_{T}-\hat{\bm{\sigma}}^{\text{E}}_{\mathcal{D}}\right\rVert_{\infty}}_{I_{2}}$

Here, $(i)$ follows from the triangle inequality for the supremum norm and
$(ii)$ from Hölder’s inequality
$|\bm{x}^{\top}\bm{y}|\leq\left\lVert\bm{x}\right\rVert_{1}\left\lVert\bm{y}\right\rVert_{\infty}$.
The first term is readily bounded by

$I_{1}=\left\lVert(1-\gamma)\mathbb{E}_{\pi}\sum_{t=T+1}^{\infty}\gamma^{t}\bm{\Phi}(s_{t},a_{t})\right\rVert_{\infty}\leq\gamma^{T+1}R\leq\gamma^{T}R,$
(109)

where
$R\vcentcolon=\max_{s,a}\left\lVert\bm{\Phi}(s,a)\right\rVert_{\infty}=\max_{s,a}\max_{\left\lVert\bm{w}\right\rVert_{1}\leq
1}\left|{\bm{r}_{\bm{w}}(s,a)}\right|$. Thus in order to have
$I_{1}\leq\varepsilon_{1}$ it suffices to choose
$T=\left\lceil\log(\varepsilon_{1}/R)/\log(\gamma)\right\rceil$.

For the second term $I_{2}$, we make use of Hoeffding’s inequality

######  Lemma D.1 (Hoeffding).

Consider iid random variables $X_{1},\ldots,X_{N}$ with $X_{i}\in[a,b]$ and
let $\bar{X}_{N}:=\frac{1}{N}(X_{1}+\ldots+X_{N})$. Then,

$\Pr\left({\left|{\bar{X}_{N}-\mathbb{E}X_{i}}\right|\geq t}\right)\leq
2\exp\left({-\dfrac{2t^{2}N}{(b-a)^{2}}}\right).$ (110)

Since
$\left|{\left({\hat{\bm{\sigma}}^{\text{E}}_{\mathcal{D}}}\right)_{j}}\right|\leq
R$ and
$\mathbb{E}\left({\hat{\bm{\sigma}}^{\text{E}}_{\mathcal{D}}}\right)_{j}=\left({\bm{\sigma}^{\text{E}}_{T}}\right)_{j}$
for all $j=1,\ldots,d$, we get by Hoeffding’s inequality

$\Pr\left({\left|{\left({\hat{\bm{\sigma}}^{\text{E}}_{\mathcal{D}}-\bm{\sigma}^{\text{E}}_{T}}\right)_{j}}\right|\geq\varepsilon_{2}}\right)\leq
2\exp\left({-\dfrac{\varepsilon_{2}^{2}N}{2R^{2}}}\right),\quad j=1,\ldots,d.$
(111)

Using the union bound

$\displaystyle\Pr\left({I_{2}<\epsilon_{2}}\right)$
$\displaystyle=1-\Pr\left({I_{2}\geq\epsilon_{2}}\right)$ (112)
$\displaystyle\geq
1-\left({\Pr\left({\left|{\left({\hat{\bm{\sigma}}^{\text{E}}_{\mathcal{D}}-\bm{\sigma}^{\text{E}}_{T}}\right)_{1}}\right|\geq\varepsilon_{2}}\right)+\ldots+\Pr\left({\left|{\left({\hat{\bm{\sigma}}^{\text{E}}_{\mathcal{D}}-\bm{\sigma}^{\text{E}}_{T}}\right)_{d}}\right|\geq\varepsilon_{2}}\right)}\right),$

yields $I_{2}\leq\varepsilon_{2}$ with probability at least $1-\delta$ when
$N\geq\log\left({2d/\delta}\right)2R^{2}/\varepsilon_{2}^{2}$. Thus choosing
$N=\left\lceil\log\left({2d/\delta}\right)32R^{2}/\varepsilon^{2}\right\rceil$
and $T=\left\lceil\log(\varepsilon/(8R))/\log(\gamma)\right\rceil$ it holds
with probability at least $1-\delta$

$\max_{\bm{\mu}\in\mathcal{F},\left\lVert\bm{w}\right\rVert_{1}\leq
1}\left|{L(\bm{\mu},\bm{w})-\hat{L}(\bm{\mu},\bm{w})}\right|\leq\varepsilon/2.$
(113)

Step 2:
Next we use the fact that for two real-valued functions $g$ and $h$ we always
have $\left|{\max g-\max h}\right|\leq\max\left|{g-h}\right|$ and similarly
$\left|{\min g-\min h}\right|\leq\max\left|{g-h}\right|$. Therefore, it holds
for all
$\hat{\bm{r}}\in\mathsf{IRL}_{\mathcal{R},\mathcal{F}}(\hat{\bm{\mu}}^{\text{E}}_{\mathcal{D}})$
with probability at least $1-\delta$:

$\displaystyle\left|{L(\bm{\mu}^{\text{E}},\bm{r}^{\text{E}})-\hat{L}(\hat{\bm{\mu}},\hat{\bm{r}})}\right|$
$\displaystyle=\left|{\min_{\left\lVert\bm{w}\right\rVert_{1}\leq
1}\max_{\bm{\mu}\in\mathcal{F}}L(\bm{\mu},\bm{w})-\min_{\left\lVert\bm{w}\right\rVert_{1}\leq
1}\max_{\bm{\mu}\in\mathcal{F}}\hat{L}(\bm{\mu},\bm{w})}\right|$ (114)
$\displaystyle\leq\max_{\left\lVert\bm{w}\right\rVert_{1}\leq
1}\left|{\max_{\bm{\mu}\in\mathcal{F}}L(\bm{\mu},\bm{w})-\max_{\bm{\mu}\in\mathcal{F}}\hat{L}(\bm{\mu},\bm{w})}\right|$
$\displaystyle\leq\max_{\bm{\mu}\in\mathcal{F},\left\lVert\bm{w}\right\rVert_{1}\leq
1}\left|{\hat{L}(\bm{\mu},\bm{w})-L(\bm{\mu},\bm{w})}\right|\leq\varepsilon/2.$

Due to Assumption 4.1 there is $\bm{w}^{*}$ with
$\left\lVert\bm{w}^{*}\right\rVert_{1}\leq 1$ such that
$\bm{r}_{\bm{w}^{*}}=\bm{r}^{\text{E}}$. Furthermore, let $\hat{\bm{w}}$ with
$\left\lVert\hat{\bm{w}}\right\rVert_{1}\leq 1$ be such that
$\bm{r}_{\hat{\bm{w}}}\in\mathsf{IRL}_{\mathcal{R},\mathcal{F}}(\hat{\bm{\mu}}^{\text{E}}_{\mathcal{D}})$.
Then, $L$ has a saddle-point in $(\bm{\mu}^{\text{E}},\bm{w}^{*})$ and
$\hat{L}$ in $(\hat{\bm{\mu}},\hat{\bm{w}})$. Choosing $N$ and $T$ as in
(113), it holds with probability at least $1-\delta$

$\displaystyle
J(\bm{\mu}^{\text{E}},\bm{r}^{\text{E}})-J(\hat{\bm{\mu}},\bm{r}^{\text{E}})$
$\displaystyle=L(\bm{\mu}^{\text{E}},\bm{w}^{*})-L(\hat{\bm{\mu}},\bm{w}^{*})$
(115)
$\displaystyle\stackrel{{\scriptstyle(i)}}{{\leq}}L(\bm{\mu}^{\text{E}},\bm{w}^{*})-\hat{L}(\hat{\bm{\mu}},\bm{w}^{*})+\varepsilon/2$
$\displaystyle\stackrel{{\scriptstyle(ii)}}{{\leq}}L(\bm{\mu}^{\text{E}},\bm{r}^{\text{E}})-\hat{L}(\hat{\bm{\mu}},\hat{\bm{w}})+\varepsilon/2$
$\displaystyle\stackrel{{\scriptstyle(iii)}}{{\leq}}\varepsilon,$

where $(i)$ is a consequence of (113), inequality $(ii)$ follows since
$\hat{L}$ has a saddle point in $(\hat{\bm{\mu}},\hat{\bm{w}})$, and $(iii)$
from (114). This proves the first result in (105). The second follows
analogously from a similar series of inequalities

$\displaystyle
J(\hat{\bm{\mu}},\hat{\bm{r}})-J(\bm{\mu}^{\text{E}},\hat{\bm{r}})$
$\displaystyle=L(\hat{\bm{\mu}},\hat{\bm{w}})-L(\bm{\mu}^{\text{E}},\hat{\bm{w}})$
(116)
$\displaystyle\leq\hat{L}(\hat{\bm{\mu}},\hat{\bm{w}})-L(\bm{\mu}^{\text{E}},\hat{\bm{w}})+\varepsilon/2$
$\displaystyle\leq\hat{L}(\hat{\bm{\mu}},\hat{\bm{w}})-L(\bm{\mu}^{\text{E}},\bm{w}^{*})+\varepsilon/2$
$\displaystyle\leq\varepsilon.$

Step 3:
To prove the inequality (106), we use the fact that since $f$ is $L$-strongly
convex,
$J(\bm{\mu},\bm{r}^{\text{E}})={\bm{r}^{\text{E}}}^{\top}\bm{\mu}-f(\bm{\mu}_{k})$
is $L$-strongly concave in $\bm{\mu}$ i.e. it holds

$J(\bm{\mu},\bm{r}^{\text{E}})\leq
J(\bm{\mu}^{\text{E}},\bm{r}^{\text{E}})+\nabla_{\bm{\mu}}J(\bm{\mu}^{\text{E}},\bm{r}^{\text{E}})^{\top}\left({\bm{\mu}-\bm{\mu}^{\text{E}}}\right)-\dfrac{L}{2}\left\lVert\bm{\mu}-\bm{\mu}^{\text{E}}\right\rVert^{2}.$
(117)

By optimality of $\bm{\mu}^{\text{E}}$ it holds
$\nabla_{\bm{\mu}}J(\bm{\mu}^{\text{E}},\bm{r}^{\text{E}})^{\top}\left({\bm{\mu}-\bm{\mu}^{\text{E}}}\right)\leq
0$ for all $\bm{\mu}\in\mathcal{F}$. Rearranging terms yields and taking the
square root yields

$\left\lVert\bm{\mu}-\bm{\mu}^{\text{E}}\right\rVert\leq\sqrt{\dfrac{2}{L}\left({J(\bm{\mu}^{\text{E}},\bm{r}^{\text{E}})-J(\bm{\mu},\bm{r}^{\text{E}})}\right)},\quad\forall\mu\in\mathcal{F}.$
(118)

Combining this with (105) yields the desired result

$\left\lVert\hat{\bm{\mu}}-\bm{\mu}^{\text{E}}\right\rVert\leq\sqrt{\dfrac{2\varepsilon}{L}}.$
(119)

Although
$f(\bm{\mu})=-\beta\,\mathbb{E}_{(s,a)\sim\bm{\mu}}\left[{H\left(\bm{\pi}^{\bm{\mu}}(\cdot|s)\right)}\right]$
is in general not strongly convex (although it is strictly convex), we can
make use of the following result for entropy-regularized MDPs.

######  Lemma D.2 ((Mei et al., 2020)).

For any occupancy measure $\bm{\mu}\in\mathcal{F}$ it holds

$J(\bm{\mu}^{*},\bm{r})-J(\bm{\mu},\bm{r})=\beta\sum_{s}\bm{\nu}(s)D_{\text{KL}}\left({\bm{\pi}^{\bm{\mu}}(\cdot|s)||\bm{\pi}^{\bm{\mu}^{*}}(\cdot|s)}\right),$

where $\bm{\mu}^{*}\in\mathsf{RL}_{\mathcal{F}}(\bm{r})$ and
$\bm{\nu}(s)=\sum_{a}\bm{\mu}(s,a)$ is the state occupancy measure.

Making use of the following lower bound on the KL-divergence (Cover, 1999)

$D_{\text{KL}}(\bm{q}||\bm{p})\leq\dfrac{1}{2}\left\lVert\bm{q}-\bm{p}\right\rVert_{1}^{2},$
(120)

we arrive at

$\displaystyle
J(\hat{\bm{\mu}},\hat{\bm{r}})-J(\bm{\mu}^{\text{E}},\hat{\bm{r}})$
$\displaystyle\geq\beta\sum_{s}\bm{\nu}^{\text{E}}(s)D_{\text{KL}}\left({\bm{\pi}^{\bm{\mu}^{\text{E}}}(\cdot|s)||\bm{\pi}^{\hat{\bm{\mu}}}(\cdot|s)}\right)$
$\displaystyle\geq\dfrac{\beta}{2}\sum_{s}\bm{\nu}^{\text{E}}(s)\left\lVert\bm{\pi}^{\bm{\mu}^{\text{E}}}(\cdot|s)-\bm{\pi}^{\hat{\bm{\mu}}}(\cdot|s)\right\rVert_{1}^{2}$
$\displaystyle\geq\dfrac{\beta}{2}\left({\sum_{s}\bm{\nu}^{\text{E}}(s)\left\lVert\bm{\pi}^{\bm{\mu}^{\text{E}}}(\cdot|s)-\bm{\pi}^{\hat{\bm{\mu}}}(\cdot|s)\right\rVert_{1}}\right)^{2}$
$\displaystyle=\dfrac{\beta}{2}\left({\mathbb{E}_{(s,a)\sim\bm{\mu}^{\text{E}}}\left[{\left\lVert\bm{\pi}^{\bm{\mu}^{\text{E}}}(\cdot|s)-\bm{\pi}^{\hat{\bm{\mu}}}(\cdot|s)\right\rVert_{1}}\right]}\right)^{2},$
(121)

where the last inequality follows from Jensen’s inequality. Making use of
(105) and rearranging terms yields

$\mathbb{E}_{(s,a)\sim\bm{\mu}^{\text{E}}}\left[{\left\lVert\bm{\pi}^{\bm{\mu}^{\text{E}}}(\cdot|s)-\bm{\pi}^{\hat{\bm{\mu}}}(\cdot|s)\right\rVert_{1}}\right]\leq\sqrt{\dfrac{2\varepsilon}{\beta}}.$
(122)

∎

##  Appendix E Algorithm

We present the policy based algorithm for entropy regularization as used in
the experiments. To this end, recall the min-max problem (23)

$\min_{\bm{r}\in\mathcal{R}}\max_{\bm{\mu}\in\mathcal{F}}\;\;\bm{r}^{\top}\left({\bm{\mu}-\hat{\bm{\mu}}^{\text{E}}_{\mathcal{D}}}\right)-f(\bm{\mu}),$
(123)

and the entropy regularization
$f(\bm{\mu})=-\beta\,\mathbb{E}_{(s,a)\sim\bm{\mu}}\left[{H\left(\bm{\pi}^{\bm{\mu}}(\cdot|s)\right)}\right]$.
Applying Proposition 3.4 this is equivalent to

$\displaystyle\min_{\bm{\xi}\geq\bm{0},\bm{r}\in\mathcal{R}}\max_{\bm{\mu}\in\mathcal{M}}\bm{r}^{\top}\left({\bm{\mu}-\hat{\bm{\mu}}^{\text{E}}_{\mathcal{D}}}\right)-f(\bm{\mu})+\bm{\xi}^{\top}\left({\bm{b}-\bm{\Psi}^{\top}\bm{\mu}}\right).$
(124)

Using the one-to-one mapping between policies and occupancy measures and the
linear reward class

$\displaystyle\mathcal{R}\vcentcolon=\left\\{{\bm{r}_{\bm{w}}=\bm{\Phi}\bm{w}:\bm{\Phi}\in\mathbb{R}^{nm\times
d},\left\lVert\bm{w}\right\rVert\leq c}\right\\},$ (125)

this can be rewritten as

$\displaystyle\min_{\bm{\xi}\geq\bm{0},\bm{r}\in\mathcal{R}}\max_{\bm{\pi}\in\Pi}\bm{r}^{\top}\left({\bm{\mu}^{\bm{\pi}}-\hat{\bm{\mu}}^{\text{E}}_{\mathcal{D}}}\right)-f(\bm{\mu}^{\bm{\pi}})+\bm{\xi}^{\top}\left({\bm{b}-\bm{\Psi}^{\top}\bm{\mu}^{\bm{\pi}}}\right)$
(126) $\displaystyle=$
$\displaystyle\min_{\bm{\xi}\geq\bm{0},\left\lVert\bm{w}\right\rVert\leq
c}\max_{\bm{\pi}\in\Pi}\bm{w}^{\top}\bm{\Phi}^{\top}\left({\bm{\mu}^{\bm{\pi}}-\hat{\bm{\mu}}^{\text{E}}_{\mathcal{D}}}\right)-f(\bm{\mu}^{\bm{\pi}})+\bm{\xi}^{\top}\left({\bm{b}-\bm{\Psi}^{\top}\bm{\mu}^{\bm{\pi}}}\right)$
$\displaystyle=$
$\displaystyle\min_{\bm{\xi}\geq\bm{0},\left\lVert\bm{w}\right\rVert\leq
c}\max_{\bm{\pi}\in\Pi}L_{\mathcal{D}}(\bm{\pi},\bm{w},\bm{\xi}),$

with the Lagrangian
$L_{\mathcal{D}}(\bm{\pi},\bm{w},\bm{\xi})\vcentcolon=\bm{w}^{\top}\bm{\Phi}^{\top}\left({\bm{\mu}^{\bm{\pi}}-\hat{\bm{\mu}}^{\text{E}}_{\mathcal{D}}}\right)-f(\bm{\mu}^{\bm{\pi}})+\bm{\xi}^{\top}\left({\bm{b}-\bm{\Psi}^{\top}\bm{\mu}^{\bm{\pi}}}\right)$.
Motivated by recent advances in min-max optimization (Daskalakis & Panageas,
2018), we suggest to use a gradient descent-ascent method, where policy and
reward are updated simultaneously within a single optimization loop.

Algorithm 1 Gradient Descent Ascent for Constrained Entropy-Regularized IRL

Input: Expert data $\mathcal{D}$, learning rate $\eta$.

Initialize $\bm{\pi}\in\Pi,\bm{w}=\bm{0},\bm{\xi}=\bm{0}$.

for $i=1$ to $N_{\text{episodes}}$ do

$\bm{r}\leftarrow\bm{\Phi}\bm{w}-\bm{\Psi}\bm{\xi}$

$\bm{\pi}\leftarrow$ NPG$(\bm{\pi},\bm{r},\eta)$

$\bm{w}\leftarrow
P_{B_{c}}\left({\bm{w}-\eta\nabla_{\bm{w}}L_{\mathcal{D}}(\bm{\pi},\bm{w},\bm{\xi})}\right)$

$\bm{\xi}\leftarrow
P_{[\bm{0},\bm{\infty})}\left({\bm{\xi}-\eta\nabla_{\bm{\xi}}L_{\mathcal{D}}(\bm{\pi},\bm{w},\bm{\xi})}\right)$

end for

Return: $\bm{w},\bm{\pi}$.

Here, NPG$(\bm{\pi},\bm{r},\eta)$ refers to a single entropy-regularized
Natural Policy Gradient step with softmax parametrization of the policy (Cen
et al., 2022). Note that for the step size $\eta=(1-\gamma)/\beta$ this policy
gradient step is completely equivalent to soft policy iteration (Haarnoja et
al., 2018). Furthermore, the gradients for $\bm{w}$ and $\bm{\xi}$ are given
by

$\displaystyle\nabla_{\bm{w}}L_{\mathcal{D}}(\bm{\pi},\bm{w},\bm{\xi})$
$\displaystyle=\bm{\Phi}^{\top}\left({\bm{\mu}^{\bm{\pi}}-\hat{\bm{\mu}}^{\text{E}}_{\mathcal{D}}}\right)$
(127)
$\displaystyle\nabla_{\bm{\xi}}L_{\mathcal{D}}(\bm{\pi},\bm{w},\bm{\xi})$
$\displaystyle=\left({\bm{b}-\bm{\Psi}^{\top}\bm{\mu}^{\bm{\pi}}}\right).$

Moreover, $P_{B_{c}}$ denotes the projection onto the ball
$B_{c}\vcentcolon=\left\\{{\bm{w}:\left\lVert\bm{w}\right\rVert\leq
c}\right\\}$ and $P_{[\bm{0},\bm{\infty})}$ the trivial projection onto the
non-negative orthant. For the $1$-norm projection we use the efficient
projection algorithm by Duchi et al. (2008) and its implementation by Ong &
Lustig (2019). Algorithm 1 has been shown to provably converge for CMDPs (Ying
et al., 2022; Ding et al., 2022) (with known reward) and IRL (Zeng et al.,
2022) (without safety constraints).

Finally, we note that in our code we also provide a gradient descent ascent
primal dual method that directly optimizes the min-max problem (124) in the
occupancy measure space. Since the problem is convex-concave such algorithms
enjoy provable convergence guarantees (Daskalakis & Panageas, 2018;
Nemirovski, 2004), and work for arbitrary convex regularizations $f$. As the
focus of this paper is not on algorithmic convergence and Algorithm 1 was both
– easier to tune and converged more efficiently – we used Algorithm 1
throughout the presented experiments.