# Charged ultralong-range Rydberg trimers

Daniel J. Bosworth<EMAIL_ADDRESS>Zentrum für Optische
Quantentechnologien, Universität Hamburg,
Luruper Chaussee 149, 22761 Hamburg, Germany
The Hamburg Centre for Ultrafast Imaging, Universität Hamburg,
Luruper Chaussee 149, 22761 Hamburg, Germany
Frederic Hummel<EMAIL_ADDRESS>Max Planck Institute for the Physics of
Complex Systems,
Nöthnitzer Straße 38, 01187 Dresden, Germany  Peter Schmelcher  Zentrum für
Optische Quantentechnologien, Universität Hamburg,
Luruper Chaussee 149, 22761 Hamburg, Germany
The Hamburg Centre for Ultrafast Imaging, Universität Hamburg,
Luruper Chaussee 149, 22761 Hamburg, Germany

###### Abstract

We show that the recently observed class of long-range ion-Rydberg molecules
can be divided into two families of states, which are characterised by their
unique electronic structures resulting from the ion-induced admixture of
quantum defect-split Rydberg $n$P states with different low-field seeking
high-$l$ states. We predict that in both cases these diatomic molecular states
can bind additional ground state atoms lying within the orbit of the Rydberg
electron, thereby forming charged ultralong-range Rydberg molecules (ULRM)
with binding energies similar to that of conventional non-polar ULRM. To
demonstrate this, we consider a Rydberg atom interacting with a single ground
state atom and an ion. The additional atom breaks the system’s cylindrical
symmetry, which leads to mixing between states that would otherwise be
decoupled. The electronic structure is obtained using exact diagonalisation
over a finite basis and the vibrational structure is determined using the
Multi-Configuration Time-Dependent Hartree method. Due to the lobe-like
structure of the electronic density, bound trimers with both linear and
nonlinear geometrical configurations of the three nuclei are possible. The
predicted trimer binding energies and excitation series are distinct enough
from those of the ion-Rydberg dimer to be observed using current experimental
techniques.

##  I Introduction

Hybrid atom-ion systems serve as a test-bed for fundamental quantum physics
research Zipkes  _et al._ (2011); Tomza  _et al._ (2019); Lous and Gerritsma
(2022), enabling studies of cold collisions and chemistry Ratschbacher  _et
al._ (2012); Hall and Willitsch (2012); Härter  _et al._ (2012); Meir  _et
al._ (2016); Pérez-Ríos (2021); Oghittu and Negretti (2022) such as the
formation of cold molecular ions Côté  _et al._ (2002); Schurer  _et al._
(2017); Bosworth  _et al._ (2021). They also provide a platform for precision
measurements Schmid  _et al._ (2010); Zipkes  _et al._ (2010); Veit  _et al._
(2021) and quantum simulation Gerritsma  _et al._ (2012); Schurer  _et al._
(2016). Recent milestones in this field include the first reports of $s$-wave
atom-ion collisions Feldker  _et al._ (2020) and the observation of atom-ion
Feshbach resonances Weckesser  _et al._ (2021).
Over the last two decades, there has also been growing interest in combining
ions with Rydberg atoms in order to engineer atom-ion interactions Secker  _et
al._ (2016); Ewald  _et al._ (2019) and control cold collisions and charge
hopping Secker  _et al._ (2017); Schmid  _et al._ (2018); Côté (2000);
Dieterle  _et al._ (2021). Additionally, an ion-induced Rydberg blockade
effect has been established Engel  _et al._ (2018) and Rydberg atoms have been
used to realise hybrid atom-ion systems in the quantum regime without the need
for an ion trap Kleinbach  _et al._ (2018). Recent theoretical works predicted
bound molecular states between ions and Rydberg states of Rb and Cs Duspayev
_et al._ (2021); Deiß  _et al._ (2021) and their existence and vibrational
dynamics were observed shortly after Zuber  _et al._ (2022); Zou  _et al._
(2022). These ion-Rydberg molecules have bond lengths and energies ranging
from nm to $\mu$m and MHz to GHz, respectively.
Such extreme bonding lengths and energies are seen in another exotic type of
Rydberg molecule formed between a Rydberg atom and one or more ground state
atoms which become bound due to attractive electron scattering. These are
known as ultralong-range Rydberg molecules (ULRM) Eiles (2019); Fey  _et al._
(2020). The importance of electron scattering for describing Rydberg atoms in
atomic gases originated with Fermi Fermi (1934). Fermi’s model was later
applied within the context of ultracold atomic gases, leading to the
prediction of ULRM in 2000 Greene  _et al._ (2000). These molecules were first
observed in 2009 Bendkowsky  _et al._ (2009) and over the past decade they
have been used for probing spatial correlations in ultracold atomic gases
Manthey  _et al._ (2015); Whalen  _et al._ (2019a, b), studying low-energy
electron-atom collisions Anderson  _et al._ (2014); Saßmannshausen  _et al._
(2015); Böttcher  _et al._ (2016) and the formation of Rydberg polarons in the
high-density regime Schmidt  _et al._ (2016); Camargo  _et al._ (2018); Sous
_et al._ (2020); Kleinbach  _et al._ (2018). Furthermore, the formation of
polyatomic ULRM and Rydberg composites has been a topic of major interest
Bendkowsky  _et al._ (2010); Gaj  _et al._ (2014); Liu  _et al._ (2009); Fey
_et al._ (2016); Eiles  _et al._ (2016); Fey  _et al._ (2019a); Eiles  _et
al._ (2020); Hunter  _et al._ (2020) as well as the behaviour of ULRM in
external electric and magnetic fields Lesanovsky  _et al._ (2006); Kurz and
Schmelcher (2014); Krupp  _et al._ (2014); Gaj  _et al._ (2015); Niederprüm
_et al._ (2016); Hummel  _et al._ (2019); Engel  _et al._ (2019); Hummel  _et
al._ (2021a).

Figure 1: Schematic illustration of the two- and three-body systems. (a) a
Rydberg atom (ionic core and valence electron) interacting with an ion. (b) a
Rydberg atom interacting with an ion and a ground state atom.

