# Polaritonic ultracold reactions: cavity controlled molecular
photoassociation

Vasil Rokaj<EMAIL_ADDRESS>ITAMP, Center for Astrophysics $|$
Harvard $\&$ Smithsonian, Cambridge, USA  Department of Physics, Harvard
University, Cambridge, USA  Simeon I. Mistakidis  ITAMP, Center for
Astrophysics $|$ Harvard $\&$ Smithsonian, Cambridge, USA  Department of
Physics, Harvard University, Cambridge, USA  H. R. Sadeghpour  ITAMP, Center
for Astrophysics $|$ Harvard $\&$ Smithsonian, Cambridge, USA

###### Abstract

We introduce a prototypical model for cavity polaritonic control of ultracold
photochemistry by considering the resonant vibrational strong coupling of a
rubidium dimer ($\rm{Rb_{2}}$) to a terahertz cavity. We demonstrate that at
avoided crossings between a vibrational excitation and the vacuum photon
absorption, the resulting polaritonic states between the molecule and photons
can efficiently control the molecular vibrational Franck-Condon (FC) factors.
Due to the entanglement between light and matter, FC factor is transferred
from one polaritonic branch to the other, leading to a polariton with a
substantially enhanced FC factor. Utilizing this polariton state for
photoassociation results in the enhanced formation of ultracold molecules.
This work suggests a path to controlling photoassociation with cavity vacuum
fields, and lays the ground for the emerging subfield of polaritonic ultracold
chemistry.

###### pacs:

Valid PACS appear here

Achieving two of the holy grails of chemistry Kohler  _et al._ (1995),
manipulation of matter at the atomic level Eigler and Schweizer (1990);
Sugimoto  _et al._ (2005); Gross  _et al._ (2009), and coherent control over
chemical pathways Tannor  _et al._ (1986); Shapiro and Brumer (2003); Peirce
_et al._ (1988) has made great strides in the last few decades. Meanwhile,
quantum leap advances in cooling and trapping of atoms and molecules are
bringing these dreams closer to reality, as single atoms and molecules can be
now trapped in periodic arrays, created by optical lattices or tweezers
Anderegg  _et al._ (2019); Kaufman and Ni (2012); Anderegg  _et al._ (2018);
Park  _et al._ (2023). These advances have enabled the application of
ultracold molecules in quantum chemistry Ni  _et al._ (2010), searches for
fundamental physics Safronova  _et al._ (2018), quantum simulation DeMille
(2002); Blackmore  _et al._ (2018); Anderegg  _et al._ (2019), and sensing
Carr  _et al._ (2009); Quéméner and Julienne (2012). This wide range of
applications renders the production of ultracold molecules highly desirable.
Typically, for the enhancement of molecular photoassociation the stimulated
Raman adiabatic passage (STIRAP) is employed, which consists of a so-called
counterintuitive succession of two laser pulses driving molecular transitions
Winkler  _et al._ (2007); Vitanov  _et al._ (2017) as it is illustrated in
Fig. 1(c).

Recently, the possibility to control chemical reactions Hutchison  _et al._
(2012); Ebbesen (2016); Hutchison  _et al._ (2013); Orgiu  _et al._ (2015);
Galego  _et al._ (2016); Flick  _et al._ (2017); Schäfer  _et al._ (2019) and
material properties Bloch  _et al._ (2022); Schlawin  _et al._ (2022);
Paravicini-Bagliani  _et al._ (2019); Appugliese  _et al._ (2022); Rokaj  _et
al._ (2022, 2023a) with cavity vacuum fields has emerged. The experimental
demonstration of vibrational strong coupling and resonant modification of
reaction rates Hutchison  _et al._ (2012); Ahn  _et al._ (2023) indicates that
chemistry can be affected by the quantum hybridization of molecular states and
cavity photons, known as polaritons Ruggenthaler  _et al._ ; Campos-Gonzalez-
Angulo  _et al._ (2023); Mandal  _et al._ (2023); Garcia-Vidal  _et al._
(2021). The polaritons formed in vacuum cavities are free from external
driving fields that can lead to thermalization or destroy quantum coherence
Rider and Barnes (2021).

It may now be possible to combine some of the unique elements of ultracold
molecular physics with cavity quantum electrodynamics (cQED) to steer
photochemistry through vibrational polariton states. In this Letter, we
present a prototypical model of an ultracold reaction under strong coupling to
the cavity field. We showcase that vibrational polariton states, formed
between the $\textrm{Rb}_{2}$ dimer and the cavity photons, can be utilized as
a new paradigm for the controlled photoassociation of molecules.

Figure 1: a) Schematic depiction of the $\textrm{Rb}_{2}$ molecule coupled to
a cavity mode. The distance between the cavity mirrors is $L$ and the area of
the mirrors $A$. b) Potential energy surfaces of $\textrm{Rb}_{2}$ for the
$\Sigma^{+}_{u},\Pi_{u}$ and $\Sigma^{+}_{g}$ states. c) Gaussian pulses
applied from the open sides of the cavity for the photoassociation of
$\textrm{Rb}_{2}$ molecules to the vibrational ground state in
$\Sigma^{+}_{g}$.

We find that at the avoided crossings, the transition probabilities between
the vibrational states, i.e., the FC factors, of the polaritons are
significantly modified. More precisely, in the vicinity of an avoided
crossing, one polariton branch acquires an enhanced FC factor while the FC
factor for the other polariton branch is suppressed. We name this phenomenon
polariton FC transfer, as the population is transferred from one polariton
state to the other. We show that the polariton FC transfer is a result of
entanglement between the photon and molecule excitations. The FC factors
determine the molecular yield. By utilizing the polariton branch with the
enhanced FC factors for STIRAP, it is shown that molecular photoassociation of
$\textrm{Rb}_{2}$ is substantially enhanced when the molecules are coupled to
the cavity field.

Rubidium Dimer in the Cavity.—We consider the $\textrm{Rb}_{2}$ coupled to a
single-mode terahertz cavity Hutchison  _et al._ (2012); Ebbesen (2016);
Hutchison  _et al._ (2013); Ahn  _et al._ (2023). The molecule-photon
interacting system is described by the length gauge Hamiltonian Rokaj  _et
al._ (2018); Cohen-Tannoudji  _et al._ (1997)

$\displaystyle\mathcal{H}=\mathcal{H}_{M}+\hbar\omega
a^{\dagger}a-\sqrt{\frac{\hbar\omega}{2}}\lambda
d_{R}\left(a+a^{\dagger}\right)+\frac{\lambda^{2}d^{2}_{R}}{2}.$ (1)

Here, $\mathcal{H}_{M}=-\frac{\hbar^{2}}{2\mu}\frac{\partial^{2}}{\partial
R^{2}}+D\left[e^{-2w(R-R_{0})}-2e^{-w(R-R_{0})}\right]$ represents the
molecular Hamiltonian of $\textrm{Rb}_{2}$. All potential energy surfaces
(PESs) of $\textrm{Rb}_{2}$ are described with the Morse potential Morse
(1929) characterized by the experimental parameters reported in Ref. Park  _et
al._ (2001), thus rendering our molecular model quantitatively precise. The
reduced mass of the dimer is $\mu$, $D$ ($w$) denotes the depth (width) of the
potential well, and $R_{0}$ is the equilibrium bond distance. Further, we
assume a linear dipole operator $d_{R}=\mu_{0}R$ with $\mu_{0}$ being the
strength of the molecular dipole moment. The length gauge Hamiltonian contains
both the bilinear coupling, and the quadratic dipole self-energy,
$\sim\lambda^{2}d^{2}_{R}$, which is necessary for the stability of the
molecular system and for preserving gauge invariance Rokaj  _et al._ (2018);
De Bernardis  _et al._ (2018). We note that in Eq. 1 the trivial zero-point
energy of the cavity mode has been removed. The cavity field strength
$\lambda=\sqrt{1/\epsilon_{0}\mathcal{V}}$ depends on the cavity volume
$\mathcal{V}=AL=A\pi c/\omega$ with $A$ being the area of the cavity mirrors
having a distance $L$ and $\omega=\pi c/L$ is the frequency of the cavity
mode. Further, $c$ is the speed of light, $\epsilon_{0}$ the vacuum
permittivity, and $a,a^{\dagger}$ are bosonic operators of the cavity mode
satisfying $[a,a^{\dagger}]=1$. The independent cavity-related parameters are
the frequency $\omega$ and the light-matter coupling strength
$g=\mu_{0}\lambda/\sqrt{\omega}=\mu_{0}/\sqrt{A\pi c\epsilon_{0}}$ (see the
Supplementary Material (SM) for more details sup ).

The $\textrm{Rb}_{2}$ molecule is homonuclear, and thus in the
$\Sigma^{+}_{g}$ and the $\Sigma^{+}_{u}$ PESs it has no permanent dipole.
Consequently, the cavity field does not couple the vibrational states in these
PESs, but rather couples only with the $\Pi_{u}$ PES where $\textrm{Rb}_{2}$
has a dipole moment. To treat the photon-molecule system we expand the
Hamiltonian in the basis consisting of the vibrational eigenstates of the
$\Pi_{u}$ PES $\\{|\Phi^{\Pi_{u}}_{n}\rangle\\}$ and the eigenstates
$\\{|i\rangle\\}$ of the bosonic cavity mode Griffiths (1995). As such, we
obtain the matrix form $\mathcal{H}^{ji}_{mn}:=\langle
j|\langle\Phi^{\Pi_{u}}_{m}|\mathcal{H}|\Phi^{\Pi_{u}}_{n}\rangle|i\rangle$ of
the Hamiltonian $\mathcal{H}$. Diagonalizing numerically
$\mathcal{H}^{ji}_{mn}$ we find the vibrational polaritonic energy levels and
the hybrid polaritonic states. For convenience, the Hamiltonian $\mathcal{H}$
is set up in atomic units (a.u.), but then all observables are provided in
units of Hertz. We restrict ourselves to weak and strong light-matter coupling
strengths quantified by $g=10^{-3}\textrm{a.u.}$ and $g=10^{-2}\textrm{a.u.}$
respectively. The corresponding polariton energy spectra, for these $g$
values, are shown in Fig.2 as a function of the cavity frequency for the
vibrational levels $n=60$ and $m=61$ in the $\Pi_{u}$ PES. For both coupling
strengths at the point of resonance indicated by the vertical dashed line in
Figs.2(a) and (b), an avoided crossing emerges which manifests the
hybridization between light and matter. For weak (strong) coupling, the lowest
energetically avoided crossing presented in the inset of Fig. 2(a) (Fig. 2(b))
features the normalized Rabi splitting
$\Omega_{R}=(E_{+}-E_{-})/\hbar\omega_{\rm res}=0.0018$ ($\Omega_{R}=0.018$).
This enhanced normalized Rabi splitting, $\Omega_{R}=0.018$, indicates the
entrance to the strong coupling regime, at which the frequency window where
the avoided crossing occurs becomes larger. This implies that light-matter
entanglement persists for a wider region around the resonance point as can be
readily seen by comparing the insets in Figs. 2(a) and (b).

Figure 2: Polaritonic energy spectra for the vibrational states $n=60$ and
$m=61$ in the $\Pi_{u}$ PES coupled to a cavity mode. In a) the light-matter
coupling strength is $g=0.001\textrm{a.u.}$ and three photonic Fock states are
needed to converge the first avoided crossing (inset) which exhibits
normalized Rabi splitting $\Omega_{R}=1.8\times 10^{-3}$. In b)
$g=0.01\textrm{a.u.}$ is enhanced and 5 photonic Fock states are required to
converge the first avoided crossing (inset) which exhibits an order of
magnitude larger normalized Rabi splitting $\Omega_{R}=1.8\times 10^{-2}$.

