Vibronic interactions in trilobite and butterfly Rydberg molecules

Ultralong-range Rydberg molecules provide an exciting testbed for molecular physics at exaggerated scales. In the so-called trilobite and butterfly Rydberg molecules, the Born-Oppenheimer approximation can fail due to strong non-adiabatic couplings arising from the combination of radial oscillations and rapid energy variations in the adiabatic potential energy curves. We utilize an accurate coupled-channel treatment of the vibronic system to observe the breakdown of Born-Oppenheimer physics, such as non-adiabatic trapping and decay of molecular states found near pronounced avoided crossings in the adiabatic potential curves. Even for vibrational states localized far away from avoided crossings, a single channel model is quantitatively sufficient only after including the diagonal non-adiabatic corrections to the Born-Oppenheimer potentials. Our results indicate the importance of including non-adiabatic physics in the description of ultralong-range Rydberg molecules and in the interpretation of measured vibronic spectra.


I. INTRODUCTION
Molecular states of Rydberg atoms and atoms or ions feature exotic binding mechanisms and offer novel opportunities to investigate fundamental physics. In ultralongrange Rydberg molecules, the Rydberg atom binds to ground-state atoms via an attractive interaction induced by the scattering of the Rydberg electron from the atoms [1][2][3]. In charged Rydberg molecules, the Rydberg atom binds to an ion due to long-range electrostatic interactions leading to an avoided crossing of an attractive polarization potential and a repulsive dipole potential [4][5][6]. In Rydberg macrodimers, two Rydberg atoms bind at the avoided crossing of attractive and repulsive vander-Wals potentials [7][8][9]. All of these molecules feature bond lengths on the 0.1-10 µm scale leading to slow vibrational dynamics which can be observed in real time [10].
Due to the high density of states in Rydberg atoms and the correspondingly small energy gaps, non-adiabatic couplings play a crucial role in Rydberg molecules. They are known to modify the vibrational spectrum of macrodimers [11] and provide a decay channel for charged Rydberg molecules [12]. In ultralong-range Rydberg molecules, non-adiabatic couplings have been considered at the level of the Landau-Zener approximation to predict electronic transitions [13]. Conical intersections in their molecular potentials can be exploited to control atomic collisions and lifetimes over a wide range of principle quantum numbers and temperatures [14].
Ultralong-range Rydberg molecules have very unusual properties, such as oscillatory Born-Oppenheimer po-tentials, rapid changes of electronic character, and extreme sensitivity to external fields. Therefore, they provide a unique environment to test molecular physics on slow vibrational time scales, allowing direct observation of effects arising from the entanglement of vibrational and electronic degrees of freedom and offerering abundant opportunities to examine concepts of molecular physics in energetic and spatial regimes very different from those relevant for ground state molecules. Despite these promising opportunities, an ab initio treatment of non-adiabatic couplings in ultralong-range Rydberg molecules has not been performed yet.
In this article, we study ultralong-range Rydberg molecules in a complete coupled-channel treatment for the electronic and vibrational degrees of freedom. The arising non-adiabatic couplings lead to vibronic interactions. We focus on molecular states formed in the vicinity of the avoided crossing of the trilobite and butterfly electronic states. These polar molecules have high electronic orbital angular momentum character and their non-adiabatic coupling can become singular. We discuss several strategies for diabatization of the adiabatic representation of the total molecular wave function in order to solve the molecular Hamiltonian.
We find that non-adiabatic corrections to the molecular potentials are crucial for an accurate interpretation of molecular spectra. This sheds new light on measurements of electron scattering phase shift relying on Rydberg molecules [15][16][17][18][19], as the level shifts caused by nonadiabatic corrections could be mis-identified as stemming from inaccurate scattering properties. Furthermore, offdiagonal non-adiabatic couplings lead to the formation of molecular states without single-channel counterparts. Their structure resembles the quantum reflection states observed in single channel calculations [20,21], although they are classically trapped by a potential barrier created by non-adiabatic coupling. Other classes of molecular states decay non-adiabatically. Spectral signatures for the experimental observation of the vibronic interaction effects are discussed.
This manuscript is structured as follows. In section II, we lay out the theory and methods used throughout. Section II A introduces in detail the non-adiabatic representation before we elaborate on diabatization in II B. This is to promote and simplify adaptation of this framework for future investigations. Details of the model for ultralong-range molecules and the physical parameters are provided in II C. Section III contains the results, firstly introducing the adiabatic potential energy curves in III A, secondly considering variations in the electronic structure in III B, thirdly providing the molecular spectra in III C, and finally focusing on vibronic interactions in III D. We provide a conclusion and an outlook in section IV.