In light of recent developments, combining ULRM with ions is a natural step
forward: it constitutes a system for studying Rydberg molecules exposed to
inhomogeneous electric fields created by the ion and in which Rydberg atoms
may become bound to both ions and ground state atoms simultaneously. In
contrast to neutral ULRM where the ground state atom is bound within the
Rydberg cloud, ions bind with Rydberg atoms far outside the orbit of the
Rydberg electron. The ion-Rydberg potential wells are formed due to couplings
between neighbouring electronic states, reminiscent of Rydberg macrodimer
potential wells Hollerith  _et al._ (2019).
In this work, we begin by examining the electronic structure of a two-body
system consisting of an ion and a Rydberg atom. We focus in particular on the
electronic density distributions of ion-Rydberg bound states that exist within
different potential wells present in the adiabatic potential energy curves
(PEC). These reveal patterns of pronounced density maxima that are unique to
each potential well due to the ion-induced admixture of different high- and
low-$l$ Rydberg states. Despite their differences, we show that the electronic
density structures can be grouped into two families of states, which differ
primarily in their angular arrangement of the density maxima. We remark that
the Rydberg electron’s probability density in one of these wells has been
discussed previously within the context of the Rydberg atom’s flipping
electric dipole moment around the well minimum and the electronic oscillations
accompanying nuclear dynamics in this well Zuber  _et al._ (2022); Zou  _et
al._ (2022). Building on our considerations of the two-body system’s
electronic structure, we then introduce an additional ground state atom into
the system. We determine adiabatic potential energy surfaces (PES) near the
32P atomic Rydberg state and uncover that their local minima support weakly-
bound trimer states with unique excitation spectra for different geometrical
arrangements of the three species. The properties of the trimer states are
compared against those of the dimer for principle quantum numbers ranging from
17 to 90.
This work is organised as follows. Section II presents an overview of the
current understanding of ion-Rydberg molecules as well as a discussion of
their electronic structure. In section III, we introduce our three-body system
and analyse the corresponding adiabatic PES. Section IV discusses the
resulting vibrational states supported by these surfaces and compares their
properties to those of ion-Rydberg dimers, including possible experimental
aspects. Our conclusions are provided in section V.

##  II Electronic structure of ion-Rydberg dimers

This section first provides some background information on ion-Rydberg
molecules. We then present results for the electronic structure of different
molecular states which motivate the discussion of the three-body system in
sections III and IV.

###  II.1 Squid and snow angel states

Figure 2: Adiabatic potential energy curves and electronic density of the ion-
Rydberg dimer system. Left: adiabatic PEC near the 32P atomic Rydberg state.
The colourbar indicates the $l$-character $\braket{l}$ of the corresponding
electronic states. Energies are given relative to the field-free 32P atomic
Rydberg state. Inset: close-up of well A showing the energies of the first few
vibrational states. Right: the electronic density $|\psi(r,\theta,\phi)|^{2}$
in the $y=0$ plane at the four minima A-D shown in (a). Plots are normalised
such that $\int r^{2}\sin\theta|\psi(r,\theta,\phi)|^{2}drd\theta d\phi=1$.
Note that the ion (not shown) lies outside the Rydberg cloud.

We consider a Rydberg atom in the presence of an ion which lies outside the
Rydberg cloud, as illustrated in Fig. 1 (a).
The systems’s electronic structure is described by the Hamiltonian
$H_{e}=H_{0}+V_{ei}(R_{I})$ 111Unless stated otherwise, atomic units are
assumed throughout.. $H_{0}$ describes the Rydberg electron interacting with
the positively-charged Rydberg core and $V_{ei}(R_{I})$ is the net multipole
interaction between the ion and Rydberg atom for a given internuclear
separation $R_{I}$ Duspayev  _et al._ (2021); Deiß  _et al._ (2021)

$V_{ei}(R_{I})=-\sum_{\lambda=1}^{\infty}\sqrt{\frac{4\pi}{2\lambda+1}}\frac{r^{\lambda}}{R_{I}^{\lambda+1}}Y_{\lambda}^{0}(\theta,\phi).$
(1)

We choose our coordinate system such that the Rydberg core is located at the
origin and the ion-Rydberg internuclear axis is on the $z$-axis. The position
of the Rydberg electron relative to the Rydberg core is given by
$(r,\theta,\phi)$ in spherical coordinates. $Y_{\lambda}^{\mu}$ are the
spherical harmonics with angular momentum $\lambda$ and angular momentum
projection $\mu$, which should not be confused with the orbital angular
momentum quantum numbers defining the Rydberg state, namely $l$ and $m$. Since
we have chosen $\mathbf{R}_{I}=R_{I}\hat{z}$, $\mu$ is restricted to $\mu=0$.
The order of the multipole expansion in (1) is typically truncated at
$\lambda=6$, since higher-order terms only provide energy corrections on the
sub-MHz level which can be safely neglected here.
The time-independent Schrödinger equation for the electronic Hamiltonian reads
$H_{e}\psi_{\nu}(\mathbf{r};R)=\varepsilon_{\nu}(R)\psi_{\nu}(\mathbf{r};R)$,
which depends parametrically on the internuclear separation $R$ whilst $\nu$
labels the separate adiabatic electronic states. The Born-Oppenheimer PEC
$\\{\varepsilon_{\nu}\\}$ obtained from the exact diagonalisation of $H_{e}$
in a finite Rydberg basis are shown in the region of the Rydberg 32P state in
Fig. 2. Note that we neglect the fine and hyperfine structures in our analysis
since they are not responsible for the primary features and results obtained
here.
As predicted in Duspayev  _et al._ (2021); Deiß  _et al._ (2021), the ion-
Rydberg interaction potential $V_{ei}(R_{I})$ couples nearby Rydberg states of
different angular momentum $l$ character, leading to a series of potential
wells in the vicinity of the Rydberg $p$-state which support bound vibrational
states with a spacing on the order of 100 MHz (see inset of Fig. 2). The
colourbar in Fig. 2 encodes the $l$-character of the electronic states, which
changes in particular around the avoided crossings. In principle, these
avoided crossings should introduce non-adiabatic corrections to the Born-
Oppenheimer approximation, yet remarkably the vibrational energies calculated
using the adiabatic approximation are in excellent agreement with current
measurements Zuber  _et al._ (2022). Moreover, a recent theoretical study
determined that the non-adiabatic decay rate of ion-Rydberg molecules should
be far smaller than the radiative decay of the parent Rydberg atom Duspayev
and Raithel (2022).
As a direct result of the admixture of low- and high-$l$ states, the
electronic densities of these molecules display interesting lobe-like
patterns. These can be seen in Fig. 2, which show 2D slices of the Rydberg
electron’s probability density in the $y=0$ plane for the first four potential
wells, marked A to D. The lobes differ in the degree of their azimuthal
localisation. Additionally, the number of lobes present in the probability
density increases for potential wells lying deeper in the fan of electronic
states, thereby forming a series of unique electronic densities. We choose to
classify the series of electronic densities formed by states of $m=0$
character (Fig. 2 A and B) as ’squid states’, from their resemblance to a head
with several appendages, whilst the wing-like features of the density patterns
formed by $|m|=1$ states (Fig. 2 C and D) are reminiscent of a snow angel
pattern.

Figure 3: Adiabatic potential energy surfaces of the squid and snow angel
states for the trimer system. (a) shows the electronic density of the three-
legged squid state in the dimer system. (b) is a 2-D slice taken at
$R_{I}=0.510~{}\mu$m of the squid state’s PES with major local minima denoted
by {$\text{w}_{i}$}. (c) shows the electronic density of the snow angel state
in the $y=0$ plane. (d) is a 2-D slice taken at $R_{I}=0.488~{}\mu$m of the
snow angel state’s PES with major local minima denoted by
{$\text{w}^{\prime}_{i}$}. Energies are given relative to the field-free 32P
atomic Rydberg state.

###  II.2 From the charged dimer to the trimer