Polariton Franck-Condon Transfer.—The FC factors are decisive for the
photoassociation of ultracold molecules, as they determine the transition
probabilities between vibrational states. The polariton states in the
$\Pi_{u}$ PES have the form
$|\Psi_{p}\rangle=\sum_{n,i}c^{p}_{ni}|\Phi^{\Pi_{u}}_{n}\rangle|i\rangle$
with the coefficients $c^{p}_{ni}$ obtained from the diagonalization of
$\mathcal{H}^{ji}_{mn}$. To compute the FC factors between the vibrational
polariton states $|\Psi_{p}\rangle$ in the $\Pi_{u}$ PES, and the vibrational
states in the $\Sigma^{+}_{g}$ and $\Sigma^{+}_{u}$ PESs we need to construct
the tensor product state between the vibrational states
$|\Phi^{\Sigma^{+}_{g,u}}_{v}\rangle$ and the photonic unit vector
$|\mathbb{I}\rangle=\sum^{l}_{i=0}|i\rangle$. Here, $l$ denotes the amount of
photonic Fock states included in the computation of the vibrational
polaritons. Then, the overlap
$\langle\Phi^{\Sigma^{+}_{g,u}}_{v}|\Psi_{p}\rangle:=\langle\mathbb{I}|\langle\Phi^{\Sigma^{+}_{g,u}}_{v}|\Psi_{p}\rangle$
provides the polariton FC factors. Fig. 3 depicts the bare, $\Omega_{R}=0$,
and the respective polaritonic FC factors at the lowest avoided crossing for
$\Omega_{R}=0.0018$ and $\Omega_{R}=0.018$ with respect to the vibrational
levels $n=60$ and $m=61$ in $\Pi_{u}$. Two particular transitions are shown,
from the last vibrational state in $\Sigma^{+}_{u}$ PES,
$\Phi^{\Sigma^{+}_{u}}_{33}$, and to the vibrational ground state
$\Phi^{\Sigma^{+}_{g}}_{0}$ as a function of the cavity frequency. These
states are of specific interest since they will be deployed later for the
STIRAP photoassociation. In both panels of Fig. 3 we observe that before the
resonance point (indicated by the vertical dashed line) the upper (lower)
polariton has FC factor equal to the bare vibrational states $m=61$ ($n=60$).
However, after the resonance point $\omega/2\pi=1.3\textrm{THz}$ the situation
is inverted. The most interesting phenomenon is that at resonance the FC
factor of the lower polariton increases substantially while the one of the
upper polariton is suppressed. Therefore, through the hybridization and
entanglement between the vibrational states with the cavity photons, FC factor
is transferred from the upper polariton to the lower one. This phenomenon
occurring in the vicinity of the avoided crossing will be named in the
following as polariton Franck-Condon transfer. It is crucial for the
photoassociation of molecules because the FC factors determine the molecular
yield. To the best of our knowledge the polariton Franck-Condon transfer has
not been reported elsewhere. Importantly, in order to guarrantee that the
molecular formation through STIRAP using the polariton states will be
enhanced, the FC factors for both transitions,
$\Sigma^{+}_{u}\leftrightarrow\Pi_{u}$ and
$\Pi_{u}\leftrightarrow\Sigma^{+}_{g}$, of the same polariton branch need to
be enhanced.

Insights from the Jaynes-Cummings model.—To shed further light on the behavior
of the polariton FC factors from an analytical perspective, we employ the
Jaynes-Cummings (JC) model. The latter captures the hybridization of two
vibrational levels by a photon in the weak coupling regime, $\Omega_{R}\sim
10^{-3}$, Shore and Knight (1993). However, for strong coupling,
$\Omega_{R}\sim 10^{-2}$, the interactions beyond the rotating-wave
approximation become substantial and the JC model is no longer quantitatively
precise, see SM for further details sup . The hybrid excited states in the JC
model for single photon excitation exactly at resonance are
$|\pm\rangle=(|\Phi^{\Pi_{u}}_{n}\rangle|1\rangle\mp|\Phi^{\Pi_{u}}_{n+1}\rangle|0\rangle)/\sqrt{2}$
Scully and Zubairy (1997); Vogel and Welsch (2006). Their FC factors with the
vibrational states $|\Phi^{\Sigma^{+}_{g,u}}_{v}\rangle$ are defined as
$\textrm{FC}_{\pm}=\langle\mathbb{I}|\langle\Phi^{\Sigma^{+}_{g,u}}_{v}|\pm\rangle$
and can be expressed as

$\textrm{FC}_{\pm}=\frac{1}{\sqrt{2}}\left[\langle\Phi^{\Sigma^{+}_{g,u}}_{v}|\Phi^{\Pi_{u}}_{n}\rangle\mp\langle\Phi^{\Sigma^{+}_{g,u}}_{v}|\Phi^{\Pi_{u}}_{n+1}\rangle\right].$
(2)

It becomes evident that at resonance the FC factors of the JC states can be
significantly modified. Namely, if the vibrational overlaps
$\langle\Phi^{\Sigma^{+}_{g,u}}_{v}|\Phi^{\Pi_{u}}_{n}\rangle,\langle\Phi^{\Sigma^{+}_{g,u}}_{v}|\Phi^{\Pi_{u}}_{n+1}\rangle$
have the same sign then $\textrm{FC}_{-}$ ($\textrm{FC}_{+}$) gets enhanced
(suppressed). However, if the integrals have opposite signs the situation is
inverted. This demonstrates that around the resonance point there is a
significant transfer of FC factor from one branch of the avoided crossing to
the other. The phenomenon of FC factor exchange is a result of entanglement
between the photons and the molecule around the resonance point. Importantly,
for the STIRAP through the lower (upper) polariton to be enhanced it is
necessary the overlap integrals for both transitions to have the same
(opposite) sign. For the states $\Phi^{\Pi_{u}}_{60},\Phi^{\Pi_{u}}_{61}$ the
overlaps for the upwards (downwards) transition are both positive (negative)
[see also SM sup ]. This leads to enhanced FC factor for $\Psi_{-}$ in both
transitions.

Figure 3: Bare and polaritonic FC factors of $\rm{Rb}_{2}$ for normalized Rabi
splitting $\Omega_{R}=1.8\times 10^{-3}$ and $\Omega_{R}=1.8\times 10^{2}$. In
a) the FC factors for the upward transition from $\Phi^{\Sigma^{+}_{u}}_{33}$
to $\Phi^{\Pi_{u}}_{60,61}$ and the respective polariton states
$|\Psi_{\pm}\rangle$ in $\Pi_{u}$ are shown (see legend). In b) the FC factors
for the downward transition to the ground state $\Phi^{\Sigma^{+}_{g}}_{0}$,
starting from $\Phi^{\Pi_{u}}_{60,61}$ and the polariton states
$|\Psi_{\pm}\rangle$ are presented (see legend). As it can be seen, at
resonance the lower (upper) polariton $\Psi_{-}$ ($\Psi_{+}$) shows enhanced
(suppressed) FC factors. Importantly, for strong coupling,
$\Omega_{R}=1.8\times 10^{-2}$, polariton FC transfer is enhanced.

Photoassociation in the Cavity.—The phenomenon of FC transfer around the
avoided crossing and the selective enhancement of polariton states suggests
that molecular photoassociation of $\textrm{Rb}_{2}$ molecules to the
$\Sigma^{+}_{g}$ PES can be controlled with resonant light-molecule coupling.
To demonstrate this we perform a prototype STIRAP photoassociation for
$\textrm{Rb}_{2}$ in the cavity, employing the $\Lambda$-transition scheme
with Gaussian laser fields shown in Fig. 1(c). The first laser drives the
transition $\Sigma^{+}_{u}\leftrightarrow\Pi^{+}_{u}$, and the second one the
$\Pi_{u}\leftrightarrow\Sigma^{+}_{g}$. The populations of each PES during the
STIRAP obey the following coupled differential equations Winkler  _et al._
(2007),

$\displaystyle\textrm{i}\dot{\Sigma}_{u}(t)$ $\displaystyle=$
$\displaystyle-\textrm{i}\frac{\gamma_{\Sigma_{u}}}{2}\Sigma_{u}(t)-\frac{\Omega_{1}}{2}\Pi_{u}(t),$
$\displaystyle\textrm{i}\dot{\Pi}_{u}(t)$ $\displaystyle=$
$\displaystyle-\textrm{i}\frac{\gamma_{\Pi_{u}}}{2}\Pi_{u}(t)-\frac{1}{2}\left(\Omega_{1}\Sigma_{u}(t)+\Omega_{2}\Sigma_{g}(t)\right),$
(3) $\displaystyle\textrm{i}\dot{\Sigma}_{g}(t)$ $\displaystyle=$
$\displaystyle-\textrm{i}\frac{\gamma_{\Sigma_{g}}}{2}\Sigma_{g}(t)-\frac{\Omega_{2}}{2}\Pi_{u}(t).$

Here, $\Omega_{1,2}$ are the Rabi frequencies corresponding to the two laser
fields $\mathbf{E}_{1,2}$ respectively. Further, we use the decay rates
$\gamma_{\Sigma_{u}}=2\pi\times 0.72\;\rm{kHz}$, $\gamma_{\Pi_{u}}=2\pi\times
12\;\rm{MHz}$ for the $\Pi_{u}$ and the $\Sigma_{u}$ states respectively taken
from the experiment of Ref. Winkler  _et al._ (2007). Notice that
$\gamma_{\Sigma_{g}}=0$, because we aim at the vibrational ground state and
thus there is no lower state for the molecule to decay to. Also, in accordance
with the theoretical modeling of the experiment of Ref. Winkler  _et al._
(2007), the detunings of the lasers are assumed to be zero. Within the Born-
Oppenheimer approximation the Rabi frequencies $\Omega_{1,2}$ for the laser
fields $\mathbf{E}_{1,2}$ are a product of the FC factors of the vibrational
levels and the electronic dipole moment for the transition $5s\rightarrow 5p$.
Therefore, it holds that
$\Omega_{i}=\frac{|\mathbf{E}_{i}|}{\hbar}\langle\Phi_{m}|\Phi_{n}\rangle\langle
5s|d_{e}|5p\rangle$. The value for the electronic transition $\langle
5s|d_{e}|5p\rangle$ is obtained from the $A_{5s5p}=3.812\times
10^{7}\textrm{s}^{-1}$ Einstein coefficient Scully and Zubairy (1997) provided
in Ref. Kramida  _et al._ (2022). The decisive difference for the STIRAP in
the cavity is that the bare FC factors of $\textrm{Rb}_{2}$ are replaced by
the FC factors of the polariton states.

The time-evolution of the populations of each state following STIRAP are
illustrated in Figs. 4(a) and (b) for the lower polariton branch for weak
($\Omega_{R}=1.8\times 10^{-3}$) and strong ($\Omega_{R}=1.8\times 10^{-2}$)
coupling respectively. In both cases, as expected after the application of the
second pulse, around $t=1ms$, the system attains a steady state where the
population of the molecular ground state $\Phi^{\Sigma^{+}_{g}}_{0}$ becomes
stationary because the system is no longer driven. Most importantly, we find
that the population of molecules in the vibrational ground state,
$\Phi^{\Sigma^{+}_{g}}_{0}$, is substantially enhanced when STIRAP is
performed through the lower polariton $\Psi_{-}$. Particularly, for weak
coupling $\Omega_{R}=1.8\times 10^{-3}$ photoassociation of molecules becomes
$17.3\%$ and for strong coupling $\Omega_{R}=1.8\times 10^{-2}$ reaches
$18.7\%$. In comparison to molecular photoassociation through the bare
vibrational states $\Phi^{\Pi_{u}}_{60}$ and $\Phi^{\Pi_{u}}_{61}$ which are
$12.4\%$ and $3.9\%$ respectively (see the dashed lines in Fig. 4) we deduce
that with respect to the state $\Phi^{\Pi_{u}}_{60}$ the relative percentage
increase is $39.5\%$ ($50.8\%$) for $\Omega_{R}=1.8\times 10^{-3}$
($\Omega_{R}=1.8\times 10^{-2}$). This is a consequence of the enhanced FC
factors of the lower polariton state shown in Fig. 3. For completeness we
mention that the molecular photoassociation through the upper polariton does
not exhibit any enhancement but it is extremely suppressed, see more details
in SM sup . This happens because the FC factors of $\Psi_{+}$ decrease
drastically at the resonance [see Figs. 3(a) and (b)].

Figure 4: STIRAP photoassociation to $\Phi^{\Sigma^{+}_{g}}_{0}$ through the
lower polariton state $\Psi_{-}$ at the exact resonance point. In a) for weak
coupling $\Omega_{R}=1.8\times 10^{-3}$ the photoassociation is $17.3\%$ and
in b) for strong coupling $\Omega_{R}=1.8\times 10^{-2}$ reaches $18.7\%$. As
compared to photoassociation through the molecular states
$\Phi^{\Pi_{u}}_{60}$ (green dashed line) and the states $\Phi^{\Pi_{u}}_{61}$
(brown dashed line) the lower polariton demonstrates enhanced molecular
photoassociation. This is a consequence of the enhanced polariton FC factors.

Discussion and Outlook.—The photoassociation of ultracold molecules can be
substantially enhanced by resonant vibrational strong coupling to the vacuum
cavity. This is due to the phenomenon of polariton FC transfer that emerges
around the avoided crossings between the polaritons. This generic mechanism,
arising from the photon-molecule entanglement around a resonance, can be
understood within the standard JC paradigm. In the strong coupling regime,
however, the JC model is not adequate and the full length gauge Hamiltonian is
necessary for the accurate description of the photon-molecule interaction. The
relative percentage increase of molecular photoassociation is substantial, and
in the strong coupling regime, it approaches $\sim 50\%$.

It is important to mention however, that the photoassociation through
polaritonic states is not always advantageous for all vibrational states. The
polaritonic STIRAP demonstrates clear advantage only when the FC factors for
both branches in the same polariton state are enhanced in unison in both
directions of the $\Lambda$-scheme. Thus, observing an avoided crossing does
not directly imply that a chemical process of interest will be necessarily
modified; for a more detailed discussion, see the SM sup .