II. THEORY & METHODS
We consider two atoms, one excited to a high-lying Rydberg state, separated by a large distance R along the z axis of the coordinate system. The total Hamiltonian reads (in atomic units) where µ is the reduced mass of the diatomic system, r is the coordinate of the Rydberg electron, V c is the interaction between the Rydberg electron and the ionic core, and V ea the interaction of the Rydberg electron with the distant ground-state atom. The last three terms of equation (1) form the electronic HamiltonianĤ e (r, R) consisting of the corresponding kinetic energy and two potential terms, one for each heavy particle. To solveĤ, we employ the adiabatic representation and expand the total wave function in terms of solutions of the electronic Hamiltonian [22,23] such that the sum over R-dependent coefficients α |χ α | 2 can be interpreted as a vibrational wave function of the nuclei, and defines the adiabatic electronic eigenstates ψ α with eigenvalue V α given a fixed nuclear configuration R. V α (R) are the adiabatic potential energy curves.

A. Non-adiabatic couplings
Applying the Hamiltonian (1) to equation (2) and projecting out the electronic eigenstates leads to a general set of coupled-channel, radial Schrödinger equations for the vibrational motion [24] for a total energy E. ∇ R indicates the derivative with respect to the radial, nuclear coordinate R. The righthand side contains the non-adiabatic coupling operator which consists of the derivative coupling elements and the non-adiabatic scalar coupling elements For all practical purposes, the sum in equations (2) and (4) needs to be truncated. In the Born-Oppenheimer approximation, all non-adiabatic terms are neglected. The right-hand side of equation (4) vanishes and the resulting decoupled equations can be solved individually. In the Born-Huang approximation, only the diagonal non-adiabatic terms are considered. In many cases, these terms are crucial to capture qualitative corrections [25,26]. It is helpful to absorb the diagonal couplings in the potential energy curves through the definition Note that P αα = 0, due to anti-hermiticity of the derivative coupling operator. Since the diagonal elements Q αα are strictly negative (compare equation (12) below), this leads to a shift to larger potential energies W α ≥ V α . In principle, the non-adiabatic coupling elements are straightforwardly obtained from the adiabatic electronic states. However, numerical difficulties frequently arise in practice. With respect to the electronic states, the operator Λ αα is hermitian, while P αα is anti-hermitian. Hence, Q αα is non-hermitian and numerical inaccuracies or round-off errors in its computation can easily lead to an overall non-hermitian matrix and non-physical complex energy eigenvalues. This can be avoided by recasting and using partial integration to ensure hermiticity of the matrix elements of the non-adiabatic coupling. To this end, let us introduce a basis of the vibrational Hilbert space B i (R), such that any R-dependent function can be conveniently expanded as with constant coefficients c i . The matrix elements of the non-adiabatic coupling operator in this basis are As the first term is the product of two anti-hermitian operators in the vibrational and electronic Hilbert subspaces, respectively, and the second term is the product of two hermitian operators, this guarantees a total hermitian representation. It additionally avoids the need to compute second derivatives of the electronic wave functions, because In conclusion, we can solve equation (4), e.g. by performing an exact diagonalization in the electronic and the vibrational degree of freedom. In some circumstances, however, the derivative coupling elements may contain highly singular points such that equation (4) becomes formally ill-defined [27]. It is then reasonable to switch to a representation that circumvents such difficulties.

B. Diabatic representation
We can recast equation (4) in the following form where the first term is understood to contain a dressed kinetic energy operator that has absorbed all non-adiabatic couplings. Note that equation (13) is a matrix equation whereV is a matrix that contains the V α along the diagonal and χ is a vector of all χ α . Equation (13) is invariant under unitary transforma-tionsŜ(R). Exact diabatization is achieved whenŜ satisfies [28] (∇ RŜ +PŜ) = 0, in which case equation (13) becomes whereχ =Ŝ † χ is the total wave function expressed in the diabatic representation andÛ =Ŝ †VŜ is the diabatic potential matrix. In general,Ŝ is not unique, but as long as there is only a single adiabatic parameter (in our case R), equation (14) does possess a unique solution and can be solved analytically [29]. Considering two adiabatic electronic states,Ŝ is given as a simple rotation about the angle The solution for three and four electronic states is given in Appendix D. There exist other strategies for diabatization [30][31][32], typically related to achieving an electronic representation that provides continuity of some physical observable, such as electronic dipole moments or charge distributions. These are especially useful ifŜ is not known. We discuss dipole diabatization and a related approach in Appendices B and C.