The interaction of a Rydberg atom and a ground state atom is determined by the
scattering of the highly-excited Rydberg electron off the ground state atom.
The corresponding $s$-wave scattering term is described by the Fermi
pseudopotential $V_{ea}=2\pi
a_{s}[k(R_{A})]\delta(\mathbf{r}-\mathbf{R_{A}})$. This model predicts the
appearance of wells in the PEC for regions where the scattering length
$a_{s}[k(R_{A})]$ is negative. These wells support weakly-bound vibrational
states with binding energies in the range of MHz to GHz as first shown in
Greene  _et al._ (2000). These bound molecular states are ultralong-range
Rydberg molecules (ULRM) and examples include the high angular momentum polar
trilobite state and the low angular momentum non-polar $s$-state Eiles (2019).
Motivated by this fact, we now turn to the question of whether the same
mechanism could enable ion-Rydberg molecules to bind an additional ground
state atom within the orbit of the Rydberg electron, forming a charged ULRM
trimer. To that end, we consider an additional ground state atom within the
Rydberg cloud, as shown by the schematic in Fig. 1 (b). If this binding is
possible, the various lobes of high electron density present in the squid and
snow angel states would enable the formation of both linear and nonlinear
bound configurations of the atoms, akin to the $d$-state angular trimers
reported in Fey  _et al._ (2019a).

##  III Electronic structure of the trimer

In this section, we consider the effect of an additional ground state atom on
the first squid and snow angel states, whose unperturbed electronic densities
are shown in Fig. 2 A and C, respectively.

###  III.1 Setup and interactions

With the Rydberg core at the origin of our coordinate system, the ion and
ground state atom are located at positions $(R_{I},0,0)$ and
$(R_{A},\Theta,0)$ respectively. A schematic for this three-body system is
provided in Fig. 1 (b). Within the Born-Oppenheimer approximation, the
electronic structure of the Rb+-Rb∗-Rb system is based on the following
Hamiltonian

$H_{e}=H_{0}+V_{ei}(R_{I})+V_{ea}(R_{A})+V_{ca}(R_{A})+V_{ia}(R_{I},R_{A},\Theta).$
(2)

$V_{ei}$ describes the electron-ion interaction (1) and $V_{ea}$ represents
the $s$-wave scattering between the electron and the ground state atom.
$V_{ca}$ and $V_{ia}$ are the interactions between the Rydberg core and the
ion with the neutral atom, respectively. These interactions take the form of a
charge-induced dipole (polarisation potential) interaction, e.g.
$V_{ca}\propto-1/R_{A}^{4}$.
We obtain the adiabatic PES $\\{\varepsilon_{\nu}(R_{I},R_{A},\Theta)\\}$ from
the exact diagonalisation of the electronic Hamiltonian (2) using a finite
basis of unperturbed Rydberg states $\\{\ket{n,l,m}\\}$, where the quantum
numbers take their usual meaning. Specifically, our electronic basis includes
all states enclosed by the nearest six hydrogenic manifolds centred around the
32P atomic Rydberg state up to a maximum magnetic quantum number of $|m|=4$.
Besides introducing additional interactions, the presence of the ground state
atom breaks the dimer’s cylindrical symmetry, such that $m$ is no longer a
good quantum number. This means couplings between states of different
$m$-character are no longer prohibited, though generally this will only be
relevant where states become near-degenerate. One example of this occurs close
to the minimum of the first snow angel state, labelled as C in Fig. 2. In
principle, this requires that we enlarge our basis to account for the
additional couplings that arise between different $m$-states.

###  III.2 Adiabatic potential energy surfaces

Let us analyse in the following the relevant adiabatic PES of the trimer.
Figure 3 shows the electronic densities of the squid and snow angel states in
the dimer system alongside 2D slices of the corresponding Born-Oppenheimer PES
$\varepsilon(R_{I},R_{A},\Theta)$ for the trimer system at fixed $R_{I}$. Each
lobe in the electronic density leads to a unique local minimum along $\Theta$.
Along $R_{A}$, the surfaces exhibit a series of local minima due to the
oscillatory electronic density along each lobe. The deepest minima are found
at $R_{A}\approx 1700~{}a_{0}$ and have a depth on the order of 100 MHz,
similar to a conventional non-polar ULRM at this principle quantum number.
For the squid state we focus on three local minima labelled w1, w2 and w3 in
Fig. 3 (b). These minima occur at angles $\Theta=0$, $\Theta=0.8\pi$ and
$\Theta=\pi$, respectively which map directly to the three unique lobes in the
electronic density in Fig. 3 (a). Well w3 is the deepest of these, which is
not surprising since the electronic probability density is highly-localised at
this position. For the snow angel state, we focus on the two minima labelled
w${}^{\prime}_{1}$ and w${}^{\prime}_{2}$. The depths of these wells are
similar in magnitude to those of the squid state. Since the electronic
probability density vanishes along the $z$-axis, no local minima appear along
$\Theta=0$ and $\Theta=\pi$ and as such it should not be possible for the snow
angel to support linear trimer configurations.

###  III.3 Coupling among states of different magnetic quantum number

Figure 4: $m$-admixture in the three-body squid state. (a) angular slice of
the squid state’s PES with colourbar denoting the electronic state’s
$m$-character $\braket{}{m}{}$. (b)-(d) radial slices of the PES near the w1,
w2 and w3 local minima, respectively. (e)-(g) the squid states’ composition in
terms of Rydberg basis states grouped by their $m$ quantum number at the w1,
w2 and w3 local minima, respectively. Energies are given relative to the
field-free 32P Rydberg state. Figure 5: $m$-admixture in the three-body snow
angel state. (a) angular slice of the snow angel state’s PES with colourbar
denoting the electronic state’s $m$-character $\braket{m}$. (b),(c) radial
slices of the PES near the w${}^{\prime}_{1}$ and w${}^{\prime}_{2}$ local
minima, respectively. (d),(f) the snow angel state’s composition in terms of
Rydberg basis states grouped by their $m$ quantum number at the w1 and w2
local minima, respectively. Energies are given relative to the field-free 32P
Rydberg state.