This work introduces the element of cavity polaritonic control to ultracold
chemistry, and paves the way for future investigations at the intersection
between cQED and ultracold chemistry. A particularly exciting research
direction is to couple dipolar ultracold molecules, e.g. $\rm{KRb}$, to cavity
fields where due to the existence of a permanent dipole moment, the ground
state of the molecule can be modified through polariton formation Rokaj  _et
al._ (2023b). In addition, $\rm{KRb-KRb}$ collisions Ni  _et al._ (2010)
inside a cavity can be investigated. In the many-body regime of a gas of
ultracold molecules Zhang  _et al._ (2021); Li  _et al._ (2021); Konishi  _et
al._ (2021), coherent control of molecular orientation by the cavity field is
an enticing possibility Sidler  _et al._ (2023). Along these lines, recently,
universal behavior has been reported for pair-polariton states emerging when a
Fermi gas is coupled to a cavity Konishi  _et al._ (2021). The unparalleled
precision and controlability provided in ultracold physics, can aid to uncover
mechanisms behind the effects observed in polaritonic chemistry. Achieving
mechanistic understanding of these phenomena through ultracold molecules will
be an important leap forward.

###### Acknowledgements.

The authors thank Dan Stamper-Kurn and Giacomo Valtolina for stimulating
discussions. We acknowledge support from the NSF through a grant for ITAMP at
Harvard University.

## References

  * Kohler  _et al._ (1995) B. Kohler, J. L. Krause, F. Raksi, K. R. Wilson, V. V. Yakovlev, R. M. Whitnell, and Y. Yan, “Controlling the future of matter,” Accounts of Chemical Research 28, 133–140 (1995).
  * Eigler and Schweizer (1990) D. M. Eigler and E. K. Schweizer, “Positioning single atoms with a scanning tunnelling microscope,” Nature 344, 524–526 (1990).
  * Sugimoto  _et al._ (2005) Y. Sugimoto, M. Abe, S. Hirayama, N. Oyabu, Ó. Custance, and S. Morita, “Atom inlays performed at room temperature using atomic force microscopy,” Nature Materials 4, 156–159 (2005).
  * Gross  _et al._ (2009) L. Gross, F. Mohn, N. Moll, P. Liljeroth, and G. Meyer, “The chemical structure of a molecule resolved by atomic force microscopy,” Science 325, 1110–1114 (2009).
  * Tannor  _et al._ (1986) D. J. Tannor, R. Kosloff, and S. A. Rice, “Coherent pulse sequence induced control of selectivity of reactions: Exact quantum mechanical calculations,” The Journal of Chemical Physics 85, 5805–5820 (1986).
  * Shapiro and Brumer (2003) M. Shapiro and P. Brumer,  _Principles of the Quantum Control of Molecular Processes_ (Wiley-VCH, 2003).
  * Peirce  _et al._ (1988) A. P. Peirce, M. A. Dahleh, and H. Rabitz, “Optimal control of quantum-mechanical systems: Existence, numerical approximation, and applications,” Phys. Rev. A 37, 4950–4964 (1988).
  * Anderegg  _et al._ (2019) L. Anderegg, L. W. Cheuk, Y. Bao, S. Burchesky, W. Ketterle, K.-K. Ni, and J. M. Doyle, “An optical tweezer array of ultracold molecules,” Science 365, 1156–1158 (2019).
  * Kaufman and Ni (2012) A. M. Kaufman and K.-K. Ni, “Quantum science with optical tweezer arrays of ultracold atoms and molecules,” Nature Physics 17, 1324–1333 (2012).
  * Anderegg  _et al._ (2018) L. Anderegg, B. L. Augenbraun, Y. Bao, S. Burchesky, L. W. Cheuk, W. Ketterle, and J. M. Doyle, “Laser cooling of optically trapped molecules,” Nature Physics 14, 890–893 (2018).
  * Park  _et al._ (2023) J. J. Park, Y.-K. Lu, A. O. Jamison, and W. Ketterle, “Magnetic trapping of ultracold molecules at high density,” Nature Physics (2023), 10.1038/s41567-023-02141-0.
  * Ni  _et al._ (2010) K.-K. Ni, S. Ospelkaus, D. Wang, G. Quéméner, B. Neyenhuis, M. H. G. de Miranda, J. L. Bohn, J. Ye, and D. S. Jin, “Dipolar collisions of polar molecules in the quantum regime,” Nature 464, 1324–1328 (2010).
  * Safronova  _et al._ (2018) M. S. Safronova, D. Budker, D. DeMille, Derek F. Jackson Kimball, A. Derevianko, and Charles W. Clark, “Search for new physics with atoms and molecules,” Rev. Mod. Phys. 90, 025008 (2018).
  * DeMille (2002) D. DeMille, “Quantum computation with trapped polar molecules,” Phys. Rev. Lett. 88, 067901 (2002).
  * Blackmore  _et al._ (2018) J. A. Blackmore, L. Caldwell, P. D. Gregory, E. M. Bridge, R. Sawant, J. Aldegunde, J. Mur-Petit, D. Jaksch, J. M Hutson, B. E. Sauer, M. R. Tarbutt, and S. L. Cornish, “Ultracold molecules for quantum simulation: rotational coherences in caf and rbcs,” Quantum Science and Technology 4, 014010 (2018).
  * Carr  _et al._ (2009) L. D. Carr, D. DeMille, R. V. Krems, and J. Ye, “Cold and ultracold molecules: science, technology and applications,” New Journal of Physics 11, 055049 (2009).
  * Quéméner and Julienne (2012) G. Quéméner and P. S. Julienne, “Ultracold molecules under control!” Chemical Reviews 112, 4949–5011 (2012), pMID: 22921011.
  * Winkler  _et al._ (2007) K. Winkler, F. Lang, G. Thalhammer, P. v. d. Straten, R. Grimm, and J. Hecker Denschlag, “Coherent optical transfer of feshbach molecules to a lower vibrational state,” Phys. Rev. Lett. 98, 043201 (2007).
  * Vitanov  _et al._ (2017) N. V. Vitanov, A. A. Rangelov, B. W. Shore, and K. Bergmann, “Stimulated raman adiabatic passage in physics, chemistry, and beyond,” Rev. Mod. Phys. 89, 015006 (2017).
  * Hutchison  _et al._ (2012) J. A. Hutchison, T. Schwartz, C. Genet, E. Devaux, and T. W. Ebbesen, “Modifying chemical landscapes by coupling to vacuum fields,” Angewandte Chemie International Edition 51, 1592–1596 (2012).
  * Ebbesen (2016) T. W. Ebbesen, “Hybrid light–matter states in a molecular and material science perspective,” Accounts of Chemical Research 49, 2403–2412 (2016).
  * Hutchison  _et al._ (2013) J. A. Hutchison, A. Liscio, T. Schwartz, A. Canaguier-Durand, C. Genet, V. Palermo, P. Samorì, and T. W. Ebbesen, “Tuning the work-function via strong coupling,” Advanced Materials 25, 2481–2485 (2013).
  * Orgiu  _et al._ (2015) E. Orgiu, J. George, J. A. Hutchison, E. Devaux, J. F. Dayen, B. Doudin, F. Stellacci, C. Genet, J. Schachenmayer, C. Genes, G. Pupillo, P. Samorì, and T. W. Ebbesen, “Conductivity in organic semiconductors hybridized with the vacuum field,” Nat. Mater. 14, 1123–1129 (2015).
  * Galego  _et al._ (2016) J. Galego, F. J. Garcia-Vidal, and J. Feist, “Suppressing photochemical reactions with quantized light fields,” Nature Communications 7, 13841 (2016).
  * Flick  _et al._ (2017) J. Flick, M. Ruggenthaler, H. Appel, and A. Rubio, “Atoms and molecules in cavities, from weak to strong coupling in quantum-electrodynamics (qed) chemistry,” Proceedings of the National Academy of Sciences 114, 3026–3034 (2017).
  * Schäfer  _et al._ (2019) C. Schäfer, M. Ruggenthaler, H. Appel, and A. Rubio, “Modification of excitation and charge transfer in cavity quantum-electrodynamical chemistry,” Proceedings of the National Academy of Sciences 116, 4883–4892 (2019).
  * Bloch  _et al._ (2022) J. Bloch, A. Cavalleri, V. Galitski, M. Hafezi, and A. Rubio, “Strongly correlated electron–photon systems,” Nature 606, 41–48 (2022).
  * Schlawin  _et al._ (2022) F. Schlawin, D. M. Kennes, and M. A. Sentef, “Cavity quantum materials,” Applied Physics Reviews 9, 011312 (2022).
  * Paravicini-Bagliani  _et al._ (2019) G. L. Paravicini-Bagliani, F. Appugliese, E. Richter, S. Fallahi, F. Valmorra, J. Keller, M. Beck, N. Bartolo, C. Rössler, T. Ihn, K. Ensslin, C. Ciuti, G. Scalari, and J. Faist, “Magneto-transport controlled by landau polariton states,” Nature Physics 15, 186–190 (2019).
  * Appugliese  _et al._ (2022) F. Appugliese, J. Enkner, G. L. Paravicini-Bagliani, M. Beck, C. Reichl, W. Wegscheider, G. Scalari, C. Ciuti, and J. Faist, “Breakdown of topological protection by cavity vacuum fields in the integer quantum Hall effect,” Science 375, 1030–1034 (2022).
  * Rokaj  _et al._ (2022) V. Rokaj, M. Ruggenthaler, F. G. Eich, and A. Rubio, “Free electron gas in cavity quantum electrodynamics,” Phys. Rev. Res. 4, 013012 (2022).
  * Rokaj  _et al._ (2023a) V. Rokaj, J. Wang, J. Sous, M. Penz, M. Ruggenthaler, and A. Rubio, “Weakened topological protection of the quantum hall effect in a cavity,” Phys. Rev. Lett. 131, 196602 (2023a).
  * Ahn  _et al._ (2023) W. Ahn, J. F. Triana, F. Recabal, F. Herrera, and B. S. Simpkins, “Modification of ground-state chemical reactivity via light–matter coherence in infrared cavities,” Science 380, 1165–1168 (2023).
  * Ruggenthaler  _et al._ (0) M. Ruggenthaler, D. Sidler, and A. Rubio, “Understanding polaritonic chemistry from ab initio quantum electrodynamics,” Chemical Reviews 0, null (0), pMID: 37729114.
  * Campos-Gonzalez-Angulo  _et al._ (2023) J. A. Campos-Gonzalez-Angulo, Y. R. Poh, M. Du, and J. Yuen-Zhou, “Swinging between shine and shadow: Theoretical advances on thermally activated vibropolaritonic chemistry,” The Journal of Chemical Physics 158, 230901 (2023).
  * Mandal  _et al._ (2023) A. Mandal, M. A. D. Taylor, B. M. Weight, E. R. Koessler, X. Li, and P. Huo, “Theoretical advances in polariton chemistry and molecular cavity quantum electrodynamics,” Chemical Reviews 123, 9786–9879 (2023).
  * Garcia-Vidal  _et al._ (2021) F. J. Garcia-Vidal, C. Ciuti, and T. W. Ebbesen, “Manipulating matter by strong coupling to vacuum fields,” Science 373, eabd0336 (2021).
  * Rider and Barnes (2021) M. S. Rider and W. L. Barnes, “Something from nothing: linking molecules with virtual light,” Contemporary Physics 62, 217–232 (2021).
  * Rokaj  _et al._ (2018) V. Rokaj, D. M. Welakuh, M. Ruggenthaler, and A. Rubio, “Light–matter interaction in the long-wavelength limit: no ground-state without dipole self-energy,” Journal of Physics B: Atomic, Molecular and Optical Physics 51, 034005 (2018).
  * Cohen-Tannoudji  _et al._ (1997) C. Cohen-Tannoudji, J. Dupont-Roc, and G. Grynberg, “Photons and atoms-introduction to quantum electrodynamics,” Photons and Atoms-Introduction to Quantum Electrodynamics Wiley-VCH, February 1997. , 486 (1997).
  * Morse (1929) P. M. Morse, “Diatomic molecules according to the wave mechanics. ii. vibrational levels,” Phys. Rev. 34, 57–64 (1929).
  * Park  _et al._ (2001) S. J. Park, S. W. Suh, Y. S. Lee, and G.-H. Jeung, “Theoretical study of the electronic states of the rb2 molecule,” Journal of Molecular Spectroscopy 207, 129–135 (2001).
  * De Bernardis  _et al._ (2018) D. De Bernardis, P. Pilar, T. Jaako, S. De Liberato, and P. Rabl, “Breakdown of gauge invariance in ultrastrong-coupling cavity qed,” Physical Review A 98, 053819 (2018).
  * (44) See Supplemental Material for the details on the molecular model and Hamiltonian, the Jaynes-Cummings model, as well as the STIRAP photoassociation to different states.
  * Griffiths (1995) D. J. Griffiths,  _Introduction to Quantum Mechanics_ (Prentice Hall, 1995).
  * Shore and Knight (1993) Bruce W. Shore and Peter L. Knight, “The jaynes-cummings model,” Journal of Modern Optics 40, 1195–1238 (1993).
  * Scully and Zubairy (1997) M. O. Scully and M. S. Zubairy,  _Quantum Optics_ (Cambridge University Press, 1997).
  * Vogel and Welsch (2006) W. Vogel and D.-G. Welsch,  _Quantum Optics_ (Wiley‐VCH Verlag GmbH $\&$ Co. KGaA, 2006).
  * Kramida  _et al._ (2022) A. Kramida, Yu. Ralchenko, J. Reader, and and NIST ASD Team, NIST Atomic Spectra Database (ver. 5.10), [Online]. Available: https://physics.nist.gov/asd [2023, July 10]. National Institute of Standards and Technology, Gaithersburg, MD. (2022).
  * Rokaj  _et al._ (2023b) V. Rokaj, S. I. Mistakidis, and H. R. Sadeghpour, “Cavity induced collective behavior in the polaritonic ground state,” SciPost Phys. 14, 167 (2023b).
  * Zhang  _et al._ (2021) Z. Zhang, L. Chen, K.-X. Yao, and C. Chin, “Transition from an atomic to a molecular bose–einstein condensate,” Nature 592, 708–711 (2021).
  * Li  _et al._ (2021) J.-R. Li, W. G. Tobias, K. Matsuda, C. Miller, G. Valtolina, L. De Marco, R. R. W. Wang, L. Lassablière, G. Quéméner, J. L. Bohn, and J. Ye, “Tuning of dipolar interactions and evaporative cooling in a three-dimensional molecular quantum gas,” Nature Physics 17, 1144–1148 (2021).
  * Konishi  _et al._ (2021) H. Konishi, K. Roux, V. Helson, and J.-P. Brantut, “Universal pair polaritons in a strongly interacting fermi gas,” Nature 596, 509–513 (2021).
  * Sidler  _et al._ (2023) D. Sidler, T. Schnappinger, A. Obzhirov, M. Ruggenthaler, M. Kowalewski, and A. Rubio, “Unraveling a cavity induced molecular polarization mechanism from collective vibrational strong coupling,” (2023), arXiv:2306.06004 [quant-ph] .