C. Ultralong-range Rydberg molecules
We consider ultralong-range Rydberg molecules of rubidium. The adiabatic potential energy curves are obtained from diagonalizing the electronic Hamiltonian in a basis of atomic states φ nlm . These states are labeled by the principle and orbital angular momentum quantum numbers of the Rydberg atom n, l, and m l , and satisfy the Rydberg atom's Schrödinger equation where E nl = − 1 2(n−µ l ) 2 and µ l is an angular-momentumdependent, but energy-independent quantum defect. These quantum defects detune the atomic states with low angular momentum l < 3 from the degenerate Rydberg manifold − 1 2n 2 . We use quantum defect data from [33].
The electron-atom interaction is modeled using the Fermi pseudopotential [34,35] where a s(p) is the scattering length (volume) in the s-(p)wave channel. These depend on R via the semiclassical kinetic energy of the electron, k 2 = 2E nl + 2/R. We use scattering data from [19], where the zero-energy s-wave scattering length is a s (0) = −15.24 a 0 . To represent the electronic degrees of freedom for a Rydberg molecule formed by a Rydberg atom excited to the principal quantum number n, we include in our atomic basis all states with n and n − 1. From this, we obtain good qualitative agreement with more accurate Kirchhoff integral methods based on the Coulomb Green's function [36,37]. This descriptive level serves our goal to quantify and elucidate the importance of non-adiabatic physics, which can be obscured by the complexity of the full molecular picture. This would include quantitative corrections to the electronic structure resulting from higher partial wave scattering [38] and the electronic and nuclear spins of both atoms [39,40]. To represent the vibrational degree of freedom, we employ B splines [41]. Due to the banded structure of the involved matrices, this is numerically very efficient. Depending on the principle quantum number, we use 500 to 1500 splines of order 12 and a constant spacing of knot points to converge the vibrational energies and wave functions.
To simulate an absorption spectrum, the standard experimental observable to study Rydberg molecules, we account for the dipole transition strength of the total wave function assuming a two-photon excitation from the 5s ground state via the 6p intermediate state and for the Franck-Condon factor assuming a homogeneous and isotropic gas of ground-state atoms, i.e. a uniform distribution of the initial vibrational state: where the wave functions are defined in equation (2) and below. Specifically, only the atomic S-and D-state character of the total wave function contribute to the spectral signal. Naturally, a one-photon transition directly from the 5s ground state changes the spectral characteristics. Furthermore, a Gaussian profile for the line broadening of the excitation laser with a width of 1 MHz is assumed.

III. RESULTS & DISCUSSION
In the following, we present the characteristic properties of the electronic, the vibrational, and the spectral structure of neutral Rydberg molecules taking into account non-adiabatic couplings. To illustrate the main physical phenomena, we focus on a representative subset of principle quantum numbers n = 43 to 45, through which we can elaborate on generic n-dependent features that we found studying a larger range between n = 25 and 70.
A. Structure of the adiabatic potential curves Let us first consider the adiabatic potential curves shown in figure 1 for n = 43. The level structure of Rydberg molecules is such that these fall into two classes. The electronic states of dominant low-l character are, due to their quantum defects, only weakly perturbed by the interaction with the ground state atom. Their potential curves are, on the scale of figure 1, nearly flat. Two additional adiabatic potential energy curves split away from the degenerate Rydberg manifold. Loosely speaking, each of these can be associated with the interaction channels given by the two terms in equation (18). The resulting s-wave dominated and p-wave dominated molecular states are called trilobite and butterfly states, respectively [1,36]. These two states have a prominent avoided crossing [37] shown in the inset. Non-adiabatic effects are expected to be strongest in its vicinity.
Based on this Rydberg level structure, the adiabatic eigenstates ψ 1 and ψ 2 defining the potentials V 1 and V 2 shown in figure 1 are subject to several changes in the dominant electronic character. Three avoided crossings occur for V 1 (solid blue). Starting asymptotically at large R, V 1 has dominant F -state character. Going inward towards smaller R-values, firstly V 1 attains trilobite char- acter, secondly attains butterfly character, and lastly, at much lower energies, attains S-state character. For V 2 on the other hand (solid red) four avoided crossings occur. Starting asymptotically at large R, V 2 has dominant trilobite character. Going inward towards smaller R-values, firstly V 2 attains F -state character, secondly attains butterfly character, thirdly attains trilobite character, and lastly attains F -state character once more. We obtain converged results for the vibronic spectra in the examined energy ranges shown in the following sections by selecting only these two adiabatic states for the calculation. No deviations were found between these calculations and those performed with three and four channels, including also the asymptotic S and F adiabatic potentials shown in figure 1. Additionally, to obtain molecular states, we set a variable hard wall boundary at small R 0 and average the results over 50 different R 0 values. This is the stabilization method established for ultralong-range molecules in reference [20].