Let us now examine couplings among PES of different magnetic quantum number
$m$ that arise due to the presence of the ground state atom. Fig. 4 (a) shows
an angular cut of the squid state’s PES $\varepsilon(\Theta)$ for
$R_{I}=9756~{}a_{0}$ and $R_{A}=1653~{}a_{0}$, with a colourbar denoting the
$m$-character of the state. The surface exhibits potential wells at
$\Theta=0$, $\Theta=0.8\pi$ and $\Theta=\pi$, which correspond to the local
minima $\text{w}_{1}$, $\text{w}_{2}$ and $\text{w}_{3}$ defined in Fig. 3
(b), respectively. Along this cut, the $m$-character is pure $m=0$ and no
other surfaces are present in the energy window shown. However, as is visible
from Fig. 2, the ion-Rydberg dimer’s potential well exists above a fan of
high-$l$ Stark-split states. Fig. 4 (b)-(d) show slices of the PES along
$R_{I}$ for values of $\Theta$ and $R_{A}$ taken around the three local minima
in Fig. 4 (a). For wells $\text{w}_{1}$ and $\text{w}_{2}$, the PES are
qualitatively similar to those of the dimer and show only marginal deviation
from $m=0$ character. This is further illustrated by Fig. 4 (e) and (f), which
show the cumulative overlap coefficients
$|c_{i}|^{2}=|\braket{\psi(R_{I},R_{A},\Theta)}{\varphi_{i}}|^{2}$ of the
squid state $\psi(R_{I},R_{A},\Theta)$ for each local minimum with atomic
Rydberg basis states $\\{\varphi_{i}\\}$, grouped by their magnetic quantum
number $m$. The lack of coupling is perhaps most surprising for the case of
well w2, which corresponds to a non-linear configuration of the three atoms
and is hence the least symmetric of all three wells. In contrast, the slice
along $R_{I}$ for well $\text{w}_{3}$ in Fig. 4 (d) shows that the local
minimum is shifted to slightly lower $R_{I}$, such that it crosses into the
fan of high-$l$ states. This observed shift does not seem to modify the
$l$-character, but it does result in an admixture with states of finite $m$ at
small angles around the local minimum $\Theta=\pi+\epsilon$. This admixture is
primarily with $|m|=1$ states, as shown in Fig. 4 (g). Note that these
couplings vanish at $\Theta=\pi$, since cylindrical symmetry is restored, i.e.
the three atoms are co-linear.
As discussed already in section III.1 in the context of the dimer system, the
squid state’s PEC crosses the curve of the snow angel state near the snow
angel’s potential well (cf. well C in Fig. 2). In the dimer system, this
crossing is exact since there can be no coupling among surfaces of different
magnetic quantum number $m$. In the trimer system, this no longer holds and a
finite coupling between the $m=0$ squid surface and the $|m|=1$ snow angel
surfaces is observed in Fig. 5. This admixture is most significant near the
w${}^{\prime}_{2}$ well, for which a pronounced avoided crossing arises, as
shown in Figure 5 (a) and (c). Whilst the shallower w${}^{\prime}_{1}$ well
also shows minor admixture (Figure 5 (d)), it is orders of magnitude smaller
than for w${}^{\prime}_{2}$ and if any avoided crossing is present, its gap
cannot be discerned on the energy scale of Figure 5 (b).
From here on, we limit our analysis of the vibrational structure to those
local minima which show only relatively weak $m$-admixture (w1, w2 and
w${}^{\prime}_{1}$) because we expect the Born-Oppenheimer approximation to
continue to hold for these. We treat crossings with states of different
$m$-character near these minima as being approximately exact, though this only
applies to w${}^{\prime}_{1}$, since no crossings were observed in the
vicinity of the minima of wells w1 and w2. For w3 and w${}^{\prime}_{2}$, we
expect strong non-adiabatic couplings to exist with neighbouring PES.

##  IV Vibrational structure of the trimer

In this section, we explore the vibrational structure of the three-body system
defined by the nuclear Hamiltonian $H_{n}$ in the vicinity of the local minima
w1, w2 and w${}^{\prime}_{1}$ of the electronic squid and snow angel states.
Let us first discuss our computational approach and consequently the analysis
of the vibrational structure.

###  IV.1 Methodology and computational approach

The total Hamiltonian of our trimer system reads $H=H_{e}+H_{n}$. In the
previous section, we have discussed solutions to the electronic Hamiltonian
$H_{e}$ obtained via exact diagonalisation in a finite basis of Rydberg states
$\\{\ket{n,l,m}\\}$. Here, we focus on the vibrational motion of the nuclei,
assuming that the system possesses zero angular momentum $J=0$. Accordingly,
the full vibrational Hamiltonian Fey  _et al._ (2019b); Carter and Handy
(1982); Handy (1987) is given by

$\begin{split}H_{n}&=\frac{1}{m}\bigg{[}-\frac{\partial^{2}}{\partial
R_{I}^{2}}-\frac{\partial^{2}}{\partial
R_{A}^{2}}-\cos(\Theta)\frac{\partial}{\partial R_{I}}\frac{\partial}{\partial
R_{A}}\bigg{]}\\\
&-\frac{1}{m}\bigg{(}\frac{1}{R_{I}^{2}}+\frac{1}{R_{A}^{2}}-\frac{\cos(\Theta)}{R_{I}R_{A}}\bigg{)}\bigg{(}\frac{\partial^{2}}{\partial\Theta^{2}}+\cot(\Theta)\frac{\partial}{\partial\Theta}\bigg{)}\\\
&-\frac{1}{m}\bigg{(}\frac{1}{R_{I}R_{A}}-\frac{1}{R_{A}}\frac{\partial}{\partial
R_{I}}-\frac{1}{R_{I}}\frac{\partial}{\partial
R_{A}}\bigg{)}\bigg{(}\cos(\Theta)+\sin(\Theta)\frac{\partial}{\partial\Theta}\bigg{)}\\\
&+\varepsilon_{\nu}(R_{I},R_{A},\Theta),\\\ \end{split}$ (3)

where $m$ is the atomic mass of Rb and $\varepsilon_{\nu}(R_{I},R_{A},\Theta)$
is the $\nu^{\text{th}}$ adiabatic electronic PES satisfying
$H_{e}\psi_{\nu}(\mathbf{r};R_{I},R_{A},\Theta)=\varepsilon_{\nu}(R_{I},R_{A},\Theta)\psi_{\nu}(\mathbf{r};R_{I},R_{A},\Theta)$.
The time-independent Schrödinger equation for the vibrational Hamiltonian in
the Born-Oppenheimer approximation reads

$H_{n}\chi_{i}(R_{I},R_{A},\Theta)=E_{i}\chi_{i}(R_{I},R_{A},\Theta),$ (4)

with vibrational eigenstates $\\{\chi_{i}(R_{I},R_{A},\Theta)\\}$ and
associated eigenenergies $E_{i}$.
In order to efficiently solve the above vibrational problem, we utilise the
powerful Multi-Configuration Time-Dependent Hartree method (MCTDH) Meyer  _et
al._ (1990); Beck  _et al._ (2000); Meyer and Worth (2003); Meyer  _et al._
(2006); Doriol  _et al._ (2008); Meyer (2012); Worth  _et al._ . A thorough
introduction to MCTDH can be found in Meyer  _et al._ (2006); Meyer (2012). In
the following we provide a brief account of the approach in order to be self-
contained. MCTDH is an ab initio method for multi-mode wavepacket propagation
in high-dimensional spaces. The MCTDH representation of our vibrational
wavefunction $\chi(R_{I},R_{A},\Theta,t)$ is written as a series of Hartree
products:

$\begin{split}\chi(R_{I},R_{A},\Theta,t)&=\sum_{i_{1}=1}^{n_{1}}\sum_{i_{2}=1}^{n_{2}}\sum_{i_{3}=1}^{n_{3}}A_{i_{1},i_{2},i_{3}}(t)\\\
&\times\varphi_{i_{1}}^{(1)}(R_{I},t)\varphi_{i_{2}}^{(2)}(R_{A},t)\varphi_{i_{3}}^{(3)}(\Theta,t),\end{split}$
(5)

