A comparative analysis of binding in ultralong-range Rydberg molecules

We perform a comparative analysis of different computational approaches employed to explore the electronic structure of ultralong-range Rydberg molecules. Employing the Fermi pseudopotential approach, where the interaction is approximated by an $s$-wave bare delta function potential, one encounters a non-convergent behavior in basis set diagonalization. Nevertheless, the energy shifts within the first order perturbation theory coincide with those obtained by an alternative approach relying on Green's function calculation with the quantum defect theory. A pseudopotential that yields exactly the results obtained with the quantum defect theory, i.e. beyond first order perturbation theory, is the regularized delta function potential. The origin of the discrepancies between the different approaches is analytically motivated.


Introduction
Diatomic ultralong-range Rydberg molecules consisting of a Rydberg atom whose electron, upon frequent scattering off a ground state atom, binds the atom, localized at large distances (∼ 10 3 Bohr radii), to the ion core of the Rydberg atom were predicted theoretically in Ref. [1]. A subset of these Rydberg molecules with s-wave dominated electronic orbitals were later observed [2]. This initial experimental work was followed by a flurry of new observations of Rydberg molecules with a variety of electronic and rovibrational structures [3][4][5][6][7][8][9][10]. For electronic s-wave scattering of the Rydberg electron from the ground state atom, two types of ultralong-range molecular states can be distinguished: the non-polar, low angular momentum quantum defect states and the polar, high angular momentum "trilobite" states. Theoretically, some aspects of these high angular momentum states have been explored, notably the semiclassical nature of the rich oscillatory structure in the adiabatic potential energy surfaces [11], their large permanent electric dipole moments [12], as well as the precise control over their electronic properties and molecular orientation by static electric and magnetic fields [13,14].
Further emphasis has been laid on polyatomic systems consisting of a Rydberg atom bound to two or more ground state atoms [15,16] or to a diatomic polar perturber [17,18]. Theoretical approaches to describe ultralong-range Rydberg molecules can be divided into two categories: methods using the Fermi pseudopotential [1,19] and methods using the quantum defect approach [20]-for a recent overview, we refer the reader to [21]. Although the two approaches are expected to be equivalent, they indeed differ in numerical implementation [19] and are even used alternatively [3].
In this work we therefore aim at a concise analysis of the origin of the discrepancies obtained from these two approaches. We thereby point out the interconnection and limitations in both approaches which should be taken into consideration in future studies of ultralong-range Rydberg molecules. Specifically, we show that a pseudopotential modeled by a bare delta function potential leads to non-converging molecular potential energy curves in basis set diagonalization of the electronic Hamiltonian. Modeling the pseudopotential instead by a regularized delta function reproduces potential energy curves which agree exactly with those obtained by the quantum defect theory approach. Throughout our analysis we not only give limits of the bare delta function potential but also its validity as an approximation in first order perturbation theory. This ultimately links the two approaches.

Ultralong-range Rydberg molecules
We consider a Rydberg atom whose ionic core is located at the origin and a neutral ground state atom located at the position R within the Rydberg electron orbit. In the Born-Oppenheimer (BO) approximation, the electronic Hamiltonian for this Rydberg molecule readsĤ whereĤ 0 is the Hamiltonian describing the Rydberg electron in its ionic core potential and V (R) is the interaction between the Rydberg electron and the neutral ground state atom which we call from now on the perturber. If V is short-ranged, the corresponding electronic wave function Ψ(r) outside the range of V can be determined from the scattering phase shifts induced by V . Although a more general treatment is possible, here we focus exclusively on pure s-wave electron-perturber scattering.
In the region where the electron is near the perturber, but outside of the range of V , the electron wave function behaves as where ρ = |r − R| is the relative electron-perturber coordinate and δ(k) is the s-wave phase shift depending on the wave vector k given semiclassically as The s-wave scattering length of the potential V is obtained from a[k] = − tan(δ(k))/k, which in the low energy limit behaves as where a[0] is the zero-energy scattering length and α[0] is the zero-frequency (static) polarizability of the ground state atom (polarized in the charge-atom interaction) [2,22]. The momentum k can be considered to be small (low-energy limit) because R in (3) is typically large and close to the classical turning points. The idea behind the Fermi pseudopotential approach is to replace V by a zero-range pseudopotential that possesses the same s-wave scattering length as the original potential, while the quantum defect approach replaces the potential V by the Dirichlet boundary condition (2).