SUPPLEMENTARY MATERIAL

##  I Molecular Model

Here, we provide details about the potential energy surfaces (PESs) used for
$\textrm{Rb}_{2}$, the respective vibrational wavefunctions involved in the
STIRAP photoassociation, as well as the bare Franck-Condon (FC) factors. For
the description of the $\rm{Rb}_{2}$ we employ the well-known Morse model
Morse (1929)

$\displaystyle\mathcal{H}_{M}=-\frac{\hbar^{2}}{2\mu}\frac{\partial^{2}}{\partial
R^{2}}+D\left[e^{-2w(R-R_{0})}-2e^{-w(R-R_{0})}\right],$ (S1)

where $\mu$ is the reduced mass of the dimer, $D$ is the depth of the
potential well, $w$ the width of the potential, and $R_{0}$ denotes the
equilibrium bond distance Morse (1929). To perform the photoassociation
process in the main text we used the following PESs of $\rm{Rb}_{2}$:
$\Sigma^{+}_{u}$, $\Pi_{u}$ and $\Sigma^{+}_{g}$. The PESs are modeled with
the Morse potentials and its corresponding parameters were obtained from the
experimental data of Ref. Park  _et al._ (2001). The explicit characteristics
of these PESs are given in Table I.

State |  $D[\textrm{cm}^{-1}]$ |  $\omega_{e}[\textrm{cm}^{-1}]$ |  $R_{0}\left[\textrm{\AA}\right]$ |  $T[\textrm{cm}^{-1}]$ | Dissociation limits
---|---|---|---|---|---
$(1)^{3}\Sigma^{+}_{u}$ | 240 | 14 | 6.13 | 3587 | 5s+5s
$(1)^{3}\Pi_{u}$ | 6828 | 59 | 4.20 | 9779 | 5s+5p
$X^{1}\Sigma^{+}_{g}$ | 3827 | 59 | 4.23 | 0 | 5s+5s
Table 1: Experimentally obtained parameters for the Morse potentials
describing the potential energy surfaces of the $\textrm{Rb}_{2}$ molecule as
reported in Ref Park  _et al._ (2001). Note that $\omega_{e}=w\sqrt{2D/\mu}$.

The vibrational states were obtained through numerical diagonalization.
Characteristic examples of vibrational wavefunctions to explore
photoassociation are shown in Fig. S1. Particularly, in the main text, the
states $\Phi^{\Pi_{u}}_{60},\Phi^{\Pi_{u}}_{61},\Phi^{\Sigma^{+}_{u}}_{33}$
and $\Phi^{\Sigma^{+}_{g}}_{0}$ presented in Figs. S1(b), (c) have been
employed.

Figure S1: Vibrational wavefunctions of all the states considered in the main
text and in supplement (see legends). Several of the vibrational wavefunctions
are shifted such that all profiles are well discernible.

Further, in Fig. S2 we show several of the bare FC factors for the in
$\Pi_{u}$ PES, among which in panel (a) are the states
$\Phi^{\Pi_{u}}_{60},\Phi^{\Pi_{u}}_{61}$ which we used for the STIRAP
photoassociation in the main text. As it can be seen from Fig. S2(a) the FC
factors of $\Phi^{\Pi_{u}}_{60}$ and $\Phi^{\Pi_{u}}_{61}$ (in the red box)
for both the upwards (from $\Phi^{\Sigma^{+}_{u}}_{33}$) and the downwards (to
$\Phi^{\Sigma^{+}_{g}}_{0}$) transitions have the same sign. This explains the
fact that the corresponding FC factors of the lower polariton for both
transitions are enhanced. This leads subsequently to the enhancement of
photoassociation of molecules through the lower polariton branch as it was
discussed in the main text. Moreover, in Fig. S2(a) we present the FC factors
of the states $\Phi^{\Pi_{u}}_{57}$ and $\Phi^{\Pi_{u}}_{58}$ (in the blue
box) which are used for the photoassociation in Sec. V. Finally, in Fig. S2(b)
we provide additionally some of the FC factors of higher lying vibrational
states for the upwards (from $\Phi^{\Sigma^{+}_{u}}_{33}$) and downwards (to
$\Phi^{\Sigma^{+}_{g}}_{1}$) transitions. In the black box we highlight the FC
factors of the states $\Phi^{\Pi_{u}}_{205}$ and $\Phi^{\Pi_{u}}_{206}$ which
are used in Sec. IV.

Figure S2: FC factors of different vibrational states in the $\Pi_{u}$ PES
(see legend). The transitions to a) $\Phi^{\Sigma^{+}_{g}}_{0}$ and b)
$\Phi^{\Sigma^{+}_{g}}_{1}$ are shown. In both panels the upwards transitions
from $\Phi^{\Sigma^{+}_{u}}_{33}$ are provided which are used in the
$\Lambda$-scheme for photoassociation.

##  II Molecule-Photon Hamiltonian

The coupled molecule-photon system is described by the Pauli-Fierz Hamiltonian
in the length gauge Rokaj  _et al._ (2018); Cohen-Tannoudji  _et al._ (1997)

$\displaystyle\mathcal{H}=\mathcal{H}_{M}-\frac{\hbar\omega}{2}\frac{\partial^{2}}{\partial
q^{2}}+\frac{\hbar\omega}{2}\left(q-\frac{\lambda
d_{R}}{\sqrt{\hbar\omega}}\right)^{2}.$ (S2)

The cavity field is described by the displacement coordinate $q$ and its
conjugate momentum $\partial/\partial q$ while the cavity frequency is
$\omega=\pi c/L$. Equivalently, the cavity mode can be described in terms of
annihilation and creation operators
$a=(q+\partial_{q})/\sqrt{2},\;\textrm{and}\;a^{\dagger}=(q-\partial_{q})/\sqrt{2}$,
which satisfy bosonic commutations relations $[a,a^{\dagger}]=1$. After
expanding the quadratic term $\left(q-\frac{\lambda
d_{R}}{\sqrt{\hbar\omega}}\right)^{2}$ and formulating the Hamiltonian in
terms of annihilation and creation operators we have

$\displaystyle\mathcal{H}=\mathcal{H}_{M}+\hbar\omega\left(a^{\dagger}a+\frac{1}{2}\right)-\sqrt{\frac{\hbar\omega}{2}}\lambda
d_{R}\left(a+a^{\dagger}\right)+\frac{\lambda^{2}d^{2}_{R}}{2}.$ (S3)

The cavity field strength $\lambda=1/\sqrt{\epsilon_{0}\mathcal{V}}$ depends
on the effective cavity volume $\mathcal{V}=AL$ which is a product of the
mirror distance $L$ and the area of the cavity mirrors $A$. The cavity
frequency $\omega=\pi c/L$ and the mirror distance $L$ are also related, and
consequently the field strength $\lambda$ depends on $\omega$ as well, i.e.
$\lambda=\sqrt{\omega/\epsilon_{0}\pi cA}$. Further, $c$ is the speed of light
and $\epsilon_{0}$ the vacuum permittivity. We use the linear dipole moment
$d_{R}=\mu_{0}R$ with strength $\mu_{0}$, and express the molecule-photon
$\mathcal{H}$ in terms of the light-matter coupling constant $g$ which is
independent of the cavity frequency,

$\displaystyle\mathcal{H}=\mathcal{H}_{M}+\hbar\omega\left(a^{\dagger}a+\frac{1}{2}\right)-\omega
g\sqrt{\frac{\hbar}{2}}R\left(a+a^{\dagger}\right)+\frac{\omega
g^{2}}{2}R^{2},\;\;\textrm{where}\;\;g=\frac{\mu_{0}\lambda}{\sqrt{\omega}}=\frac{\mu_{0}}{\sqrt{A\pi
c\epsilon_{0}}}.$ (S4)

To treat the interacting molecule-photon coupled system we expand the
Hamiltonian in the basis consisting of the molecular eigenstates of the
$\Pi_{u}$ PES $\\{|\Phi^{\Pi_{u}}_{n}\rangle\\}$ and the eigenstates
$\\{|i\rangle\\}$ of the bosonic cavity mode Griffiths (1995). In this way, we
obtain the matrix form $\mathcal{H}^{ji}_{mn}:=\langle
j|\langle\Phi^{\Pi_{u}}_{m}|\mathcal{H}|\Phi^{\Pi_{u}}_{n}\rangle|i\rangle$ of
the Hamiltonian $\mathcal{H}$ which reads

$\displaystyle\mathcal{H}^{ji}_{mn}=\left[\mathcal{E}_{n}+\hbar\omega\left(i+\frac{1}{2}\right)\right]\delta_{mn}\otimes\delta_{ji}-\omega
g\sqrt{\frac{\hbar}{2}}\langle\Phi^{\Pi_{u}}_{m}|R|\Phi^{\Pi_{u}}_{n}\rangle\otimes\left[\sqrt{i}\delta_{j,i-1}+\sqrt{i+1}\delta_{j,i+1}\right]+\frac{\omega
g^{2}}{2}\langle\Phi^{\Pi_{u}}_{m}|R^{2}|\Phi^{\Pi_{u}}_{n}\rangle\otimes\delta_{ji}.$

We note that $\mathcal{E}_{n}$ are the eigenenergies of the molecular
Hamiltonian $\mathcal{H}_{M}$ described by the Morse model.

##  III Additional Details of STIRAP to $\Phi^{\Sigma^{+}_{g}}_{0}$

Let us now provide some additional data and supporting computations, about the
photoassociation of $\textrm{Rb}_{2}$ molecules to the vibrational ground
state $\Phi^{\Sigma^{+}_{g}}_{0}$ presented in the main text. Recall that in
the main text we only discussed the STIRAP to $\Phi^{\Sigma^{+}_{g}}_{0}$
through the lower polariton state $\Psi_{-}$ which exhibits the enhanced FC
factors around the avoided crossing. For completeness, in Fig. S3 we provide
the molecular yield when STIRAP is performed through bare molecular states
$\Phi^{\Pi_{u}}_{60}$ [Figs. S3(a)] and $\Phi^{\Pi_{u}}_{61}$ [Figs. S3(a)].
Further, the STIRAP through the upper polariton state $\Psi_{+}$, exhibiting
suppressed FC factors, for $\Omega_{R}=1.8\times 10^{-3}$ and
$\Omega_{R}=1.8\times 10^{-2}$ is illustrated in Figs. S3(c) and (d). In all
cases, initially ($t=0$) the state $\Phi^{\Sigma^{+}_{u}}_{33}$ is fully
populated. However, as time-evolves the population of the state
$\Phi^{\Sigma^{+}_{u}}_{33}$ decreases mostly due to the decay term
$\gamma_{\Sigma_{u}}$ of the $\Sigma_{u}$ state, see also Eq. Polaritonic
ultracold reactions: cavity controlled molecular photoassociation in the main
text. Specifically, when the Gaussian pulses are applied [Fig. 1(c) in the
main text] population is transferred from the initial state
$\Phi^{\Sigma^{+}_{u}}_{33}$ to the final one $\Phi^{\Sigma^{+}_{g}}_{0}$. For
STIRAP through $\Phi^{\Pi_{u}}_{60}$ [Fig. S3 (a)] the amount of $12.4\%$ of
rubidium atoms turn into molecules, while for the $\Phi^{\Pi_{u}}_{61}$ [Fig.
S3(b)] we observe only $3.9\%$ production of $\textrm{Rb}_{2}$. For STIRAP
through the upper polariton $\Psi_{+}$ the molecular yield is extremely
suppressed to $\sim 0.01\%$ both for weak [Fig. S3(c)] and strong [Fig. S3(d)]
coupling, as expected from the substantially suppressed FC factors given in
the main text [Fig. 3].