B. Variation of the electronic structure
Let us consider how the adiabatic potential energy curves, the off-diagonal derivative coupling element, and the diabatic potential curves, vary as a function of n. These quantities are plotted in figures 2. The derivative coupling element P αα is a prefactor to a derivative operator of the vibrational degree of freedom in the total Schrödinger equation (4). Therefore, multiplying P αα by a unit length and dividing by the molecule's reduced mass results in a unit of energy, allowing us to show it alongside the electronic adiabatic potential energy curves. Hence, we can compare the relative strengths of non-adiabatic  couplings for molecules at different n.
All of the avoided crossings discussed previously are accompanied by signatures in the derivative coupling element P 12 , shown as solid green line. Note for example the peak at the outermost avoided crossing at large R. The most unusual feature of P 12 is its relatively large value over the range where the two adiabatic states have trilobite and butterfly character at the same R values. This is also the region where P 12 has the most significant feature, namely a one or two (depending on the principle quantum number) sharp peaks around the avoided crossing. At n = 43, in the left-most subfigure, we find a symmetric double peak structure in P 12 . At n = 44, in the middle subfigure, we find a very sharp peak in P 12 that extends beyond the range of the figure. At at n = 45, in the right-most subfigure, we find an asymmetric double peak structure in P 12 . These are the main characteristics of the derivative coupling element for the two underlying adiabatic states over the full range of n between 25 and 70. The patterns repeat in increments of approximately 3 to 4 in n. In the case of a single sharp peak, the avoided crossing of the two potentials becomes extremely narrow. This is remarkable given the overall energies of the molecular system. It resembles the cases of conical intersections discussed in reference [14]. Due to the singular character of the coupling elements, methods to solve the total Hamiltonian (1) that rely on the adiabatic representation fail. It is useful to transform to the diabatic representation.
The elements of the diabatic potential matrix are shown as dashed lines in figure 2. The diagonal elements (diabatic potential curves) are shown in purple and blue and the off-diagonal coupling is depicted in beige. Notably, over all regions where the dominant character of the adiabatic states are trilobite and butterfly, all elements of the diabatic potential matrix have strong oscillatory character. As the rotation angle of the trans-formation between the adiabatic and diabatic representation is the integral of the derivative coupling element P 12 , this behavior may not seem surprising. After all, P 12 is large in that entire region. It is, however, rather unusual, as the diabatic potentials do not provide any physical intuition for the molecular system, as would be desirable for diabatic potentials. The oscillatory behavior is known, however, to occur for states which are coupled by a large dipole transition element [42]. This is certainly true for the trilobite and butterfly state. We notice that the n-dependent structure of the derivative coupling element does not have a significant influence on the structure of the diabatic potential matrix elements. A simplified model including only trilobite and butterfly and neglecting quantum defect states is studied in Appendix A reproducing these diabatic potentials.
The overall magnitude of the non-adiabatic coupling elements, both for P αα and Q αα , is approximately constant across all n. On the other hand, the energies of the adiabatic potentials scale as n −3 , in line with the scaling of the level splittings in Rydberg atoms. As a result, non-adiabatic effects play a more important role at large n. The diagonal scalar coupling terms lead to an overall energy shift of the potentials on the order of 5 MHz. At low n 30, this is more or less negligible. However, at large n 60, this shift is larger than the energy spacing of molecular states. Figure 3 shows the adiabatic potentials for n = 65 as solid lines along with the Born-Huang corrected potentials, which include the diagonal scalar coupling, as dashed lines. The lower inset shows that the Born-Huang terms give rise to a correction on the order of 5 MHz, which leads to a significant overall shift of the molecular states. Additionally, the upper inset shows the potential barrier that appears at the narrow avoided crossing due to these diagonal nonadiabatic coupling terms. This barrier is a result of the narrow energy gap of the two states at the crossing and  the correspondingly sharply peaked derivative coupling element.