where $A_{i_{1},i_{2},i_{3}}(t)$ are time-dependent coefficients and
$\\{\varphi_{i_{d}}^{(d)}\\}_{i_{d}=1}^{n_{d}}$ are the so-called single
particles functions (SPFs) for the $d^{\text{th}}$ degree of freedom, for
which a total of $n_{d}$ SPFs are used. MCTDH reduces computational effort by
employing a small time-dependent basis that evolves according to the Dirac-
Frenkel variational principle
$\braket{\delta\chi}{(i\partial_{t}-\hat{H})}{\chi}=0$. This ensures that the
basis follows the active part of the complete Hilbert space as it evolves over
time. The time-depdendent SPFs are described using a time-independent discrete
variable representation (DVR) Harris  _et al._ (1965). In this work, we use
sine DVRs for the radial degrees of freedom and a Legendre DVR for the angular
degree of freedom. SPFs for the radial degrees of freedom are normalised as
$\int dR|\varphi(R)|^{2}=1$ and the angular SPFs are normalised according to
$\int d\Theta\sin\Theta|\varphi(\Theta)|^{2}=1$.
The vibrational ground-state can be obtained by propagating the starting
wavepacket in imaginary time, whereas excited states are obtained using so-
called improved relaxation. Here, the time-dependent coefficients
$\\{A_{i_{1},i_{2},i_{3}}(t)\\}$ are set equal to an eigenvector $\mathbf{A}$
of the Hamiltonian in the instantaneous SPF basis
$\\{\varphi_{i_{d}}^{(d)}\\}_{i_{d}=1}^{n_{d}}$. The coefficients are then
kept constant whilst the SPFs are relaxed using imaginary time propagation.
This process is repeated until $\chi$ converges to a stationary state of the
Hamiltonian. For our analysis, we seek the lowest few eigenstates in the
potential wells w1, w2 and w${}^{\prime}_{1}$ of the squid and snow angel PES.
These are obtained by first performing a block improved relaxation scheme
Meyer (2012), which provides approximate results for the first 30-40 lowest
energy eigenstates. Out of this set, we select the states of interest and
relax each of them individually using improved relaxation until convergence is
achieved.
We consider a calculation to be converged when the change in energy between
time steps remains consistently less than 1 Hz over a period of 100 time
steps. Additionally, we assure that the occupation of the $n_{d}^{\text{th}}$
orbital for each degree of freedom is less than $0.1\%$ and that the
occupations $\\{A_{i_{1},i_{2},i_{3}}(t)\\}$ decrease exponentially. For the
majority of our calculations, $100$ grid points and $10$ single particle
functions for each degree of freedom were sufficient to converge all
eigenenergies.

###  IV.2 Vibrational structure

Figure 6: Selection of vibrational states of the squid state trimer’s w1 potential well. (a) angular cut through the potential energy surface (grey) near the local minimum $(R_{I},R_{A},\Theta)=(9782a_{0},1657a_{0},0)$. (b) reduced angular densities $\rho_{\nu}(\Theta)=\int dR_{I}dR_{A}|\chi_{\nu}(R_{I},R_{A},\Theta)|^{2}$ of the vibrational states $\nu=\\{0,1,2,3,10\\}$ (filled curves) offset by their energy (dashed lines). (c) and (d) show reduced radial densities $\rho_{\nu}(R_{I},R_{A})=\int d\Theta|\chi_{\nu}(R_{I},R_{A},\Theta)|^{2}$ of $\chi_{1}$, $\chi_{3}$ and $\chi_{10}$, respectively. The eigenstates are normalised as $\int dR_{I}dR_{A}d\Theta|\chi_{\nu}(R_{I},R_{A},\Theta)|^{2}=1$. All energies given relative to the field-free 32P atomic Rydberg state. |  $\Delta\braket{R_{I}}$ ($a_{0}$) |  $\Delta\braket{d_{z}}$ (Debye) |  $\Delta E_{000}$ (MHz) |  $\tilde{E}_{001}$ (MHz) |  $\tilde{E}_{010}$ (MHz) |  $\tilde{E}_{100}$ (MHz)
---|---|---|---|---|---|---
Squid $\text{w}_{1}$ | 14.1 | 18.2 | -181.9 | 21.2 | 39.6 | 91.9 (88.1)
Squid $\text{w}_{2}$ | -5.3 | 20.8 | -102.8 | 7.4 | 31.1 | 92.6 (88.1)
Angel $\text{w}^{\prime}_{1}$ | 5.9 | 19.3 | -79.9 | 3.0 | 29.3 | 86.6 (86.6)
Table 1: Contrasting properties of the dimer and trimer vibrational states.
From left to right, the first three columns give the difference in the
vibrational ground state’s (I) expected ion-Rydberg separation
$\Delta\braket{R_{I}}=\braket{R_{I}^{\text{trimer}}}-\braket{R_{I}^{\text{dimer}}}$,
(II) expected electric dipole moment of the Rydberg atom
$\Delta\braket{d_{z}}=\braket{d_{z}^{\text{trimer}}}-\braket{d_{z}^{\text{dimer}}}$
and (III) energy $\Delta
E_{000}=E_{000}^{\text{trimer}}-E_{000}^{\text{dimer}}$. The last three
columns contain the trimer’s excitation energies wrt. the trimer ground state
of the first (IV) $\Theta$ bending mode $\tilde{E}_{001}=E_{001}-E_{000}$, (V)
$R_{A}$ stretching mode $\tilde{E}_{010}=E_{010}-E_{000}$ and (VI) $R_{I}$
stretching mode $\tilde{E}_{100}=E_{100}-E_{000}$. Where applicable,
corresponding values for the dimer are given in brackets. All quantities are
rounded to the nearest decimal place. Figure 7: (a) Scaling of the squid
dimer’s bond length with principle quantum number $n$. (b) difference in ion-
Rydberg bond length and (c) potential well minima between the squid dimer and
squid trimer for the w1 potential well (cf. well Fig. 3).