Bare delta function potential
The pseudopotential usually employed to describe the s-wave binding of ultralong-range Rydberg molecules [1] is the bare delta function potential where the scattering length a[k] depends on the energy E and the perturber position R via (3) and (4). In the case that one is interested only in energies E close to some reference energy E 0 i one may approximate We now consider the energy shift of the electronic Rydberg orbit by the perturbing potential V where |ϕ i is an eigenstate of the bare Rydberg Hamiltonian with eigen For a compact notation we choose the basis states |ϕ i such that they diagonalize V inside energetically degenerate submanifolds, i.e.
Using (6), we deduce from (5) the first order perturbation theory correction to E 0 i as E which is the BO potential, valid for low angular momentum quantum defect, as well as for high angular momentum "trilobite" states.
As an example we consider the Rb(35s) + Rb(5s) state studied in [2,3]. The correspondent eigenstates ofĤ 0 are the quantum defect states |nlm that possess the where ∆ l denotes the l-dependent quantum defect. We use ∆ 0 = 3.13, ∆ 1 = 2.65, ∆ 2 = 1.35 and neglect the quantum defects of higher, i.e. ∆ l>2 = 0, angular momentum states. Here n, l and m are the usual hydrogenic quantum numbers. The properties of the perturber Rb(5s) enter in (4) via the dominating triplet s-wave scattering parameters where a 0 is the Bohr radius. The e − -Rb(5s) triplet scattering length has been calculated in [23] and the polarizability is obtained from [24].
eff are the first order approximations based on (8), where E (1) eff uses the effective zero-energy scattering length a eff [0] = −18.5a 0 . The curve E reg (R) gives the energy for the regularized delta function potential based on (20) while the three curves E j are determined via numerical diagonalization of (11) including basis sets 32 − j ≤ n * ≤ 32 + j.
The BO potential energy curve E (1) is shown (red dashed line) in Fig. 1. At R ≈ 1900a 0 , it possesses a minimum that lies 26 MHz below the dissociation limit (here set to zero) and supports several bound states. The non-differentiability around R ≈ 2050a 0 occurs due to the non-analytic behavior of k and E at the classical turning points. Additionally, the dashed orange curve Fig. 1 shows the lower-lying potential curve E (1) eff based on the effective scattering length a eff [0] = −18.5a 0 which was found out to agree better with the experimentally observed spectra [2]. As we will outline in Sec. 4, one reason for the discrepancy between a eff [0] and a[0], which has not been fully understood so far, is the effect of couplings between non-degenerate |ϕ i and |ϕ j basis states that, at the level of first order perturbation theory, are neglected in (8) but lead to deeper binding energies. This coupling becomes especially important in cases where the quantum defect of the state in question is close to an integer value.
To take the couplings into account, which means to analyze the system beyond the first order perturbation theory, we diagonalizeĤ numerically in an extended basis set. Technically this is done by fixing k as in (6) to some energy E 0 i and diagonalizingĤ in a basis set of states range-bound to a certain number of quantum defect eigenstates |nlm whose energies lie close to E 0 i . In the case of Rb(35s) + Rb(5s), we choose E 0 i = E 0 35,0,0 . Furthermore, it is convenient to group the basis states into manifolds of similar energy, that we call n * -manifolds. Each consists of the three quantum defect states |n * + 3, 0, 0 , |n * + 2, 1, 0 , |n * + 1, 2, 0 and the n * − 3 energetically degenerate states |n * , l > 2, 0 . For practical purposes and without loss of generality, we can neglect states with m = 0. The projection of the electronic angular momentum onto the internuclear R axis parallel to the z-axis is a good quantum number. States with m = 0 have vanishing density on the z-axis and hence show no interaction via the potential V .
The resulting BO potential energy curves are shown in Fig. 1: Using only the n * = 32 manifold gives the BO energy E 0 (cyan dots), including the n * -manifolds with |n * −32| ≤ 1 leads to the BO energy E 1 (light blue dots) and taking |n * −32| ≤ 5 results in the BO energy E 5 (dark blue dashed-dotted line). Here every increase in the size of the basis set leads to considerably lower potential curves. Surprisingly, even for large basis sets consisting of up to 10 n * -manifolds, the potential curves do not converge.
The following calculation shows that this non-convergence is not a numerical artifact but inherent to the bare delta function potential used for the electron-perturber interaction and the diagonalization procedure. We will illuminate this in the following. Let us consider the finite subspace W spanned by N eigenstates |ϕ i and investigate the diagonalization ofĤ inside W in more detail.
The projection of the Schrödinger equation where we made the ansatz Consequently, (12) will lead to the exact solution of (10) in the limit N → ∞. From (11) we deduce where N is a normalization constant. Reinserting this into (11) yields the characteristic equation The roots E of (14) are now exactly the eigenvalues that one would obtain by a numerical diagonalization of (11). However, in contrast to numerical diagonalization methods there is no need to perform the approximation (6) for the involved k-vector in numerical root finding algorithms. Hence, (14) is even more general.
To see the effects of diagonalization inside large basis sets, we take the limit N → ∞. The sum in (14) can then be replaced by the Green's function ofĤ 0 evaluated in the limit r = r ′ = R, where the sum above includes discrete and continuum states. However, for r → r ′ the Green's function diverges as Therefore the characteristic equation (14) possesses for N → ∞ no well defined solution E. This explains the nonconvergent behavior that we encounter with increasing size of the basis set for the diagonalization, as illustrated in Fig. 1.