Figure S3: Full time-evolution of the populations of the three states (see
legend) involved in the STIRAP. The STIRAP for a) $\Phi^{\Pi_{u}}_{60}$ and b)
$\Phi^{\Pi_{u}}_{61}$ is shown. The STIRAP for the upper polariton $\Psi_{+}$
for c) $g=0.001$ and d) $g=0.01$ is demonstrated. The time evolution of the
population during STIRAP to $\Phi^{\Sigma^{+}_{g}}_{0}$ through the upper
polariton is also depicted. The upper polariton features very low
photoassociation of $\textrm{Rb}_{2}$ molecules due to the very strongly
suppressed FC factors in both transitions.

###  III.1 Comparison to the Jaynes-Cummings model

In the main text we discussed that in the weak coupling regime,
$\Omega_{R}=1.8\times 10^{-3}$, the hybridization between the vibrational
states and the cavity photons, but also the corresponding FC factors of the
polariton states can be captured using the Jaynes-Cummings model Vogel and
Welsch (2006); Scully and Zubairy (1997). Indeed, for the case where the
interaction term between the two-level system and light is
$\mathcal{H}_{\rm{int}}=-(\hbar\Omega_{R}/2)\sigma^{+}\hat{a}+\textrm{h.c.}$
it can be proven that the JC states are Vogel and Welsch (2006)

$\displaystyle|+\rangle=\alpha_{+}|\Phi^{\Pi_{u}}_{n}\rangle|1\rangle-\alpha_{-}|\Phi^{\Pi_{u}}_{n+1}\rangle|0\rangle\;\textrm{and}\;|-\rangle=\alpha_{-}|\Phi^{\Pi_{u}}_{n}\rangle|1\rangle+\alpha_{+}|\Phi^{\Pi_{u}}_{n+1}\rangle|0\rangle\;\textrm{where}\;\alpha_{\pm}=\frac{1}{\sqrt{2}}\sqrt{1\pm\frac{\Delta}{\sqrt{\Delta^{2}+\Omega^{2}_{R}}}}.$
(S6)

Here, $\Delta$ is the detuning between the cavity mode and the two-level
transition. The minus sign in the interaction term of the JC model stems from
the form of the length gauge Hamiltonian and the fact that the transition
dipole matrix element is positive in our example
$\langle\Phi^{\Pi_{u}}_{60}|R|\Phi^{\Pi_{u}}_{61}\rangle>0$. Using the
analytic expression for the JC states in Fig. S4 we show the FC factors of the
JC states $|\pm\rangle$ for the two values of the light-matter coupling
strength and as a function of the detuning $\Delta$. The FC factors for the
upwards transition from $\Phi^{\Sigma^{+}_{u}}_{33}$ to the states in the
$\Pi_{u}$ PES are illustrated in Fig. S4(a), while the FC factors for the
downwards transition to the vibrational ground state
$\Phi^{\Sigma^{+}_{g}}_{0}$ are provided in Fig. S4(b). For both transitions
we observe that the FC factor of the lower polariton is enhanced around the
resonance point $\Delta=0$, while the FC factor of the upper polariton is
suppressed. This behavior is in agreement with the polaritonic FC factors
obtained from the full length gauge Hamiltonian [Eq. (1) in the main text]
which includes the counter-rotating terms and the dipole self-energy. For the
stronger light-matter coupling where $\Omega_{R}=1.8\times 10^{-2}$, however,
the JC model fails to capture the further modification of the FC factors at
the resonance point ($\Delta=0$). The JC states at $\Delta=0$ are independent
of the coupling strength and subsequently they become insensitive to a further
increase of the coupling. This is a consequence of the rotating-wave
approximation.

Figure S4: Bare and polaritonic FC factors of $\rm{Rb}_{2}$ for
$\Omega_{R}=1.8\times 10^{-3}$ and $\Omega_{R}=1.8\times 10^{-2}$ computed
within the JC model. a) FC factors for the downwards transition to
$\Phi^{\Sigma^{+}_{u}}_{0}$ from $\Phi^{\Pi_{u}}_{60},\Phi^{\Pi_{u}}_{61}$ and
the respective JC states $|\pm\rangle$ (see legend). b) FC factors regarding
the upward transition from $\Phi^{\Sigma^{+}_{u}}_{33}$ to
$\Phi^{\Pi_{u}}_{60},\Phi^{\Pi_{u}}_{61}$ and the JC states $|\pm\rangle$ (see
legend). In both panels at resonance $|-\rangle$ shows enhanced FC factors,
while the FC factors of $|+\rangle$ are suppressed. The JC model for weak
coupling $\Omega_{R}=1.8\times 10^{-3}$ is in agreement with the full
Hamiltonian [Eq. (1)]. However, the JC model for strong coupling
$\Omega_{R}=1.8\times 10^{-2}$ disagrees with the full model as no further
effect on the FC factors is observed at the resonance point.

##  IV STIRAP Photoassociation to $\Phi^{\Sigma^{+}_{g}}_{1}$

For completeness, next we provide a second example where polariton formation
is advantageous for the photoasociation of $\textrm{Rb}_{2}$ molecules. In
this case, we consider higher-lying vibrational states in the $\Pi_{u}$ PES,
namely $\Phi^{\Pi_{u}}_{205}$ and $\Phi^{\Pi_{u}}_{206}$, aiming to
photoassociate $\textrm{Rb}_{2}$ to the first excited vibrational state in
$\Sigma^{+}_{g}$, i.e. $\Phi^{\Sigma^{+}_{g}}_{1}$. The respective energy
spectrum for the vibrational states $\Phi^{\Pi_{u}}_{205}$ and
$\Phi^{\Pi_{u}}_{206}$ strongly coupled to the cavity with
$g=10^{-2}\textrm{a.u.}$ using the length gauge Hamiltonian is provided in
Fig. S5 taking into account five photon Fock states. A magnified version of
the lowest avoided crossing featuring a normalized Rabi splitting
$\Omega_{R}=4.14\times 10^{-2}$ which confirms the strong light-matter
interaction is illustrated in the inset of Fig. S5.

Figure S5: Polaritonic energy spectra for the vibrational states $n=205$ and
$m=206$ in the $\Pi_{u}$ PES hybridized with the vacuum cavity photons. The
system is considered under strong coupling, i.e. $g=10^{-2}$, and 5 Fock
states are required to converge the lowest avoided crossing shown in the inset
which features a normalized Rabi splitting $\Omega_{R}=4.14\times 10^{-2}$.

Subsequently, we compute the polaritonic FC factors for the upwards
transition, namely from $\Phi^{\Sigma^{+}_{u}}_{33}$ to the states in the
$\Pi_{u}$ PES, and the transition from the $\Pi_{u}$ PES to
$\Phi^{\Sigma^{+}_{g}}_{1}$ shown in Fig. S6. As it can be seen from Fig.
S6(a), throughout the whole frequency range the FC factor of the upper
polariton $\Psi_{+}$ is significantly enhanced, while the lower polariton
exhibits strong suppression, especially at resonance indicated by the vertical
dashed line. Furthermore, inspecting Fig. S6(b) it turns out that at the
resonance point (vertical dashed line) and before it, the FC factor of the
upper polariton is enhanced as compared to the bare state
$\Phi^{\Pi_{u}}_{205}$, while shortly after the resonance is getting
suppressed. On the other hand, the FC factor of the lower polariton $\Psi_{-}$
has always a lower value than the bare state $\Phi^{\Pi_{u}}_{205}$. Thus, it
can be concluded that for both transitions the upper polariton demonstrates an
enhanced FC factor at the resonance point (vertical dashed line) as compared
to the bare vibrational states $\Phi^{\Pi_{u}}_{205,206}$. This fact will be
crucial for the STIRAP photoassociation as we will argue below. It is
important to mention that in this example it is the upper polariton that
exhibits enhancement because the transition dipole matrix element
$\langle\Phi^{\Pi_{u}}_{205}|R|\Phi^{\Pi_{u}}_{206}\rangle<0$ is negative.

Figure S6: Bare and polaritonic FC factors for the downwards transition a)
from the states in the $\Pi_{u}$ PES to the $\Phi^{\Sigma^{+}_{g}}_{1}$ (see
legend), and b) for the upwards transition from the state
$\Phi^{\Sigma^{+}_{u}}_{33}$ to the states in the $\Pi_{u}$ PES (see legend).
For both transitions the upper polariton exhibits enhanced FC factors around
the resonance while the FC factor of the lower polariton is strongly
suppressed. The system is in the strong coupling regime with
$\Omega_{R}=4.14\times 10^{-2}$.

Having obtained the polaritonic FC factors of the vibrational states
$\Phi^{\Pi_{u}}_{205}$ and $\Phi^{\Pi_{u}}_{205}$ coupled to the cavity we
continue by performing STIRAP utilizing the $\Lambda$-scheme described in the
main text, see Eq. (Polaritonic ultracold reactions: cavity controlled
molecular photoassociation). In the $\Lambda$-scheme the state
$\Phi^{\Sigma^{+}_{u}}_{33}$ is initially fully populated and it is
photoassociated to $\Phi^{\Sigma^{+}_{g}}_{1}$. The molecular production using
the bare vibrational states $\Phi^{\Pi_{u}}_{205}$ and $\Phi^{\Pi_{u}}_{205}$
is shown in Figs. S7 (a) and (b) where the percentage of $12.5\%$ and $28.7\%$
of molecular photoassociation respectively is observed. The corresponding
populations via molecular photoassociation through the polaritonic states are
presented in Fig. S8. It is found that the upper polariton exhibits a
significantly enhanced molecular photoassociation which reaches $41.4\%$ [Fig.
S8(a)]. On the other hand, the photoassociation through the lower polariton
[Fig. S8(b)] is extremely suppressed and in particular much lower than $1\%$.
The significantly enhanced photoassociation of the upper polariton and the
respective suppression of the molecular yield in the lower polariton are
traced back to the behavior of their respective FC factors shown in Fig. S6.
More precisely, in Figs. S6(a) and (b) it was shown that at the resonance
point (vertical dashed line) the upper polariton state $\Psi_{+}$ acquired an
enhanced FC factor as compared to the bare vibrational states
$\Phi^{\Pi_{u}}_{205}$ and $\Phi^{\Pi_{u}}_{205}$.

Figure S7: Temporal evolution of the populations of the three states (see
legend) involved in the STIRAP. In a) the STIRAP through the state
$\Phi^{\Pi_{u}}_{205}$ is shown and we find $12.5\%$ molecular
photoassociation. In b) the STIRAP is performed via the state
$\Phi^{\Pi_{u}}_{206}$ and $28.7\%$ molecular photoassociation is identified.

Figure S8: Time-evolution for the STIRAP photoassociation performed through
the polariton states. In a) the STIRAP is performed through $\Psi_{+}$
resulting in $41.4\%$ molecular photoassociation, while in b) the STIRAP is
made through $\Psi_{-}$ leading to photoassociation below $1\%$.

##  V Example of No Polaritonic Advantage in Photoassociation

In the main text we argued that STIRAP photoassociation through the hybrid
polariton states has a clear advantage only when the FC factors of the same
polariton branch are enhanced for both transitions following the
$\Lambda$-scheme. Here, we provide a characteristic example where the
polariton states do not exhibit an advantage for the photoassociation of
molecules. We consider the molecular states $\Phi^{\Pi_{u}}_{57}$ and
$\Phi^{\Pi_{u}}_{58}$ which are strongly coupled to the cavity with light-
matter coupling strength $g=10^{-2}\textrm{a.u.}$. In this case, the lowest
avoided crossing of the system (not shown) exhibits normalized Rabi splitting
$\Omega_{R}=1.78\times 10^{-2}$. The resulting FC factors of the upper and
lower polariton states that participate in the first avoided crossing are
depicted in Fig. S9. It can be readily seen from Fig. S9(a) that after the
resonance point, indicated by the vertical dashed line, the FC factor of the
upper polariton is larger than the bare ones of the molecular states, while
the lower polariton has a much lower FC factor. Before the resonance point the
FC factors of both polariton states are lower than the vibrational state
$\Phi^{\Pi_{u}}_{57}$. Exactly at resonance, both polaritons have lower FC
factor than the molecular state $\Phi^{\Pi_{u}}_{57}$. In Fig. S9(b) the
situation is simpler and we observe that throughout the whole frequency window
where the avoided crossing occurs, the lower polariton $\Psi_{-}$ has a
substantially larger FC factor than both bare molecular states
$\Phi^{\Pi_{u}}_{57,58}$, while the upper polariton is significantly
suppressed. At resonance in Fig. S9(b) we see the FC factors of the polariton
states are maximally modified.