C. Vibrational structure and molecular spectrum
Let us now turn to the energy eigenfunctions of the total Hamiltonian (1). Figure 4 shows the simulated spectral signal obtained from the total wave functions depending on the underlying model: In blue, the single channel model using the adiabatic potential energy curves V 1 and V 2 , the Born-Oppenheimer spectrum. In orange, the single channel model using the Born-Huang corrected potential energy curves W 1 and W 2 , which include the diagonal scalar couplings Q 11 and Q 22 , the Born-Huang spectrum. And in green (mirrored for visibility), the two-channel model in the diabatic representation, the exact spectrum.
Two features are immediately striking. Firstly, both the Born-Huang spectrum and the exact spectrum are shifted in energy compared to the Born-Oppenheimer spectrum. Secondly, the Born-Huang spectrum reproduces the exact spectrum quite accurately over most of the energy range. This clearly shows the importance of including at least the diagonal non-adiabatic couplings for investigations of ultralong-range Rydberg molecules. That said, we do find signals in the exact spectrum that have no counterpart in the single channel spectra. At n = 43, in the left-most subfigure, around E = −3.68 GHz there is a broad peak next to a sharp peak with rather strong signal. The single channel spectra only show a single asymmetric peak of weaker signal. Around E = −3.51 GHz, there is a weak peak in the exact signal not found in the single channel spectra. At n = 44, in the middle subfigure, the Born-Huang spectrum is astonishingly accurate. Some of the peaks with weaker signal are slightly shifted towards higher energies compared to the exact spectrum. At n = 45, in the right-most subfigure, around E = −3.18 GHz the exact spectrum features a double peak structure significantly shifted from a similar but weaker in signal double peak in the Born-Huang spectrum. Interestingly, the Born-Oppenheimer spectrum predicts a single peak at the correct energy. Between E = −3.03 and E = −3.07 GHz, the single channel spectra each predict a peak completely absent in the exact spectrum.