Regularized delta function potential
In the literature it is well known that the usefulness of a delta function potential in three dimensions is restricted since its action on irregular wave functions is not well defined and its spectrum is unbounded from below [25,26]. We employ therefore the regularized delta function potential whose solutions fulfill the condition (2) exactly [27]. Using again the basis functions ϕ i we obtain Similar to (14) and as shown in [28] we derive that a nontrivial solution of (18) exists only for energies satisfying In contrast to (14) this equation is not restricted to finite regions of the Hilbert space and yet well defined. It is crucial to carry out the summation before taking the derivative in (18) and (19) as the opposite, i.e. shifting the derivative into the sum, would again lead to irregularities, due to (16). To this end, we insert explicitly the Green's function (15) into (19) which yields where This is exactly the result obtained alternatively via the quantum defect theory in [20,29] and agrees nicely with the expectation that the pseudopotential approach and the quantum defect approach should result in equivalent potential curves. Eq. (20) can be evaluated numerically by expressing the Green's function in (15) in terms of Whittaker functions, which is possible for arbitrary quantum defects [30]. The BO potential energy curve E reg (R) obtained in this manner for the Rb(35s)+Rb(5s) state with the parameters (9) is shown (orange dashed curve) in Fig 1. Compared to the first order perturbation theory results, the outer minimum is around 30% deeper and shifted towards smaller R. This is an indication of the contribution of higher l-states and explains why the observed binding energies in [2] were larger than expected from the first order perturbation approximation (8). In [2], this was compensated partially by the introduction of an effective scattering length a eff [0]. More analytical insight is gained by performing first order perturbation theory with (17). Additional care is necessary when using the regularized delta function potential in standard perturbation theory, due to the non-commutativity of summation and differential operations [27]. Nevertheless, the first order corrections E (1) i to the energy E 0 i can be obtained in the normal manner as where energetically degenerate ϕ i are chosen such that they diagonalize V in the degenerate manifold. Hence the first order energies of the regularized delta function potential (22) coincide with the first order energies of the bare delta function potential (8). Therefore, even though we showed that the bare delta function potential is in general not adequate, it serves as a valid first order approximation for the regularized delta function potential.
More generally, Fig. 1 suggests that the bare delta function potential may even approximate the regularized delta function potential beyond the first order perturbation theory, by including only a few more, but not too many, adjacent states. For example the BO potential energy curve E 0 obtained by including only states adjacent to the n = 32 manifold, lies close to E reg (R). Here a more rigorous comparison between the gain in accuracy and the error accumulated when increasing the size of the basis set, would be desirable and should be carried out in the future. This would be especially useful as the regularized delta function potential can not be implemented in a numerical diagonalization procedure that relies on the basis functions ϕ i which are regular at R. We see this from (18), which, in contrast to (11), can not be converted into a system of linear equations for the coefficients α i because the sum has to be carried out before taking the differential.

Summary & Conclusion
Aspects of the numerical implementation of the Fermi pseudopotential approach, relevant to the calculation of the BO potential energy curves of the ultralong-range Rydberg molecules bound by pure s-wave scattering are delineated. It is shown that the bare delta pseudopotential leads to unphysical results that diverge with increasing basis set size in any diagonalization scheme. This behavior is described analytically. The convergent BO potential curves, that agree exactly with the results obtained within the quantum defect theory can be produced by employing a regularized delta function potential. Although the bare delta function potential still yields the correct binding energies in the first order perturbation theory, the example of the Rb(35s) + Rb(5s) molecule indicates that beyond the first order effects are not negligible. By including some of the neighboring Rydberg basis states, some of these corrections can be reproduced. However it remains to be seen how many additional states can rigorously be included in exact diagonalization schemes and whether the zero-range pseudopotential for the p-wave-interaction [19], will require similar regularization as well.