Fig. 6 shows results for the reduced angular and radial densities of the four
lowest vibrational states in the squid state’s w1 potential well. Table 1
summarises properties of the lowest vibrational states in all three trimer
wells and compares these with the vibrational states of the squid and snow
angel dimer.
In Fig. 6 (a) and (b), we observe that the energetically lowest states in the
potential well w1 are localised around $\Theta=0$ such that the three atoms
are arranged along a common internuclear axis, i.e. in a collinear
configuration. The well’s ground state is approximately 180 MHz lower in
energy than the ground state of the squid state dimer and exhibits a positive
shift of the ion-Rydberg bond length and corresponding electric dipole moment,
as shown in Table 1. From Fig. 6 (b), we find that the reduced angular density
of the first excited state exhibits an additional node, indicating that this
is a bending mode. The excitation energy of this mode is about 20 MHz above
the ground state energy. In contrast, the first excited state of the dimer is
approximately 88.1 MHz. The second excited state is also a bending mode, with
an additional node along the reduced angular density in Fig. 6 (b). The third
excited state’s energy lies close to that of the second excited state, yet its
reduced radial density in Fig. 6 (d) confirms that it is an excitation in
$R_{A}$, a stretching mode. The energy gap to the first excitation in $R_{I}$
is 91.6 MHz, which is of the same order of magnitude as the energy gap between
vibrational states of the dimer.
We now examine some global trends of properties among the three wells evident
from Table 1. Surprisingly, the ground state atom’s presence appears to
contribute to a consistent positive shift of 20 Debye in the Rydberg atom’s
electric dipole moment $\braket{d_{z}}$, independent of which well it inhabits
and regardless of whether the expected ion-Rydberg binding length
$\braket{R_{I}}$ has increased or decreased with respect to the dimer. On the
other hand, the energy gap to the first $R_{I}$ stretching mode is essentially
unchanged by the presence of the perturber. The energy gaps of the bending and
$R_{A}$ stretching modes do vary among the different potential wells. In
particular, in the w${}^{\prime}_{1}$ well 3 MHz are necessary to excite its
first bending mode, whilst the same gap in the w1 well is an order of
magnitude larger. From the 2D graphs of the PES in Fig. 3 (b) and (d), we see
that w1 and w${}^{\prime}_{1}$ are reasonably extended along $\Theta$ which
implies a weak angular confinement in both cases. Nevertheless, the w1 well is
considerably deeper than the w${}^{\prime}_{1}$ well which leads the greater
excitation gap.
In an experimental setting, it should be possible to distinguish trimer states
from dimer states due to the trimer’s deeper binding energies and
distinguishable excitation series for bending and stretching modes depending
on the geometric configuration of the nuclei. Moreover, mass spectroscopy
techniques such as that employed in Zuber  _et al._ (2022) could be used to
distinguish the two based on their mass difference. In principle, the trimer’s
shift in ion-Rydberg binding length $R_{\text{min}}$ could serve as a further
indicator for the presence of trimer states. The shifts predicted in Table 1
are however on the order of a few Bohr radii, which is clearly too small to be
observed with current techniques. As first shown in Deiß  _et al._ (2021), the
equilibrium separation in the ion-Rydberg dimer $R_{\text{min}}$ increases
with principle quantum number $n$ and in Fig. 7 (a) we find this scaling to be
$R_{\text{min}}\propto n^{2.80\pm 0.18}$. However, we see no clear trend of
$\Delta
R_{\text{min}}=R_{\text{min}}^{\text{trimer}}-R_{\text{min}}^{\text{dimer}}$
for the well w1 given in Fig. 7 (b).
On the other hand, there is a clear trend in the difference of the local
energy minima
$\Delta\epsilon_{\text{min}}=\epsilon_{\text{min}}^{\text{dimer}}-\epsilon_{\text{min}}^{\text{trimer}}$
between the PES of the dimer and trimer. Fig. 7 (c) shows
$\Delta\epsilon_{\text{min}}$ for the trimer w1 well and we find
$\Delta\epsilon_{\text{min}}\propto n^{-5.97\pm 0.06}$. This can be understood
as follows. To first order, the interaction of the ground state atom with the
Rydberg electron contributes an energy shift proportional to the electron’s
probability density $\Delta\epsilon\propto|\psi(R_{A})|^{2}$, whose sign
depends on the s-wave scattering length $a_{s}[k(R_{A})]$. The electron’s
probability density scales with the radius of the Rydberg orbit as
$|\psi(r)|^{2}\propto 1/r^{3}$, which in turn scales quadratically with $n$.
This gives $\Delta\epsilon\propto n^{-6}$, which is close to our fitted value.
Since $\Delta\epsilon_{\text{min}}$ decays with $n$, the trimer’s vibrational
spectrum will lie energetically closer to the dimer’s at higher principle
quantum numbers, making them harder to distinguish. By modulating the density
of the background gas of atoms, signals of trimer states can be suppressed or
enhanced.

##  V Conclusions

In this work, we first studied the case of an ion interacting with a Rydberg
atom and determined the electronic structure of bound states which exist in
various potential wells in the spectrum of adiabatic potential energy curves.
For atomic species like Rb and Cs, the ion’s presence mixes the quantum
defect-split $p$-state with its neighbouring quasi-degenerate manifold of
high-$l$ states. This results in the formation of a twofold series of
potential wells, one for each $|m|$ sub-level of the $p$-state. We found that
the Rydberg electron’s probability density exhibits several lobes, which vary
in their size and number between the wells. Due to their distinctive shapes,
we classify the electronic density patterns along the $m=0$ series as squid
states, whilst those along the $|m|=1$ series are termed snow angel states.
We predict that the lobes in the electronic density admit the binding of
additional ground state atoms via the same binding mechanism as conventional
polar and non-polar ULRM. We have demonstrated this for a three-body system of
a Rydberg atom interacting with a ground state atom and a single ion, in the
limit of non-overlapping charge distributions. Bound vibrational states exist
among the three species, forming a charged ultralong-range Rydberg trimer. The
ground state atom can become bound within different lobes in the electronic
density, such that linear as well as non-linear geometric configurations of
the nuclear framework are possible.
The ground state atom leads to a significant shift in the energy of the
trimer’s vibrational ground state and the additional nuclear degrees of
freedom give rise to distinct excitation series for bending and stretching
modes, which vary between the local minima associated with different lobes in
the electronic structure. Using current spectroscopic techniques, it should
thus be possible to distinguish the trimer and dimer states. Similarly to the
ion-Rydberg dimer, we expect nonadiabatic couplings between neighbouring
potential energy surfaces to be a significant channel for molecular decay in
our system. Indeed, the effect may even be stronger in the trimer system since
crossings between states of different $m$ are no longer symmetry-protected.
Corresponding investigations are left to future studies.
Our work lays the foundation for exploring ULRM in inhomogeneous electric
fields. Future work could consider the interaction of an ion with polyatomic
ULRM and wavepacket dynamics in the ion’s field. Additionally, conical
intersections can occur between adiabatic potential energy surfaces in
triatomic systems and thus one could examine the role of beyond Born-
Oppenheimer physics in our system in the spirit of recent work in this
direction Hummel  _et al._ (2021b).

## Acknowledgements

D. J. B. thanks Maxim Pyzh for many fruitful discussions. This work is funded
by the Cluster of Excellence “Advanced Imaging of Matter” of the Deutsche
Forschungsgemeinschaft (DFG)-EXC 2056, Project ID No. 390715994.