D. Off-diagonal non-adiabatic couplings
Some of these spectral features have their origin in strong off-diagonal non-adiabatic couplings that cannot be described in a single channel model. In order to gain a better understanding of these situations, we now turn to study the molecular wave functions directly. Figure 5 shows selected molecular eigenstates after tracing over the electronic degree of freedom, shifted according to their eigenenergies. The corresponding squared wave functions can therefore be interpreted as vibrational probability distributions. Here, states that are strongly suppressed in the spectra due to their vanishing Franck-Condon overlaps are not shown. This holds also for states resembling box states, which are artifacts related to the inner hard-wall boundary. The energy window is the same as that for the spectra plotted in figure 4. The gray curves represent the adiabatic potentials V 1 and V 2 shown in figure 2. Additionally, the black curves represent the Born-Huang corrected potentials W 1 and W 2 . In line with the color scheme used in figure 4, the eigenstates are obtained from the different models: Born-Oppenheimer states in blue, Born-Huang states in orange, and exact states in green. As the spectra suggest, most of the molecular wave functions obtained from the different models do not differ strongly. For example, the states that localize in the individual potential wells of the trilobite potential only deviate in energy. These are the states that correspond to strong signals in figure 4. Molecular states that have no counterpart in a different model localize in the vicinity of the avoided crossing where non-adiabatic couplings are particularly strong; these are located in the left-most region of each subfigure.
For n = 43, in the left-most subfigure, most notably around E = −3.68 GHz, a molecular state localizes at around R = 1650 a 0 . Although this state appears to be trapped in the Born-Huang potential well, the lack of a corresponding Born-Huang state shows that the potential well alone is not suitable to support such a state. In other words, this molecular state is trapped non-adiabatically by the off-diagonal derivative couplings. The exact state has a very small but finite amplitude towards small R, resembling states trapped by internal quantum reflection, discussed extensively throughout the literature on ultralong-range Rydberg molecules [19][20][21]. At higher   energy, around E = −3.55 GHz, and at similar internuclear distance, there is a Born-Oppenheimer state localizing in a potential well that is shifted upwards by the Born-Huang correction. This state cannot be identified in the spectrum because it is degenerate in energy with molecular states localizing at large internuclear distances, around R = 2500 a 0 , where non-adiabatic couplings are generally much weaker. The corresponding Born-Huang state is at approximately E = −3.53 GHz, barely visible in the figure. There is no exact counterpart to these state. In other words, these states are subject to rapid non-adiabatic decay into the continuum of the lower potential. At n = 44, in the middle subfigure, the only region with significant differences between different models is again the potential formed by the avoided crossing, where single channel states localize, but exact states do not occur. At n = 45, in the right-most subfigure, the overall structure is very similar to n = 43.
To summarize the most important observations here, we find that the potential well formed by the avoided crossing of the trilobite and butterfly does not support localized molecular bound states as would be suggested by single channel calculations. This is in marked contrast to other types of Rydberg molecules, most prominently Rydberg macrodimers, where a potential well is formed by the crossing of induced dipole potentials, and in charged Rydberg molecules, where a potential well is formed by the crossing of an attractive polarization potential and a repulsive dipole potential. Both of these cases are known to support bound molecular states [6,11]. In ultralong-range Rydberg molecules, non-adiabatic couplings are too strong for such states to exist. On the other hand, non-adiabatic couplings can lead to the formation of quantum reflection type states at the same radial po-sition, but at energies much lower than the potential well. Counterintuitively, both the observed effects, nonadiabatic decay and non-adiabatic trapping are more pronounced at principle quantum numbers where the derivative coupling element features a double peak structure. In the case of a single narrow peak of P 12 , such effects are less significant. We attribute this to the fact that nonadiabatic couplings are strongly pronounced over a much wider range of R in double peak situations. Opposed to this, the single peak can become arbitrarily narrow. Its integral converges at π 2 corresponding to a complete inversion of the adiabatic character.
Overall, the vibronic effects described in this section are not limited to the quantum numbers n = 43 to 45. On the contrary, we find similar patterns over a wide range of n = 25 to 70. Independent of n, the total energy shift between exact molecular states and states obtained within the Born-Oppenheimer approximation is approximately 5 MHz. For states localizing in the vicinity of the avoided crossing, the shift can be much larger. Vibronic interactions effects are less pronounced at principle quantum numbers featuring highly singular P -matrix elements resembling conical intersections. These special situations occur repeatedly throughout the range of n between 25 and 70.

IV. CONCLUSIONS & OUTLOOK
We have studied vibronic interaction effects in ultralong-range Rydberg molecules focusing on trilobite and butterfly states. The avoided crossing between these two states gives rise to non-adiabatic trapping of molecular states including internal quantum reflection resonances. Non-adiabatic decay is also observed. Overall, non-adiabatic corrections to the adiabatic molecular structure lead to energy shifts of spectral lines on the MHz scale. Such energy shifts may have been falsely attributed to corrections of electron phase shifts in the past. Generally, in ultralong-range Rydberg molecules, a larger scattering lengths corresponds to deeper wells in the adiabatic potential energy curves. This is why global energy shifts caused by diagonal non-adiabatic couplings can easily be mistaken for for a smaller scattering length or volume.
This work opens a pathway to study non-adiabatic interaction effects in Rydberg molecules. We consider multiple scenarios to be of high relevance: Firstly, nonadiabatic decay of molecular states correlated to Rydberg quantum defect states. This is important to understand molecular lifetimes. Secondly, non-adiabatic effects are naturally more pronounced in lighter molecules such as Li*Li, Li*Rb, or Na*Na, where they might lead to more striking effects. Thirdly, the methods developed for this study can be applied to other species of Rydberg molecules such as charged molecules and macrodimers. Lastly, exposing molecules to additional external fields is a promising perspective, especially as the gap between trilobite and butterfly curves can be controlled by an external electric field [43]. It may lead to more control of molecular properties, while causing more complexity, such as many-state interaction dynamics.
Notably, in the scenario considered here, the molecular channels are coupled but closed. In ultralong-range Rydberg molecules, however, scenarios can be investigated involving open and closed channels. This introduces coupling to continuum states and can lead for example to Feshbach resonances. In general, Rydberg molecules are highly asymmetric and imbalanced objects. Similar systems are encountered in semi-conductor materials based on Rydberg excitons. Their non-adiabatic interaction effects may lead to unexpected material response properties. The goal is to find a representation of the electronic basis that is dipole decoupled. DiagonalizingD(R) equips us with a transformationS(R) between the adiabatic and the dipole diabatic representation. It is defined viâ whered(R) is a diagonal matrix. The Schrödinger equation becomes whereχ =S † χ,Ū =S †VS , andΛ =S †ΛS . While this is formally equivalent to equations (4) and (15), the point of this representation is that residual couplingsΛ are expected to be small, such that the right hand side of this equation is set to zero. In the following, we will show that this is not an accurate assumption in our setup.