Figure S9: Bare and polaritonic FC factors for a) the upwards transition from
the state $\Phi^{\Sigma^{+}_{u}}_{33}$ to the states in the $\Pi_{u}$ PES (see
legend), and b) for the downwards transition from the states in the $\Pi_{u}$
PES (see legend) to the state $\Phi^{\Sigma^{+}_{g}}_{0}$. In a) the upper
polariton acquires an enhanced FC factor while in b) the situation is
inverted.

Having at hand the polaritonic FC factors of the vibrational states
$\Phi^{\Pi_{u}}_{57}$ and $\Phi^{\Pi_{u}}_{58}$ we continue by performing
STIRAP. The corresponding $\Lambda$-scheme is initiated with the state
$\Phi^{\Sigma^{+}_{u}}_{33}$ fully populated, and we photoassociate to
$\Phi^{\Sigma^{+}_{g}}_{0}$. The emergent molecular production using the bare
vibrational states $\Phi^{\Pi_{u}}_{57}$ and $\Phi^{\Pi_{u}}_{58}$ is shown in
Figs. S10 (a) and (b) which leads to $3.3\%$ and $\sim 0.5\%$ of molecular
photoassociation respectively. Importantly, by considering the STIRAP
photoassociation through the polaritonic states we observe $\sim 0.02\%$
molecular yield for the upper polariton [Fig. S10(c)] and $\sim 1.5\%$ for the
lower polariton [Fig. S10(d)]. Thus, neither of the polariton states shows an
advantage for the photoassociation of $\textrm{Rb}_{2}$ molecules. This is the
case because none of the polariton branches exhibits enhanced FC factors for
both transitions in the $\Lambda$-scheme.

Figure S10: Time-evolution of the populations of the three states (see
legends) involved in the STIRAP. In a) the STIRAP through
$\Phi^{\Pi_{u}}_{57}$ is shown demonstrating $3.3\%$ photoassociation. In b)
the STIRAP photoassociation through $\Phi^{\Pi_{u}}_{58}$ is $0.5\%$. In c)
and d) the STIRAP for the upper ($\Psi_{+}$) and the lower ($\Psi_{-}$)
polariton are shown for $\Omega_{R}=1.78\times 10^{-2}$ reaching $0.02\%$ and
$\sim 1.5\%$ photoassociation respectively. None of the polaritons shows
enhanced molecular photoassociation as compared to $\Phi^{\Pi_{u}}_{57}$.


# Analytic Estimation of Region of Attraction of an LQR Controller for Torque
Limited Simple Pendulum

Lukas Gross1,2, Lasse Maywald1, Shivesh Kumar1, Frank Kirchner1, and Christoph
Lüth2 This work was supported by the VeryHuman (Grant Number: 01IW20004)
project funded by the German Aerospace Center (DLR) with federal funds from
the Federal Ministry of Education and Research (BMBF).1 Robotics Innovation
Center, DFKI, 28359 Bremen, Germany.2 Cyber-Physical Systems, DFKI, 28359
Bremen, Germany.© 2022 IEEE. Personal use of this material is permitted.
Permission from IEEE must be obtained for all other uses, in any current or
future media, including reprinting/republishing this material for advertising
or promotional purposes, creating new collective works, for resale or
redistribution to servers or lists, or reuse of any copyrighted component of
this work in other works.

###### Abstract

Linear-quadratic regulators (LQR) are a well known and widely used tool in
control theory for both linear and nonlinear dynamics. For nonlinear problems,
an LQR-based controller is usually only locally viable, thus, raising the
problem of estimating the region of attraction (ROA). The need for good ROA
estimations becomes especially pressing for underactuated systems, as a
failure of controls might lead to unsafe and unrecoverable system states.
Known approaches based on optimization or sampling, while working well, might
be too slow in time critical applications and are hard to verify formally. In
this work, we propose a novel approach to estimate the ROA based on the
analytic solutions to linear ODEs for the torque limited simple pendulum. In
simulation and physical experiments, we compared our approach to a Lyapunov-
sampling baseline approach and found that our approach was faster to compute,
while yielding ROA estimations of similar phase space area.

##  I Introduction

It is well known that systems with linear dynamics can be stabilized optimally
w.r.t. a quadratic cost using a Linear-Quadratic Regulator (LQR) controller
[1], [2]. While they are straightforward to construct, LQR-controllers and the
variant TVLQR proved to be useful tackling problems involving nonlinear
dynamics. Possible applications range from stabilizing a fixed point using the
linear approximations LQR to stabilizing full trajectories by computing
multiple locally viable LQR-controllers along the trajectory [3, 4].

A sub-problem arising from LQR-methods is the estimation of the respective
Region of Attraction (ROA), i.e. the subset of the statespace in which the
system is guaranteed to evolve towards the desired final state [5]. This ROA
is not only constrained by the precision of the linear approximation that is
used, but also by underactuation due to unactuated degrees of freedom and
force or torque limitations.

Commonly, ROA estimation uses the Lyapunov criteria for stability. As known in
the study of dynamical [6] as well as power systems [7, 8, 9], for an
uncontrolled system the total energy is usually a very good candidate for a
Lyapunov function, as system losses ensure its time-derivative to be negative.
However, this does not translate to controlled systems, as a controller might
induce additional energy into the system. Because of this, the need to
construct suitable candidates for Lyapunov functions arises.

For an LQR-controller the naive ansatz is the cost-to-go, which is computed
during the LQR-controllers derivation, as it is at least locally valid around
the origin of the nonlinear dynamics. More refined methods include using sums-
of-squares (SOS) or convex optimization [10, 11, 12]. Any sublevel set of this
function, in which all Lyapunov conditions hold, when force/ torque
limitations are taken into account, is an inner estimation of the actual ROA.
The largest sublevel set can be determined via sampling [13, 14] or, again, by
using SOS optimization [15, 16].

The methods portrayed above, while working well in practice, have two main
weaknesses. Firstly, cost-to-go is quadratic and thus the resulting sublevel
sets are of ellipsoidal form. Lyapunov functions derived by SOS optimization
can yield different shapes, but in practice, it is often quadratic as well.
There is, however, no reason to expect the actual ROA to be ellipsoidal.
Hence, the resulting ROA potentially severly underestimates the actual region.
Secondly, ROA analysis is of special interest in safety-critical applications,
where, more often than not, formal verification is essential, but since
optimization is an intricate numerical approach and sampling depends on
statistics, they tend to be hard to grasp using formal methods [14].

(a)

(b)

Figure 1: (a) schematic of the pendulum, (b) long exposure shot of the
physical system

#### Contribution

As a first step towards an alternative which does not suffer from the above
shortcomings, we present an analytic approach to ROA estimations for the
stabilization of a torque limited simple pendulum in its unstable fixed point
(upright configuration). We chose a pendulum for its simple yet nonlinear
dynamics, which exhibit features of underactuation. We found that our approach
covered a ROA with area similar to that of the baseline, with a tendency to
outperform it for tight torque limitations. In addition, we found that
computational time required to compute the ROA with our proposed analytic
method is multiple orders of magnitude faster than the baseline method. The
approach has also been validated on a physical setup (see Fig. 1b) in
experiments.

#### Organization

In Section II, the closed loop dynamics of an LQR-controlled torque limited
pendulum are described in detail. In Section III, we use an analytic solution
to the initial value problem of the linearized system to derive sufficient
conditions for a system state to be part of the ROA. These findings are
evaluated in Section IV in simulation as well as on a physical setup and their
performance is compared to a Lyapunov-sampling approach as baseline. Section V
concludes the paper and highlights the future work.

##  II Closed Loop Dynamics of an LQR-controlled Torque Limited Simple
Pendulum

Consider a pendulum consisting of a rigid arm of length $l$ with an actuator
at its axis and a weight of mass $m$ attached to its end. For simplicity and
without loss of generality, consider the arm to be weightless. The equation of
motion for such a pendulum and its first order approximation around the
upright unstable fixed point $(\theta=\omega=0)$ are given as

$\displaystyle\ddot{\theta}$
$\displaystyle=\frac{1}{ml^{2}}\left(mgl\sin(\theta)-b\dot{\theta}+u\right)$
(1)
$\displaystyle\simeq\frac{1}{ml^{2}}\left(mgl\theta-b\dot{\theta}+u\right),$
(2)

where $g$ is the gravitational acceleration, $b$ a damping constant and $u$
the torque actuation. The direction of angle $\theta$, angular velocity
$\omega=\dot{\theta}$ and the torque are defined as in Fig. 1(a). Defining the
variables that way, the origin of the phase space corresponds to the upright
fixed point. With $\mathbf{x}=(\theta,\dot{\theta})^{T}$ the linear
approximation can be decomposed into a set of first order ordinary
differential equations (ODE) as

$\dot{\mathbf{x}}=\underbrace{\left(\begin{array}[]{cc}0&1\\\
\frac{g}{l}&-\frac{b}{ml^{2}}\\\
\end{array}\right)}_{\mathbf{A}}\mathbf{x}+\underbrace{\left(\begin{array}[]{c}0\\\
\frac{1}{ml^{2}}\\\ \end{array}\right)}_{\mathbf{B}}u.$ (3)

For a system of the form $\dot{\mathbf{x}}=\mathbf{A}\mathbf{x}+\mathbf{B}u$
and a quadratic cost function of the form
$\int_{0}^{\infty}\left(\mathbf{x}^{T}\mathbf{Q}\mathbf{x}+uRu\right)dt$, the
LQR-method yields an optimal controller $u(\mathbf{x})=-\mathbf{K}\mathbf{x}$,
which stabilizes the systems origin $\mathbf{x}^{*}=(0,0)^{T}$, where
$\mathbf{K}=\mathbf{R}^{-1}\mathbf{B}^{T}\mathbf{S}$ and $\mathbf{S}$ is the
solution of the algebraic Riccati equation
$0=\mathbf{SA}+\mathbf{A}^{T}\mathbf{S}-\mathbf{SBR}^{-1}\mathbf{B}^{T}\mathbf{S}+\mathbf{Q}.$
In this case, $\mathbf{K}$ is of the form
$\mathbf{K}=\left(K_{0},K_{1}\right)^{T}$. Inserting the LQR-controller into
the equations (1) and (2) yields the closed-loop dynamics and its
linearization:

$\displaystyle\ddot{\theta}$
$\displaystyle=-\frac{1}{ml^{2}}\left(mgl\sin(\theta)-K_{0}\theta-(b+K_{1})\dot{\theta}\right)$
(4) $\displaystyle\simeq-\frac{1}{ml^{2}}\left((mgl-
K_{0})\theta-(b+K_{1})\dot{\theta}\right).$ (5)

Furthermore, consider the actuators torque output to be limited, i.e.
$|u|\leq\bar{u}$. This results in the following equations of motion:

$\ddot{\theta}=\begin{cases}\frac{1}{ml^{2}}\left(mgl\sin(\theta)-b\dot{\theta}-\bar{u}\right)&u(\mathbf{x})<-\bar{u}\\\
\frac{1}{ml^{2}}\left(mgl\sin(\theta)-b\dot{\theta}+\bar{u}\right)&u(\mathbf{x})>\bar{u}\\\
\frac{1}{ml^{2}}\left(mgl\sin(\theta)-K_{0}\theta-(b+K_{1})\dot{\theta}\right)&\text{else}\end{cases}$
(6)

While we know that the origin is globally stable for (5), for the nonlinear
torque limited dynamics (6) the stability holds only locally. Thus, we need to
estimate the ROA $\mathcal{S}$, in which all states evolve to the origin as
time goes to infinity:

$\mathcal{S}=\left\\{\mathbf{\gamma}\,\middle|\,\mathbf{x}_{\mathbf{\gamma}}(t\to\infty)\to
0\right\\}.$ (7)

##  III Analytically Derived Region of Attraction Estimation

In this section we introduce our approach to estimate the ROA analytically. We
start by motivating an inner estimation of $\mathcal{S}$, which depends on the
time evolution $u_{\mathbf{\gamma}}(t)$ with initial value $\mathbf{\gamma}$.
We then approximate $u_{\mathbf{\gamma}}(t)$ by solving the linear dynamics
(5) and, to compensate for the approximation, add a heuristic constraint.

###  III-A Effects of Torque Limitations on the Region of Attraction

An LQR-controller is oblivious to torque limitations. Therefore we expect it
to perform well within the limits, but in general it is unclear how the system
will behave once the torques are clipped. Especially for safety-critical
applications this would be unacceptable. Thus we want to take a conservative
route and approximate the ROA as

$\tilde{\mathcal{S}}=\left\\{\mathbf{\gamma}\,\middle|\,\forall
t>0:|u_{\mathbf{\gamma}}(t)|\leq\bar{u}\right\\}\cap\mathcal{S}_{\text{unlim}},$
(8)

where $u_{\mathbf{\gamma}}(t)$ is the time evolution of the torque for a
system that with initial value $\mathbf{\gamma}$ and
$\mathcal{S}_{\text{unlim}}$ is the ROA of the system without torque
limitations (4). Simulations of (4) show that for the simple pendulum
$\mathcal{S}_{\text{unlim}}$ seems to be the entire phase space. Because of
this we can assume $\tilde{\mathcal{S}}\subset\mathcal{S}_{\text{unlim}}$ and
omit $\mathcal{S}_{\text{unlim}}$ from the following analysis. In the
following we will approximate $u_{\mathbf{\gamma}}(t)$ by explicitly solving
the initial value problem for (5) and inserting the solution into
$u(\mathbf{x})$.