## References

  * Zipkes  _et al._ (2011) C. Zipkes, L. Ratschbacher, S. Palzer, C. Sias, and M. Köhl, J. Phys.: Conf. Ser. 264, 012019 (2011).
  * Tomza  _et al._ (2019) M. Tomza, K. Jachymski, R. Gerritsma, A. Negretti, T. Calarco, Z. Idziaszek, and P. S. Julienne, Rev. Mod. Phys. 91, 035001 (2019).
  * Lous and Gerritsma (2022) R. S. Lous and R. Gerritsma, in  _Advances In Atomic, Molecular, and Optical Physics_, Advances in Atomic, Molecular, and Optical Physics, Vol. 71, edited by L. F. DiMauro, H. Perrin, and S. F. Yelin (Academic Press, 2022) pp. 65–133.
  * Ratschbacher  _et al._ (2012) L. Ratschbacher, C. Zipkes, C. Sias, and M. Köhl, Nature Phys 8, 649 (2012).
  * Hall and Willitsch (2012) F. H. J. Hall and S. Willitsch, Phys. Rev. Lett. 109, 233202 (2012).
  * Härter  _et al._ (2012) A. Härter, A. Krükow, A. Brunner, W. Schnitzler, S. Schmid, and J. H. Denschlag, Phys. Rev. Lett. 109, 123201 (2012).
  * Meir  _et al._ (2016) Z. Meir, T. Sikorsky, R. Ben-shlomi, N. Akerman, Y. Dallal, and R. Ozeri, Phys. Rev. Lett. 117, 243401 (2016).
  * Pérez-Ríos (2021) J. Pérez-Ríos, Mol. Phys. 119, e1881637 (2021).
  * Oghittu and Negretti (2022) L. Oghittu and A. Negretti, Phys. Rev. Research 4, 023069 (2022).
  * Côté  _et al._ (2002) R. Côté, V. Kharchenko, and M. D. Lukin, Phys. Rev. Lett. 89, 093001 (2002).
  * Schurer  _et al._ (2017) J. M. Schurer, A. Negretti, and P. Schmelcher, Phys. Rev. Lett. 119, 063001 (2017).
  * Bosworth  _et al._ (2021) D. J. Bosworth, M. Pyzh, and P. Schmelcher, Phys. Rev. A 103, 033303 (2021).
  * Schmid  _et al._ (2010) S. Schmid, A. Härter, and J. H. Denschlag, Phys. Rev. Lett. 105, 133202 (2010).
  * Zipkes  _et al._ (2010) C. Zipkes, S. Palzer, C. Sias, and M. Köhl, Nature 464, 388 (2010).
  * Veit  _et al._ (2021) C. Veit, N. Zuber, O. A. Herrera-Sancho, V. S. V. Anasuri, T. Schmid, F. Meinert, R. Löw, and T. Pfau, Phys. Rev. X 11, 011036 (2021).
  * Gerritsma  _et al._ (2012) R. Gerritsma, A. Negretti, H. Doerk, Z. Idziaszek, T. Calarco, and F. Schmidt-Kaler, Phys. Rev. Lett. 109, 080402 (2012).
  * Schurer  _et al._ (2016) J. M. Schurer, R. Gerritsma, P. Schmelcher, and A. Negretti, Phys. Rev. A 93, 063602 (2016).
  * Feldker  _et al._ (2020) T. Feldker, H. Fürst, H. Hirzler, N. V. Ewald, M. Mazzanti, D. Wiater, M. Tomza, and R. Gerritsma, Nat. Phys. 16, 413 (2020).
  * Weckesser  _et al._ (2021) P. Weckesser, F. Thielemann, D. Wiater, A. Wojciechowska, L. Karpa, K. Jachymski, M. Tomza, T. Walker, and T. Schaetz, Nature 600, 429 (2021).
  * Secker  _et al._ (2016) T. Secker, R. Gerritsma, A. W. Glaetzle, and A. Negretti, Phys. Rev. A 94, 013420 (2016).
  * Ewald  _et al._ (2019) N. V. Ewald, T. Feldker, H. Hirzler, H. A. Fürst, and R. Gerritsma, Phys. Rev. Lett. 122, 253401 (2019).
  * Secker  _et al._ (2017) T. Secker, N. Ewald, J. Joger, H. Fürst, T. Feldker, and R. Gerritsma, Phys. Rev. Lett. 118, 263201 (2017).
  * Schmid  _et al._ (2018) T. Schmid, C. Veit, N. Zuber, R. Löw, T. Pfau, M. Tarana, and M. Tomza, Phys. Rev. Lett. 120, 153401 (2018).
  * Côté (2000) R. Côté, Phys. Rev. Lett. 85, 5316 (2000).
  * Dieterle  _et al._ (2021) T. Dieterle, M. Berngruber, C. Hölzl, R. Löw, K. Jachymski, T. Pfau, and F. Meinert, Phys. Rev. Lett. 126, 033401 (2021).
  * Engel  _et al._ (2018) F. Engel, T. Dieterle, T. Schmid, C. Tomschitz, C. Veit, N. Zuber, R. Löw, T. Pfau, and F. Meinert, Phys. Rev. Lett. 121, 193401 (2018).
  * Kleinbach  _et al._ (2018) K. S. Kleinbach, F. Engel, T. Dieterle, R. Löw, T. Pfau, and F. Meinert, Phys. Rev. Lett. 120, 193401 (2018).
  * Duspayev  _et al._ (2021) A. Duspayev, X. Han, M. A. Viray, L. Ma, J. Zhao, and G. Raithel, Phys. Rev. Research 3, 023114 (2021).
  * Deiß  _et al._ (2021) M. Deiß, S. Haze, and J. Hecker Denschlag, Atoms 9, 34 (2021).
  * Zuber  _et al._ (2022) N. Zuber, V. S. V. Anasuri, M. Berngruber, Y.-Q. Zou, F. Meinert, R. Löw, and T. Pfau, Nature 605, 453 (2022).
  * Zou  _et al._ (2022) Y.-Q. Zou, M. Berngruber, V. S. V. Anasuri, N. Zuber, F. Meinert, R. Löw, and T. Pfau, (2022), arXiv:2208.07776 [cond-mat, physics:physics] .
  * Eiles (2019) M. T. Eiles, J. Phys. B: At. Mol. Opt. Phys. 52, 113001 (2019).
  * Fey  _et al._ (2020) C. Fey, F. Hummel, and P. Schmelcher, Mol. Phys. 118, e1679401 (2020).
  * Fermi (1934) E. Fermi, Nuovo Cim 11, 157 (1934).
  * Greene  _et al._ (2000) C. H. Greene, A. S. Dickinson, and H. R. Sadeghpour, Phys. Rev. Lett. 85, 2458 (2000).
  * Bendkowsky  _et al._ (2009) V. Bendkowsky, B. Butscher, J. Nipper, J. P. Shaffer, R. Löw, and T. Pfau, Nature 458, 1005 (2009).
  * Manthey  _et al._ (2015) T. Manthey, T. Niederprüm, O. Thomas, and H. Ott, New J. Phys. 17, 103024 (2015).
  * Whalen  _et al._ (2019a) J. D. Whalen, S. K. Kanungo, R. Ding, M. Wagner, R. Schmidt, H. R. Sadeghpour, S. Yoshida, J. Burgdörfer, F. B. Dunning, and T. C. Killian, Phys. Rev. A 100, 011402 (2019a).
  * Whalen  _et al._ (2019b) J. D. Whalen, R. Ding, S. K. Kanungo, T. C. Killian, S. Yoshida, J. Burgdörfer, and F. B. Dunning, Mol. Phys. 117, 3088 (2019b).
  * Anderson  _et al._ (2014) D. A. Anderson, S. A. Miller, and G. Raithel, Phys. Rev. Lett. 112, 163201 (2014).
  * Saßmannshausen  _et al._ (2015) H. Saßmannshausen, F. Merkt, and J. Deiglmayr, Phys. Rev. Lett. 114, 133201 (2015).
  * Böttcher  _et al._ (2016) F. Böttcher, A. Gaj, K. M. Westphal, M. Schlagmüller, K. S. Kleinbach, R. Löw, T. C. Liebisch, T. Pfau, and S. Hofferberth, Phys. Rev. A 93, 032512 (2016).
  * Schmidt  _et al._ (2016) R. Schmidt, H. R. Sadeghpour, and E. Demler, Phys. Rev. Lett. 116, 105302 (2016).
  * Camargo  _et al._ (2018) F. Camargo, R. Schmidt, J. D. Whalen, R. Ding, G. Woehl, S. Yoshida, J. Burgdörfer, F. B. Dunning, H. R. Sadeghpour, E. Demler, and T. C. Killian, Phys. Rev. Lett. 120, 083401 (2018).
  * Sous  _et al._ (2020) J. Sous, H. R. Sadeghpour, T. C. Killian, E. Demler, and R. Schmidt, Phys. Rev. Research 2, 023021 (2020).
  * Bendkowsky  _et al._ (2010) V. Bendkowsky, B. Butscher, J. Nipper, J. B. Balewski, J. P. Shaffer, R. Löw, T. Pfau, W. Li, J. Stanojevic, T. Pohl, and J. M. Rost, Phys. Rev. Lett. 105, 163201 (2010).
  * Gaj  _et al._ (2014) A. Gaj, A. T. Krupp, J. B. Balewski, R. Löw, S. Hofferberth, and T. Pfau, Nat. Commun. 5, 4546 (2014).
  * Liu  _et al._ (2009) I. C. H. Liu, J. Stanojevic, and J. M. Rost, Phys. Rev. Lett. 102, 173001 (2009).
  * Fey  _et al._ (2016) C. Fey, M. Kurz, and P. Schmelcher, Phys. Rev. A 94, 012516 (2016).
  * Eiles  _et al._ (2016) M. T. Eiles, J. Pérez-Ríos, F. Robicheaux, and C. H. Greene, J. Phys. B: At. Mol. Opt. Phys. 49, 114005 (2016).
  * Fey  _et al._ (2019a) C. Fey, J. Yang, S. T. Rittenhouse, F. Munkes, M. Baluktsian, P. Schmelcher, H. R. Sadeghpour, and J. P. Shaffer, Phys. Rev. Lett. 122, 103001 (2019a).
  * Eiles  _et al._ (2020) M. T. Eiles, C. Fey, F. Hummel, and P. Schmelcher, J. Phys. B: At. Mol. Opt. Phys. 53, 084001 (2020).
  * Hunter  _et al._ (2020) A. L. Hunter, M. T. Eiles, A. Eisfeld, and J. M. Rost, Phys. Rev. X 10, 031046 (2020).
  * Lesanovsky  _et al._ (2006) I. Lesanovsky, P. Schmelcher, and H. R. Sadeghpour, J. Phys. B: At. Mol. Opt. Phys. 39, L69 (2006).
  * Kurz and Schmelcher (2014) M. Kurz and P. Schmelcher, J. Phys. B: At. Mol. Opt. Phys. 47, 165101 (2014).
  * Krupp  _et al._ (2014) A. T. Krupp, A. Gaj, J. B. Balewski, P. Ilzhöfer, S. Hofferberth, R. Löw, T. Pfau, M. Kurz, and P. Schmelcher, Phys. Rev. Lett. 112, 143008 (2014).
  * Gaj  _et al._ (2015) A. Gaj, A. T. Krupp, P. Ilzhöfer, R. Löw, S. Hofferberth, and T. Pfau, Phys. Rev. Lett. 115, 023001 (2015).
  * Niederprüm  _et al._ (2016) T. Niederprüm, O. Thomas, T. Eichert, C. Lippe, J. Pérez-Ríos, C. H. Greene, and H. Ott, Nat. Commun. 7, 12820 (2016).
  * Hummel  _et al._ (2019) F. Hummel, C. Fey, and P. Schmelcher, Phys. Rev. A 99, 023401 (2019).
  * Engel  _et al._ (2019) F. Engel, T. Dieterle, F. Hummel, C. Fey, P. Schmelcher, R. Löw, T. Pfau, and F. Meinert, Phys. Rev. Lett. 123, 073003 (2019).
  * Hummel  _et al._ (2021a) F. Hummel, K. Keiler, and P. Schmelcher, Phys. Rev. A 103, 022827 (2021a).
  * Hollerith  _et al._ (2019) S. Hollerith, J. Zeiher, J. Rui, A. Rubio-Abadal, V. Walther, T. Pohl, D. M. Stamper-Kurn, I. Bloch, and Christian Gross, Science 364, 664 (2019).
  * Duspayev and Raithel (2022) A. Duspayev and G. Raithel, Phys. Rev. A 105, 012810 (2022).
  * Fey  _et al._ (2019b) C. Fey, F. Hummel, and P. Schmelcher, Phys. Rev. A 99, 022506 (2019b).
  * Carter and Handy (1982) S. Carter and N. Handy, Mol. Phys. 47, 1445 (1982).
  * Handy (1987) N. Handy, Mol. Phys. 61, 207 (1987).
  * Meyer  _et al._ (1990) H.-D. Meyer, U. Manthe, and L. S. Cederbaum, Chem. Phys. Lett. 165, 73 (1990).
  * Beck  _et al._ (2000) M. H. Beck, A. Jäckle, G. A. Worth, and H.-D. Meyer, Phys. Rep. 324, 1 (2000).
  * Meyer and Worth (2003) H.-D. Meyer and G. A. Worth, Theor. Chem. Acc. 109, 251 (2003).
  * Meyer  _et al._ (2006) H.-D. Meyer, F. Le Quéré, C. Léonard, and F. Gatti, Chem. Phys. 329, 179 (2006).
  * Doriol  _et al._ (2008) L. J. Doriol, F. Gatti, C. Iung, and H.-D. Meyer, J. Chem. Phys. 129, 224109 (2008).
  * Meyer (2012) H.-D. Meyer, Wiley Interdiscip. Rev.: Comp. Mol. Sci. 2, 351 (2012).
  * (73) G. A. Worth, M. H. Beck, A. Jäckle, O. Vendrell, and H.-D. Meyer, The MCTDH Package, Version 8.2, (2000). H.-D. Meyer, Version 8.3 (2002), Version 8.4 (2007). O. Vendrell and H.-D. Meyer Version 8.5 (2013). Versions 8.5 and 8.6 contains the ML-MCTDH algorithm. Used version: 8.6.1 (May 2022). See http://mctdh.uni-hd.de/.
  * Harris  _et al._ (1965) D. O. Harris, G. G. Engerholm, and W. D. Gwinn, J. Chem. Phys. 43, 1515 (1965).
  * Hummel  _et al._ (2021b) F. Hummel, M. T. Eiles, and P. Schmelcher, Phys. Rev. Lett. 127, 023003 (2021b).