Appendix C: Hybridized basis
The dipole diabatic representation is closely related to another approach originating in results found in references [44,45] and developed previously in references [46][47][48]. Here, a hybridized electronic basis employing analytic functions of the trilobite and butterfly states with two (R-dependent) normalization constants T 2 = β |φ β | 2 and B 2 = β |φ β | 2 is used to expand the total wavefunction. This ansatz suffers from the fact that the two basis states are not orthogonal. In fact, for large R, they become linearly dependent. The offdiagonal overlap matrix element here is 1|2 = M 2 T B , where M 2 = φ (R)φ(R). Still, the resulting potentialand non-adiabatic coupling terms provide some intuitive understanding of the physical situation. Employing the hybridized basis, the total Schrödinger equation contains a non-diagonal potential matrix where κ = {1, 2} labels the two basis states, and O −1 κκ is the matrix element of the inverse overlap matrix. The potential curves inυ almost exactly reproduce the dipole diabatic potentials for the trilobite and butterfly states inŪ . We attribute the slight differences merely to the fact that dipole diabatization provides an orthogonal basis. The elements ofυ are shown in figure 7 on the left. Although this approach reduces the singularity of non-adiabatic couplings, they remain nonzero in this approach. The derivative coupling matrix is not anti-hermitian, due the non-orthogonality of the underlying basis states. The additional terms encountered in this expression areM 2 = φ (R)φ(R) and N 2 = φ (R)φ (R). For the same reason, the scalar coupling matrixρ does not obey the usual relation (9). In fact, given this non-orthogonal basis, it is not clear how to analytically obtain a transformation to the adiabatic representation. After all, this transformation is not unitary at large R. While the analytic matrix elements of the potential and couplings are not straight forward, they all result in smooth functions of R and n, which are shown in figure 7. Notably, the derivative coupling elements (middle) become large asymptotically at large R.

Appendix D: Adiabatic to diabatic transformations
The unitary transformation between adiabatic and diabatic representation can be found analytically, as long as there is only a single adiabatic parameter (in our case R). As these transformations are often useful, however, to our knowledge typically not provided in the literature, we reproduce here the explicit transformations for three and four channel systems. Considering three adiabatic states, the unitary transformation is a general rotation in three dimensions that can be expressed via three angles, e.g. the Euler angles. Here, we choosê S 3D (R) =R Z (α) ·R X (π/2 − β) ·R Z (γ), whereR i (ϕ) is a rotation in three dimensions around the axis i by the angle ϕ. Solving equation (14) for a general three-state P matrix leads to a set of coupled differential equations for the three angles α, β, and γ: α = p 12 − tan β (p 13 cos α + p 23 sin α), (D3) β = p 13 sin α − p 23 cos α, (D4) γ = sec β (p 13 cos α + p 23 sin α). (D5) This approach extends to higher dimensional systems. A general rotation in d dimensions can be characterized by d(d − 1)/2 angles governing simple rotations that leave a (d − 2)-dimensional subspace unchanged. In four dimensions, general rotations are parameterized by six angles. The derivation and solution of the corresponding differential equations gets increasingly involved for more states. Let's call the four axes in four dimensions {w, x, y, z}, then simple rotations leave a plane unchanged. SayR ij (θ) is a rotation of angle θ in the i-j-plane, e.g.R and like-wise for δ such that the following set of differential equations is obtained for the six angles α = p 12 + a − 3 tan β + a − 4 sec β tan γ, (D17) δ = a − 3 sec β + a − 4 tan β tan γ − d − sec γ tan , (D20) = −d + sec γ,