###  III-B Analytic Solution to the Linearized Dynamics

Solutions to a linear ODE are linear combinations of exponential functions.
Their exponents depend on the solution of the characteristic polynomial. In
the case of (5) that is

$\kappa^{2}+\left(\frac{K_{1}+b}{ml^{2}}\right)\kappa+\left(\frac{K_{0}}{ml^{2}}-\frac{g}{l}\right)=0,$
(9)

with solutions

$\kappa_{0,1}=\frac{1}{2}\left(-\frac{K_{1}+b}{ml^{2}}\pm\sqrt{D}\right),$
(10)

where

$D=\left(\frac{K_{1}+b}{ml^{2}}\right)^{2}-4\left(\frac{K_{0}}{ml^{2}}-\frac{g}{l}\right).$
(11)

Assuming $D>0$ and $\kappa_{0,1}<0$ 111Due to the LQR-Controller being stable
and optimal for the linear dynamics, we expect exponentially decreasing
functions and thus negative real-valued exponents., the time evolution of
$\mathbf{x}$ is

$\mathbf{x}(t)=\left(\begin{array}[]{rr}1&1\\\ \kappa_{0}&\kappa_{1}\\\
\end{array}\right)\left(\begin{array}[]{rr}C_{0}e^{\kappa_{0}t}\\\
C_{1}e^{\kappa_{1}t}\\\ \end{array}\right),$ (12)

where

$\left(\begin{array}[]{rr}C_{0}\\\ C_{1}\\\
\end{array}\right)=\frac{1}{\sqrt{D}}\left(\begin{array}[]{rr}-\kappa_{1}&1\\\
\kappa_{0}&-1\\\ \end{array}\right)\mathbf{x}(0).$ (13)

Inserting $\mathbf{x}(t)$ into $u(\mathbf{x})$ yields

$u_{\text{lin},\mathbf{x}(0)}(t)=-\left(K_{0}+K_{1}\kappa_{0}\right)C_{0}e^{\kappa_{0}t}-\left(K_{0}+K_{1}\kappa_{1}\right)C_{1}e^{\kappa_{1}t}.$
(14)

Let $A_{i}=-\left(K_{0}+K_{1}\kappa_{i}\right)\kappa_{i}C_{i}$.
Differentiating $u_{\text{lin},\mathbf{x}(0)}(t)=0$ and solving
$\dot{u}_{\text{lin},\mathbf{x}(0)}(t^{*})=0$ for $t^{*}$, shows that $u(t)$
will have an extremum at

$t^{*}=-\frac{\ln\left(-\frac{A_{0}}{A_{1}}\right)}{\sqrt{D}},$ (15)

if $-A_{0}/A_{1}\geq 0$.

Therefore we can approximate $\tilde{\mathcal{S}}$ as

$\displaystyle\begin{split}\tilde{\mathcal{S}}\simeq\quad&\left\\{\mathbf{\gamma}\,\middle|\,|u_{\text{lin},\mathbf{\gamma}}(0)|\leq\bar{u}\right\\}\\\
\cap&\left\\{\mathbf{\gamma}\,\middle|\,t^{*}_{\mathbf{\gamma}}>0\implies|u_{\text{lin},\mathbf{\gamma}}(t^{*}_{\mathbf{\gamma}})|\leq\bar{u}\right\\}.\end{split}$
(16)

###  III-C Viability of the Linear Approximation

The linear approximation in the derivations above is only locally viable
around the origin, i.e. the pendulums upright position. As a heuristic to do
this we consider the difference of the gravity induced torques of the actual
and the linearized ODE and compare it to the torque limit. If that difference
is smaller than the torque limit, it is reasonable to assume that the
controller is still capable to compensate for the error:
$mgl|\sin(\theta)-\theta|\leq\bar{u}$. This heuristic only takes the angle as
an argument, as the approximation error in this case does not increase with
increasing angular velocity.

To conclude, our estimation of the ROA is characterized by three conditions on
a given phase space coordinate:

  * •

It must obey the heuristic $mgl|\sin(\theta)-\theta|\leq\bar{u}$

  * •

The initial torque must be within the limits
$|u_{\text{lin},\mathbf{\gamma}}(0)|<\bar{u}$.

  * •

If the applied torques will exhibit another extremum in the future $t^{*}>0$,
check whether this extremum violates the torque limitations
$|u_{\text{lin},\mathbf{\gamma}}(t^{*})|<\bar{u}$.

The resulting region in the phase space is given as

$\displaystyle\begin{split}\mathcal{S}_{\text{analytic}}=\quad&\left\\{\mathbf{\gamma}\,\middle|\,mgl|\sin(\theta)-\theta|\leq\bar{u}\right\\}\\\
\cap&\left\\{\mathbf{\gamma}\,\middle|\,|u_{\text{lin},\mathbf{\gamma}}(0)|\leq\bar{u}\right\\}\\\
\cap&\left\\{\mathbf{\gamma}\,\middle|\,t^{*}_{\mathbf{\gamma}}>0\implies|u_{\text{lin},\mathbf{\gamma}}(t^{*}_{\mathbf{\gamma}})|\leq\bar{u}\right\\}.\end{split}$
(17)

##  IV Results and Discussion

In the following section, we describe and discuss the simulations and
experiments in detail. We start by giving an overview of the pendulum’s
parameters, then we explain the simulations and physical experiments as well
as the software implementation222The open source software implementation can
be found at: https://github.com/dfki-ric-underactuated-
lab/torque_limited_simple_pendulum/tree/master/software/python/simple_pendulum/controllers/lqr/analytic_roa_estimation.
Afterwards the results will be presented and interpreted.

###  IV-A Pendulum Parameters

The physical system we used is described in detail in [17], a long exposure
shot of it is shown in Fig. 1(b). It has a total mass of $m\approx 0.676$ kg
and the distance from the axis of rotation to the center of mass is $l\approx
0.45$ m. Given these parameters the maximum torque that gravity can exert on
the system, i.e. at a $90^{\circ}$ angle, is $\tilde{u}:=mlg\approx 2.98$ Nm.
The actual motors maximum torque is $6$ Nm. In both physical experiments, we
simulated the underactuation by limiting the controllers torque outputs to
$\tilde{u}/2\approx 1.49$ Nm.

For the simulations, we chose three sets of parameters. The ‘normal’ set uses
the parameters of the real system. For the other two sets, we varied the
length and mass of the pendulum, while keeping $\tilde{u}$ fixed. That way we
kept the actuator torques in a reasonable range, but varied the inertial
properties of the system. We considered a ‘long’ system with $l/m=0.1$ and a
‘short’ one with $l/m=10$. The actual parameters and the resulting moments of
inertia are given in Tab. I. We also considered three different values of
torque limitation $\bar{u}\in\left\\{\tilde{u}/2,\ \tilde{u}/4,\
\tilde{u}/8\right\\}$. For both the simulations and physical experiments, the
LQRs cost parameters we used are $Q=\text{diag}(1,1)$ and $R=1$, additionally,
we assumed the damping to be around $b\approx 0.1\ \text{Nm s/rad}$. We did
not take Coulomb friction into account. We also considered three different
values of torque limitation $\bar{u}\in\left\\{\tilde{u}/2,\ \tilde{u}/4,\
\tilde{u}/8\right\\}$.

TABLE I: Paremeters for Simulation | mass [kg] | length [m] | inertia [kg m2]
---|---|---|---
normal | $0.676$ | $0.45$ | $0.137$
long | $0.174$ | $1.744$ | $0.531$
short | $1.744$ | $0.174$ | $0.0531$

###  IV-B Setup Description

To evaluate the validity and performance of the analytic ROA estimation
$\mathcal{S}_{\text{analytic}}$ as in (17), we compared it to a Lyapunov-
sampling ROA estimation $\mathcal{S}_{\text{Lyapunov}}$ as a baseline and also
put both in the context of the actual ROA $\mathcal{S}$ (see (7)) and the
region of all points that are stabilized without reaching the torque limits
$\tilde{\mathcal{S}}$ (see (8)).

Since $\mathcal{S}$ and $\tilde{\mathcal{S}}$ cant be computed directly, to
illustrate them, we ran $100000$ simulations for each parameter set and with
initial values randomly chosen in the range of $-\pi<\theta<\pi$ and $-10\
\text{rad/s}<\omega<10\ \text{rad/s}$. To integrate (6) numerically we used a
standard Runge-Kutta algorithm with step width $\Delta t=.1$ s and final time
$t_{\text{final}}=10$ s. We considered the initial values to be in
$\mathcal{S}$ if the integration converged to the origin within 5 digit
precision. If additionally along the trajectory the torque limits were never
exceeded, we considered the initial values to also be in
$\tilde{\mathcal{S}}$.

To run similar tests on the physical system, we used a PD-controller to
prepare the pendulum at a randomly chosen angle, then applied a torque
randomly chosen between $-5$ Nm and $5$ Nm, for $0.5$ s. Afterwards, we
switched to the LQR-controller and let the system evolve for $10$ s. We
considered the systems state at the time of switching as well as $50$ evenly
spaced (timewise) states along the trajectory as initial values
$\mathbf{\gamma}$. Again we recorded if the controller was able to stabilize
the pendulum upright ($\mathcal{S}$) and whether or not it met the torque
limits during the run ($\tilde{\mathcal{S}}$). In total we conducted $570$
runs resulting $28500$ data points.

In order to analyze the analytical approach, we implemented (17) as an oracle
style function $\mathbf{x}\mapsto\texttt{true},\text{iff}\
\mathbf{x}\in\mathcal{S}_{\text{analytic}}$ by numerically deriving
$\mathbf{K}$ and using it to compute $u_{\text{lin},\mathbf{\gamma}}$ and
$t^{*}$. We used this oracle function to estimate for each initial value from
simulation and experiment whether or not it lies inside
$\mathcal{S}_{\text{analytic}}$. In addition, to investigate the heuristic
used in the previous section, we also implemented an oracle function for
$\mathcal{S}_{\text{unbound}}$ as in (16) in an analogous way to the above.

For $\mathcal{S}_{\text{Lyapunov}}$, we took
$\mathbf{x}^{T}\mathbf{S}\mathbf{x}$ as a Lyapunov function333This is the
cost-to-go function of the system and known to be a viable Lyapunov function.
and computed the boundary $\rho$ of the largest sublevel set via sampling
states and checking the Lyapunov conditions [13]. Again we used this to
implement an oracle function: $\mathbf{x}\mapsto\texttt{true},\text{iff}\
\mathbf{x}^{T}\mathbf{S}\mathbf{x}\leq\rho$.

To demonstrate the usage of our ROA estimation, we set up a swing-up
experiment. We used a combination of an energy-shaping controller to bring the
system close to the phase space origin and, once the system is inside the
estimated ROA, switch to the LQR to stabilize the pendulum:

$u(\mathbf{x})=\begin{cases}-\mathbf{K}\mathbf{x},&\text{if}\
\mathbf{x}\in\mathcal{S}_{\text{analytic}}\\\ -c\dot{\theta}\Delta
E(\mathbf{x})+b\dot{\theta},&\text{else}\end{cases},$ (18)

with a positive constant $c$ and the difference in total energy w.r.t. origin
$\Delta E(\mathbf{x})=mgl(\cos(\theta)-1)-\frac{1}{2}ml^{2}\dot{\theta}^{2}$.

###  IV-C Validation of Heuristic

To evaluate the validity of the heuristic added in (17), we compared the
simulation results to the analytic ROA estimation with and without taking the
heuristic into account. Fig. 2 illustrates the results. Here we plotted all
four regions for the ‘long’ configuration with $\bar{u}=\tilde{u}/2$ as an
example. The grey and black regions depict $\mathcal{S}$ and
$\tilde{\mathcal{S}}$ respectively, the blue region is the analytic result
$\mathcal{S}_{\text{analytic}}$ and in red the region
$\mathcal{S}_{\text{unbound}}$ without the heuristic bound is drawn.

Figure 2: Effect of the heuristic on the ROA. Depicted are the analytic ROA
estimation both with ($\mathcal{S}_{\text{analytic}}$) and without heuristic
bounds ($\mathcal{S}_{\text{unbound}}$). As reference also $\mathcal{S}$ and
$\tilde{\mathcal{S}}$ from the simulation are shown.

Comparing $\mathcal{S}_{\text{unbound}}$ with $\mathcal{S}$ in the upper left
and lower right parts of the plot, we find that a significant number of
initial values are within $\mathcal{S}_{\text{unbound}}$ but have failed to
stabilize in simulation that means $\mathcal{S}_{\text{unbound}}$ is not a
conservative estimate of the ROA. Therefore, a heuristic to bound our
estimation is indeed necessary.

To ensure the applied heuristic is sufficient, we searched in all initial
values evaluated by simulation for false positives of the kind
$\mathbf{\gamma}\in\mathcal{S}_{\text{analytic}}\setminus\mathcal{S}$ and
found none. Thus, it seems to be a sufficient choice, but Fig. 2 also shows
that the heuristic bound reduces the area of the estimation by a lot.

###  IV-D Simulation Results

As an overview of the results from simulation, in Fig. 3, we plotted
$\mathcal{S}_{\text{analytic}}$ as blue, $\mathcal{S}_{\text{Lyapunov}}$ as
yellow and, again, $\mathcal{S}$ and $\tilde{\mathcal{S}}$ as black and grey
areas respectively. The subplots show the results for each set of parameters
and each torque limit considered. From looking at Fig. 3, we find that neither
the analytically estimated ROA is a superset of the opimization based ROA nor
the other way around. So none of the approaches seem to be clearly superior
regarding the ROA’s area. Due to their distinct difference in shape, they do
not overlap too much, therefore it might in fact be practical to combine both
methods for a larger overall ROA.

Figure 3: Comparison of analytic ROA ($\mathcal{S}_{\text{analytic}}$) to the
Lyapunov baseline ($\mathcal{S}_{\text{Lyapunov}}$). As reference also
$\mathcal{S}$ and $\tilde{\mathcal{S}}$ from the simulation are shown. The
systems mass and length vary between subfigures according to Tab. I and the
torque limit $\bar{u}$ varies as noted in the subcaptions.

For a quantitative comparison of the areas, we estimated the relative area of
$\mathcal{S}_{\text{analytic}}$ from the numbers of initial values inside each
ROA as $A_{rel}\simeq N_{analytic}/N_{Lyapunov}$ for each parameter set and
torque limit. The results are given in Tab. II. They show that the analytic
ROA’s area is bigger than the Lyapunov based one in four cases, namely the
‘normal’ and the ‘long’, where $\bar{u}\leq\tilde{u}/4$. For the ‘short’
system the Lyapunov-sampling estimation is strictly better. The numbers
suggest that the analytic approach gets better as the torque limits get more
severe. In our experiments we found only one exception to that rule. From
‘normal’ $\tilde{u}/4$ to $\tilde{u}/8$ the relative area decreased, but this
might be explained with statistical errors increasing with the smaller sample
sizes in smaller regions.

TABLE II: Relative Volume of the Analytic ROA Estimation | $\tilde{u}/2$ | $\tilde{u}/4$ | $\tilde{u}/8$
---|---|---|---
normal | $0.99$ | $1.183$ | $1.150$
long | $0.723$ | $1.152$ | $1.721$
short | $0.718$ | $0.72$ | $0.767$

As an additional aspect, we measured computation times for the oracle
functions of both $\mathcal{S}_{\text{analytic}}$ and
$\mathcal{S}_{\text{Lyapunov}}$, i.e. the time it took to prepare the function
for a new set of parameters, not the evaluation for a given initial value. The
mean of the results are:

$\displaystyle T_{\text{analytic}}=0.00136\text{ s}$ $\displaystyle
T_{\text{Lyapunov}}=12.8\text{ s}$

To put these numbers in perspective: computing $\mathcal{S}$ and
$\tilde{\mathcal{S}}$ took about $5$ hours for a single set of parameters. And
to check the stability of states that were not numerically computed, the use
of convex hull approximation of the point cloud would be required.

It took about $10^{5}$ times as long to compute the Lyapunov based prediction
as it did for the analytic method. This could be very relevant for robotic
applications where inertial parameters vary, such as grasping or carrying
external objects. For $\mathcal{S}_{\text{Lyapunov}}$ we sampled $500$ times.
Thus, even reducing this number to where the sampling stops converging
properly, a difference of $3$ to $4$ orders of magnitude between
$T_{\text{analytic}}$ and $T_{\text{Lyapunov}}$ would remain. However, one has
to keep in mind that the analytic method only yields an oracle style
prediction function, whereas the Lyapunov method yields a conservative
geometrically meaningful description of the estimated ROA in addition, which
is useful in many cases.

###  IV-E Experimental Results

The experimental results draw a similar picture to the simulations. As can be
seen in Fig. 4, the experimentally derived regions look very similar to the
ones from Fig. 3 (a). Most notably, all states inside an estimated ROA, both
analytic and Lyapunov-sampling, were stabilized by the controller. This
suggests that the systems dynamics were modeled well, and thus, it is
reasonable to use our ROA estimation on a real system.

Figure 4: Experimental evaluation of estimated ROA. Depicted are the analytic
ROA estimation ($\mathcal{S}_{\text{analytic}}$), the Lyapunov based method
($\mathcal{S}_{\text{Lyapunov}}$), and the experimentally derived
$\mathcal{S}$ and $\tilde{\mathcal{S}}$. The ‘normal’ parameters with
$\bar{u}=\tilde{u}/2$ were used.

In Fig. 5 the angle, velocity, and torque during the swing-up experiment are
shown. We see that the controller (18) managed to move the pendulum from the
lower fixed point to the upright position. In fact close inspection shows that
the angle did not quite reach the origin. This, as well as the non vanishing
torques is probably due to friction canceling out the controller torque.

The system reached the estimated ROA at $t_{\textbf{switch}}\approx 1.45$ s,
$\mathbf{x}\approx(0.416\pi\text{ rad},-4.4\text{ rad/s})$. A look at Fig. 3
(a) shows that the controller would not have switched at that time, if the
Lyapunov based approach would have been used instead. Thus, at least for this
single run, our ROA estimation provided a better prediction.

Figure 5: Swing-up of pendulum with imposed torque limit of
$\bar{u}=\tilde{u}/2$. A combination of an energy-shaping for swing-up and an
LQR-controller for stabilizing the upright position was used. The swing-up
controller is active in the brown and the LQR in the gray area. The switching
behavior was implemented using the analytically estimated ROA.

##  V Conclusion and Outlook

To summarize, we introduced an analytic approach to ROA estimation for an LQR-
controlled simple pendulum with torque limitations. We implemented this
approach and compared it to a baseline Lyapunov-sampling method in both
simulation and physical experiments. In the first case we varied the inertial
parameters of the pendulum as well as the enforced torque limitation. In
addition, we successfully demonstrated a swing-up controller that switches
between energy-shaping and LQR based on our approach as an exemplary use-case.

As far as our simulation results go, the new analytically estimated ROA is
between $0.7$ and $1.7$ times the size of the Lyapunov baseline. So generally
the performance is similar. Our results suggest that the analytic approach
becomes better with more severe torque limitations. Maybe more importantly the
two estimations overlap is comparatively small, so it might be of interest to
experiment with combining both methods. In addition to that, given a new set
of parameters, it was about $10^{5}$ times faster to compute an oracle
prediction function analytically than with the Lyapunov approach. The
experimental ROAs showed similar results. Most notably, even on the real
system, all states estimated to be inside the ROA were in fact stabilized by
the LQR-controller.

One rather obvious limitation of our approach is that all calculations done in
Section III are to some degree specific to the pendulums dynamics and higher
dimensional systems might be more intricate to solve. However, linear ODEs are
quite well studied and solutions to higher orders are known as well [18],
[19], so additional research might lead to a more general version of this
analytic method. Open problems regarding more complex systems include finding
a heuristic where the linear approximation breaks, as well as estimating the
ROA of the system without torque limitations $\mathcal{S}_{\text{unlim}}$.

Another possible extension would be to find solutions to the linearization of
the composed ODE from (6). In principle one could solve the components
individually and then try to find the states and times of the transitions from
the initial value problem. This way one could include initial values into the
ROA estimation, which start out in the area of saturated torque, as long as
they end up in the controllable area.

As mentioned in the introduction, one motivation for this work was to enable
the application of formal verification methods to the problem of ROA
estimation. A possible first step towards that goal would be to translate this
analytic approach, or parts of it, into a rigorous logic framework. A
promising candidate here would be differential dynamic logic [20] and the
associated prover keymaeraX [21]. Here, the use of a heuristic to bound the
linear approximations viability might become problematic. Further
investigations towards a better understanding of this is probably necessary.

## References

  * [1] B. D. O. Anderson and J. B. Moore, _Optimal control: linear quadratic methods_. One Lake Street Upper Saddle River, NJ, United States: Prentice-Hall, Inc., 1990.
  * [2] R. Tedrake, _Underactuated Robotics_ , 2022, ch. 8. [Online]. Available: http://underactuated.mit.edu
  * [3] W. Li and E. Todorov, “Iterative linear quadratic regulator design for nonlinear biological movement systems.” vol. 1, 01 2004, pp. 222–229.
  * [4] R. Tedrake, I. R. Manchester, M. Tobenkin, and J. W. Roberts, “Lqr-trees: Feedback motion planning via sums-of-squares verification,” _The International Journal of Robotics Research_ , vol. 29, no. 8, pp. 1038–1052, 2010\.
  * [5] H. K. Khalil, _Nonlinear systems; 3rd ed._ Upper Saddle River, NJ: Prentice-Hall, 2002. [Online]. Available: https://cds.cern.ch/record/1173048
  * [6] H.-D. Chiang and L. Fekih-Ahmed, “Quasi-stability regions of nonlinear dynamical systems: optimal estimations,” _IEEE Transactions on Circuits and Systems I: Fundamental Theory and Applications_ , vol. 43, no. 8, pp. 636–643, 1996.
  * [7] G. Gless, “Direct method of liapunov applied to transient power system stability,” _IEEE Transactions on Power Apparatus and Systems_ , no. 2, pp. 159–168, 1966.
  * [8] H. Kwatny, L. Bahar, and A. Pasrija, “Energy-like lyapunov functions for power system stability analysis,” _IEEE Transactions on Circuits and Systems_ , vol. 32, no. 11, pp. 1140–1149, 1985.
  * [9] M. A. Pai, _Energy function analysis for power system stability_. Springer Science & Business Media, 2012.
  * [10] T. A. Johansen, “Computation of lyapunov functions for smooth nonlinear systems using convex optimization,” _Automatica_ , vol. 36, no. 11, pp. 1617–1626, 2000. [Online]. Available: https://www.sciencedirect.com/science/article/pii/S0005109800000881
  * [11] G. Chesi, A. Garulli, A. Tesi, and A. Vicino, “Lmi-based computation of optimal quadratic lyapunov functions for odd polynomial systems,” _International Journal of Robust and Nonlinear Control_ , vol. 15, no. 1, pp. 35–49, 2005. [Online]. Available: https://onlinelibrary.wiley.com/doi/abs/10.1002/rnc.967
  * [12] B. Tibken and K. Dilaver, “Computation of subsets of the domain of attraction for polynomial systems,” in _Proceedings of the 41st IEEE Conference on Decision and Control, 2002._ , vol. 3, 2002, pp. 2651–2656 vol.3.
  * [13] E. Najafi, R. Babuška, and G. A. D. Lopes, “A fast sampling method for estimating the domain of attraction,” _Nonlinear Dynamics_ , vol. 86, no. 2, pp. 823–34, Oct. 2016. [Online]. Available: https://doi.org/10.1007/s11071-016-2926-7
  * [14] T. Mancini, I. Melatti, and E. Tronci, “Any-horizon uniform random sampling and enumeration of constrained scenarios for simulation-based formal verification,” _CoRR_ , vol. abs/2109.03330, 2021. [Online]. Available: https://arxiv.org/abs/2109.03330
  * [15] G. Valmórbida, S. Tarbouriech, and G. Garcia, “Region of attraction estimates for polynomial systems,” in _Proceedings of the 48h IEEE Conference on Decision and Control (CDC) held jointly with 2009 28th Chinese Control Conference_ , 2009, pp. 5947–5952.
  * [16] R. Tedrake, I. R. Manchester, M. Tobenkin, and J. W. Roberts, “Lqr-trees: Feedback motion planning via sums-of-squares verification,” _The International Journal of Robotics Research_ , vol. 29, no. 8, pp. 1038–1052, 2010\. [Online]. Available: https://doi.org/10.1177/0278364910369189
  * [17] F. Wiebe, J. Babel, S. Kumar, S. Vyas, D. Harnack, M. Boukheddimi, M. Popescu, and F. Kirchner, “Torque-limited simple pendulum: A toolkit for getting familiar with control algorithms in underactuated robotics,” _Journal of Open Source Software (JOSS)_ , (submitted).
  * [18] G. Birkhoff and G. Rota, _Ordinary Differential Equations_. Wiley, 1991. [Online]. Available: https://books.google.de/books?id=YBjEQgAACAAJ
  * [19] J. C. Robinson, _Higher order linear equations_. Cambridge University Press, 2004, pp. 150–156.
  * [20] A. Platzer, “A complete uniform substitution calculus for differential dynamic logic,” _J. Autom. Reas._ , vol. 59, no. 2, pp. 219–265, 2017.
  * [21] N. Fulton, S. Mitsch, J.-D. Quesel, M. Völp, and A. Platzer, “KeYmaera X: An axiomatic tactical theorem prover for hybrid systems,” in _CADE_ , ser. LNCS, A. P. Felty and A. Middeldorp, Eds., vol. 9195. Springer, 2015, pp. 527–538.