Quantum reactive scattering in the long-range ion-dipole potential

An ion and a polar molecule interact by an anisotropic ion-dipole potential scaling as $- \alpha \cos (\theta)/r^2$ at large distances. Due to its long-range character, it modifies the properties of angular wave functions, which are no longer given by spherical harmonics. In addition, an effective centrifugal potential in the radial equation can become attractive for low angular momenta. In this paper, we develop a general framework for ion-dipole reactive scattering, focusing on the regime of large $\alpha$. We introduce modified spherical harmonics as solutions of angular part of the Schr\"odinger equation and derive several useful approximations in the limit of large $\alpha$. We present the formula for the scattering amplitude expressed in terms of the modified spherical harmonics and we derive expressions for elastic and reactive collision rates. The solutions of the radial equation are given by Bessel functions, and we analyse their behaviour in two distinct regimes corresponding, basically, to attractive and repulsive long-range centrifugal potentials. Finally, we study reactive collisions in the universal regime, where the short-range probability of loss or reaction is equal to unity.

Another powerful platform for fundamental research in quantum physics are ultracold gases of molecules [39,40]. Trapping of ultracold polar molecules in optical lattices leads to a variety of novel quantum phases or can be applied to perform quantum computations [39,40]. So far, quantum degenerate regime has been achieved only for bialkali dimers [41][42][43][44][45][46][47][48]. Bialkali molecules in the ro-vibrational ground state can be classified into reactive and non-reactive ones [49]. While reactive collisions can be explained by relatively simple quantum scattering models based on the properties of the long range potential [50,51], collisions of nonreactive molecules are far more complicated, as the scattering is affected by the presence of a dense spectrum of overlapping resonances, leading to, so called, sticky collisions [52]. Their theoretical treatment is based on methods derived from randommatrix theory [53,54]. Ultracold chemical reactions of molecules can be controlled by external fields [55], internal spin states [56], or by aligning them in optical lattice structures of reduced dimensions [57][58][59][60][61][62]. So far they have been studied experimentally in KRb [55,56,63], NaLi [48] and triplet Rb 2 [64]. Modern techniques in manipulation of single atoms in optical tweezers, have allowed to assembly ultracold molecules directly from two atoms in a single, controlled chemical reaction [65].
Recently, first steps have been done towards combining of cold polar molecules with cold molecular ions in a single experimental setup [66]. Motivated by these attempts, in this work we study quantum scattering of an ion with a polar molecule in the low-energy regime. Here, we consider only collisions in the long-range part of the interaction given by the ion-dipole potential and assume fixed orientation of the electric dipole moment in the course of the collision. We note, that in general spatial orientation of a molecule varies in time, and even ultracold polar molecules in the ground state of rotational motion rotate having no net electric dipole moment. In this sense, our study is a necessary prerequisite before performing more elaborated analysis including rotational degrees of freedom, and the effects of molecule polarization by the ion's charge or an external electric field. Hence, the full description of the scattering problem would require solving a set of close-coupled equations expanded in the basis of rotational states of a molecule. In this context, our solutions derived in the current paper can be useful as expansion basis of the relative motion, describing long-range behaviour of the wave function components.
Another situation when our treatment is directly applicable is the collision of a very light charged particle like electron or positron with a heavy molecule. In such collisions the electron (positron) energy is much higher than rotational constant, and the scattering calculations can be done for a fixed orientation of a polar molecule. For such systems, ion-dipole collisions have been systematically studied, in particular in the low-energy regime [67][68][69][70][71][72][73]. In this context, it is known that long-range iondipole potential −α cos(θ)/r 2 modifies properties of angular momentum wave function introducing correction to the centrifugal potential in the radial equation [74]. It was shown that for dipole moments larger than the critical value α cr = 1.279, the potential becomes too at-tractive (at least in some directions) [75], and the collapse to the center takes place [68].
In collisions of atomic or molecular ions and polar molecules, typical values of α parameter are very large, and such a regime requires a separate analysis. So far, ion-molecule collisions have been studied by means of a classical dynamics, semi-classical approximation or variational methods [76][77][78][79][80]. In this paper we study the scattering problem for ion-dipole potential, focusing on the regime of very large α parameter: α 1. In such a case, the wave functions at r → 0 is singular, and one needs to impose some supplemental boundary conditions, defining the short-range behaviour of the wave function. This could be done, for instance, in the spirit of the quantumdefect theory (QDT) [81][82][83][84][85] , where one introduces some short-range parameters, that weakly depend both on the collision energy and on the angular momentum of the relative motion [24].
In general, such a treatment can be extended to the case of reactive scattering, where apart from the phase parameters, one additionally introduces amplitude of the short-range reaction processes [50,51]. For realistic collisions, the short-range QDT parameters, depend on the details of the short-range potential of a specific system. In this paper, we perform analysis of reactive scattering in the universal regime, when the reaction probability is equal to unity at short range. In this very special case, there is no outgoing probability flux at small distances, and the phase of the wave function at small distances is not important, hence, there is no need to include any additional QDT parameters.
The paper is organized as follows. In Section II we show how the Schödinger equation separates into the radial and angular parts. In Section III we derive some useful properties of modified spherical harmonics, which are later used in Sec. IV to calculate the scattering amplitude. The general formulas for elastic and reactive collision rates are derived in Sec. V and Sec. VI, respectively. In Sec. VII we investigate properties of modified spherical harmonics, while Sec. VIII is devoted to analysis of the radial solutions. Reactive scattering in the universal limit is discussed in Sec. IX. We conclude in Sec. X presenting some final remarks.

II. SEPARATION OF THE SCHÖDINGER EQUATION
We consider scattering of a molecule with a permanent electric dipole moment, and a charged particle, which could be for instance monoatomic ion as well as electron or positron. We assume that the dipole moment orientation is fixed in space, which is equivalent to solving the equations of motion in a body-fixed frame related to the polar molecule. The Schrödinger equation describing the wave function of the relative motion in the ion-dipole where µ is the reduced mass of the particles, V is the interaction potential between the ion with charge q and a polar molecule with dipole moment d. The square of the angular momentum operator is given bŷ The interaction potential is given by the scalar product of the dipole moment d and the electric field E of the ion: Denoting the angle between the vectors d and E by θ, we have Introducing E = 2 k 2 /2µ, the Schrödinger equation can be written as The dimensionless parameter α is given by: In physical systems composed of a polar molecule and a monoatomic ion, a permanent dipole moment is of the order of 1D (Debye). For such systems, α is usually a large number. From the point of view of recent experiment on ultracold systems, the most relevant are polar molecules of two-alkali metal atoms and alkaline earth metal ions. In such a case α ranges from α = 4.43 × 10 3 for LiNa-9 Be + up to α = 6.07 × 10 5 for LiCs-174 Yb + . In contrast, for electron/positron scattering on a molecule with a permanent dipole moment, α is typically of the order of one.
The system possess cylindrical symmetry, so the quantum number m associated with the z-component of the angular momentum is conserved. In general the wave function ψ(r) can be decomposed into the radial R(r) and angularỸ m (θ, φ) parts The centrifugal barrier and the dipol-ion potential fall off with the distance according to the same power law. Therefore the radial and angular part of the wave functions can be solved independently. We introduce the op-eratorÛ that describes the angular part of the stationary states Y m (θ, φ). We put an additional index to that function, which numbers the eigenvalues ofÛ . The solution of the eigenvalue problem gives the spectrum λ ,m and the corresponding eigen-functionsỸ ,m (θ, φ), which in the rest of the paper will be referred as modified spherical harmonics. We choose the numbering of the eigenvalues so that = |m|, |m| + 1, |m| + 2, . . .. Using this convention, in the limiting case of vanishing α, we recover λ , m → ( + 1) and the modified spherical harmonicsỸ ,m (θ, φ) reduce to standard spherical harmonics Y ,m (θ, φ). Similarly to standard spherical harmonic, we choose their normalization in the way, they form an orthogonal basis The next step is to solve the radial part of the Schödinger equation given by On sufficiently large distances, when the short range potential can be entirely neglected and only dipol-ion interaction is present, the radial part can be expressed in terms of spherical Bessel functions of the order given by a real or purely imaginary number.

III. RESOLUTION OF PLANE WAVE IN MODIFIED SPHERICAL HARMONICS
Before we analyse the scattering problem, we present two formulas that are used in derivation of the scattering amplitude. The first is the resolution of the identity operator: where n i are unit vectors along θ i and φ i (i = 1, 2). The vectors n 1 and n 2 can be interchanged on each side without affecting the sums. Also, the position of the complex conjugate is unimportant. This formula represents the fact that eigenvectors of operatorÛ form a complete orthogonal basis (for fixed value of α). The second formula is the expansion of the plane wave e ik.r in the basis of the modified spherical harmonics Y ,m . We start from the familiar expansion (see e.g. [75]) where j (x) is the spherical Bessel function of the first kind,k andr are unit vectors (denoted by hat) directed along k and r. For large values of r = |r| this expansion takes the following form: This can be rewritten as Here, we have used the parity of the spherical harmonics (−1) Y ,m (k) = Y ,m (−k). Emploing Eq. (12) in the second line in the above formula, we arrive at Note, that equivalently the minus sign in front ofk in the second term could be put in the front ofr. This formula cannot be further simplified, because modified spherical harmonicsỸ ,m (r) in general do not have specified parity.

IV. SCATTERING PROBLEM
In the scattering problem we solve the Schrödiner equation (1) with the boundary conditions where k i is the initial wavevector of the incident particle, k = kr and f (k i → k) is the scattering amplitude describing the scattering process. Note, that the above wave function is normalized such that the probability current of the incident particle (plane wave e iki·r ) is k i /µ, which is its velocity. Below, we express the amplitude f in modified spherical harmonicsỸ ,m and elements S ,m of the scattering S-matrix defined as follows. After solving the radial part of the Schrödinger equation (11), we find its asymptotic form for large r: where A ,m are given by the boundary conditions of the considered physical problem, and S ,m are the elements of the scattering S-matrix. Now, we will find the scattering amplitude. We start with the paricles wave function written in the basis R ,m (r) andỸ ,m (r) where R ,m (r) is given by (17). Setting the coefficients A ,m = 4π(−1) Ỹ * ,m (−k i ), we find the relation (16), where the scattering amplitude is given by In the following sections we present the formulas for elastic and reactive collision rate constants K el and K re .
The radial coordinate j r of the probability current for the term f e ikr /r, describing the scattering wave, is given by j r = k|f | 2 /(µr 2 ). The differential elastic cross section dσel/dΩ is defined by the following relation: is the differential part of the total cross section, that contributes to the scattering into the solid angle dΩ. By equating the number of particles scattered into the solid angle dΩ per unit time: j r dΩr 2 , with the number of particles in the incident particle probability flux is the velocity of incident particles, and we assume that probability density in the incident wave function is normalized to unity: ρ = |ψ| 2 = 1, in accordance with normalization of the wave function assumed in (16). Hence, dσ el (k i → k) = |f | 2 dΩ. The total elastic cross section is obtained by integrating dσ el (k i → k)/dΩ over all possible directionsk of the scattered particle, Inserting here the formula (19), we arrive at This final expression for the total cross section depends on the direction of the incident particlek i . It is thus natural to consider the cross section averaged over all possible directions of incidence. Therefore, the averaged elastic cross section is given by Note that theσ el does not depend only on S ,m but also on the form ofỸ ,m . This comes from the fact that the potential V id (r) is not spherically symmetric. This is mathematically expressed by the fact that the scalar product Ỹ * ,m (−k i )Ỹ ,m (k i )dΩ i in general depends on α, whereas in the case of spherical harmonics we have . Nevertheless, an important simplification comes from the presence of the cylindrical symmetry. Namely, the function σ el (k i ) depends only on angle θ i between the dipole moment d and k i . This can be seen from (20), where the phase φ enters only as a phase e imφi which is unimportant after taking the modulus squared. Consequently, in Eq. (21) only the integration over one variable θ i has to be performed.
Here, we will give an expression for the elastic collision rate constant K el , which is by definition given by the averaged elastic cross section and the velocity of the incident particle: An alternative formulation involves the probability flux j scatt r of the scattered particle. It is given by taken at the limit r → ∞. Following the normalization assumed in (16), we have assumed that the flux of the incident particles is equal to the velocity v i = k i /µ, i.e. one particle per unit area per unit time.
It can be shown that for a pure ion-dipole potential, the total scattering cross section is infinite, which is a consequence of its long-range character and the anisotropy. Mathematically, it is related to the fact, that in the expression for the cross section, given by Eq. (20), a mixed termỸ ,m (−k i )Ỹ * ,m (k i ) appears. Without the anisotropy, the modified spherical harmonics are equal to the standard ones, and this term is exactly (−1) after integration over directions ofk i . With the anisotropy, the term approaches the standard value, but not sufficiently fast, to make the integral convergent (see the details in Appendix C).

VI. REACTIVE COLLISION RATE K re
The reactive rate constant K re is most easily obtained by formulating it as a lost probability flux averaged over all directions. We rewrite Eq. (16) in terms of amplitudes f out (r) and f in (r): The radial component of probability flux corresponding to the outgoing and incoming particle are given by j out r = k|f out | 2 /(µr 2 ) and j in r = k|f in | 2 /(µr 2 ), respectively. The difference between these currents integrated over all possible directions of the scattered particle describes the rate of the probability loss due to the reactions during the scattering process. The reactive rate is given by that loss of the probability averaged over all directions of incidence (calculated at the limit r → ∞): It should be noted that the above formula resembles the equation (23). Using the formula for the scattering amplitude found above, see Eq. (19), the reactive rate can be expressed in terms of the scattering matrix elements S ,m according to the formula: Unlike the elastic rate, this equation is similar to the usual relation for the reactive rate expressed in spherical harmonics.

VII. SOLUTION OF THE ANGULAR PART
Below we present general features of the solutions for the angular part of the wave functions. The equation that we consider is expressed by Eq. (9). To proceed, we decomposeÛ in the basis of spherical harmonics. We expressỸ ,m as a linear combination of spherical harmon- Note that all the harmonics have the same quantum number m. The coefficients c follow from the diagonalization of the operatorÛ .
In order to solve the eigenvalue problem, we need to evaluate the matrix elements of this operator. The operatorl 2 is diagonal in the basis of Y ,m . The matrix elements of the cos θ are given by: with β m, = ( 2 − m 2 )/(4 2 − 1); the matrix elements are zero if l |m|.
Denoting now the eigenvalues ofl 2 by V = ( + 1), and the contribution from the dipole by D = −αβ m, , the matrixÛ takes the following form: The eigenproblem for this matrix with |m| is expressed by the equation Below we present the results of numerical calculations of the modified spherical harmonicsỸ , m in the regimes α ∼ 1 and α 1.

A. Angular orbitals for intermediate values of α
In this section we investigate the properties of angular wave functions numerically. To this end we diagonalize the matrixÛ given in Eq. (30), truncating it for final max = 4000 and evaluate the modified spherical harmonics for α = 1 and α = 10, assuming m = 0 .
The results are plotted in Figs. 1 and 2, where we display three-dimensional plots of |Ỹ ,m (θ, φ)| 2 as a function of θ and φ for = 0, 1, 2, 3. Since all the functionỸ ,m are proportional to e imφ , the shown plots possess cylindrical symmetry with respect to rotations about z-axis.
We observe that for α = 1 modified spherical harmonic look similarly to standard spherical harmonics, but they are slightly shifted toward upper half-plane, which is due to the attractive part of the ion-dipole potential. However, for α = 10 the orbitals are highly anisotropic, and they are no longer symmetric with respect to reflections z → −z. In the figure, the gray XY surface marks z = 0 plane, and the intersection of the orbital with that plane is shown by a white contour. In particular, the orbitals = 0, m = 0 for both α = 1 and α = 10 are strongly shifted to upper half-space.
In Fig. 3 we present dependence of the orbital = 0 and m = 0 on α. We note that for increasing value of α, orbitals become more anisotropic. Even for relatively low values α ≈ 1, the displacement of the orbital is significant compared to the isotropic case α = 0. The numerical calculations show that for large values of α, the lowest orbitalsỸ , m are localized close θ = 0. In such a case one can expand Eq. (9) around that point, and obtain approximate form of angular orbitals analytically. To this end, we introduce the function Θ(θ) by writing:Ỹ The problem of finding eigenvalues and eigenvectors ofÛ reduces to finding a solution for Θ ,m (θ) of the following equation: where λ ,m is an eigenvalue ofÛ . Note that in the equation m parameter enters only in m 2 term, so λ ,m and Θ ,m (θ) depend only on |m|.
We exploit now the fact that the functions are localized around θ ≈ 0 and expand Eq. (33) for small θ substituting: cos θ ≈ 1 − 1 2 α 2 and sin θ ≈ θ. This leads to The solutions that satisfy appropriate boundary conditions, i.e. is finite for θ = 0 and vanish for large θ, are given in terms of confluent hypergeometric function: where the quantum number n = − |m|=0,1,2,. . . , indexes the number of nodes of the angular wave function. The corresponding eigenvalues are given by: The spectrum ofÛ operator start at −α + √ 2α(|m| + 1). The lowest lying eigenvalues are evenly distributed with interval 2 √ 2α. The eigenvalues increase with growing |m|. For completeness, we give here the normalization constant, which is given by From the form of solution (35) it can be seen that the function Θ ,m (θ) is negligible for θ (2/α) 1/4 , and combining this with the condition θ 1 we obtain the necessary condition for the validity of the presented approximation If this parameter is small, the approximation that led to (35) and (36) is applicable. This condition, however, is not sufficient, because for large , the wave function extends to larger θ, and we violate the condition θ 1. From quasi-classical considerations (see subsection VII B 3) we obtain that the region of nonvanishing wave function is of the order of θ (1 + λ , m /α). This gives the required condition for The low-lying part of the spectrum, as indicated by Eq. (36), is linear in . In Fig. 4 we present the spectrum ofÛ for α = 3.65 × 10 4 for different values of m calculated numerically by solving Eq. (31). This value of α roughly corresponds to collisions of KRb polar molecule with 86 Sr + or 87 Rb + ions. We observe at small values of , that the spectrum exhibits linear behaviour, which is well described by Eq. (36) for m = 0 and m = 50. For m = 100 the slope of the linear dependence at small is slightly different than predicted by Eq. (36). At large values of m the spectrum is no longer linear with , and interestingly it becomes universal, not depending on the value of m.

Continuous-l approximation.
In this section we exploit the fact that for large α 1 the coefficients c in Eq. (27), which enter the recurrence relation (31), change smoothly with . The details of the derivation can be found in Appendix A, here, we merely state the final results.
First, we introduce the small parameter of the expansion, which is = α −1/4 . Then we drop indices and m, rename λ , m to λ and write c order terms being neglected): The solution that is finite both at x → 0 and for large x is given in terms of the confluent hypergeometric function: where n = 0, 1, 2, 3, . . . is the quantum number labelling the solutions, and the eigenvalues are given by We notice that eigenvalues in Eq. (42) are identical as in Eq. (36).

Quasi-classical approximation.
The last approach that we apply to solve the angular part of the Schrödinger equation in the regime α 1 is the quasi-classical approximation. Our starting point is Eq. (33). We introduce the function χ(θ) defined by with normalization condition π 0 |χ , m (θ)| 2 dθ = 1.
This new function χ satisfies the following exact equation, that resembles the Schrödinger equation: with the effective potential: where the eigenvalue ε , m = λ , m + α + 1 4 . Below, where no ambiguity arises, we drop the indices and m in χ and ε.
Eq. (45) is the stationary Schrödinger equation for a particle moving in potentialṽ(θ) with energy ε. In the quasi-classical approximation we introduce left and right classical turning points, denoted further by θ L and θ R , respectively, defined by k cl (θ L ) = k cl (θ R ) = 0, and θ L < θ R . The classical wavevector k cl is given by Applying the quasi-classical method to the radial Schrödinger equation requires inclusion of the so-called Langer correction [86], which basically boils down to dropping of 1/4 in the second term of Eq.
Within the quasi-classical approximation, the wave functions for the eigenstates in the classically accessible region are given by: for θ L < θ < θ R , and far from the ends of that interval. In the classically inaccessible region, the wave function are given by for θ > θ R , and for θ < θ L . For completeness we give the expression for the normalization constant, The eigenvalues ε can be obtained from the Bohr-Sommerfeld's quantization rule, given by: where n is a positive integer number indexing the n-th eigenvalue in the potential v(θ). Note that n = − |m|.
We have checked numerically that inclusion of the correction −1/4 significantly improves the accuracy of the approximate solutions. The procedure of solving the eigenproblem for the angular part, defined by equation (33), is now straightforward. We first find the shifted eigenvalues ε n indexed by a non-negative integer n by solving Born-Sommerfeld's quantization conditions (52). The wave functions are then given by equations: (49) -classically accessible region, (50) -right classically inaccessible region, and (51) -left classically inaccessible region.
In Fig. 5 we plot the angular parts of the wave function |χ , m | 2 as a function of θ for m = 20 and m = 100, and different = 35, 115, 120, 200. These functions were calculated numerically from Eq. (27), with the expansion coefficients c (m) , obtained by numerical solving of recurrence relation (31). Solutions are rescaled and shifted in order to fit into the quasi-classical potential v(θ). In Appendix B we derive formulas for the eigenvalues that are obtained from the Bohr-Sommerfeld's quantizations rule (52) for small |m|. In particular, we show that for small we recover the spectrum given by Eq. (36). For 0 ε < 2α we obtain the following equation deter-mining the eigenvalues ε with θ R = arccos(1 − ε/α). In the other regime, for the eigenvalues ε 2α we arrive at the following formula: Moreover, we obtain a closed formula for large : √ α, that is independent of m: This formula is depicted in Fig. 4 with a dotted line. We observe that for large eigenvalues corresponding to different m indeed reduce to a single m-independent curve given with a good approximation by Eq. (55). In Fig. 6 we display the spectrum λ , m for m = 20 and α = 3.65 × 10 4 . We compare the numerically calculated values from Eq. (31) to the ones obtained by solving the quasi-classical quantization rule, Eq. (52). We also display eigenvalues obtained by approximate quasi-classical quantization conditions, Eqs. (53) and (54). All the solution are in very good agreement with the exact numerical result. The relative error , which is shown in the inset, remains below 10 −3 % for the full quasi-classical formula and is larger for the quasi-classical solution without inclusion of the Langner correction, and for asymptotic formula (55).
In Fig. 7 we show the spectrum and relative errors in the case m = 100. Here the condition |m| α is not satisfied, and the potential v(θ) is strongly affected by the presence of the term m 2 / sin 2 θ in Eq. (48). As a consequence it cannot be neglected, and so the formulas given by Eqs. (53) and (54) are not as accurate as for m = 20.
Figs. 8 and 9 compare the wave functions |χ , m | 2 obtained from Eqs. (27) and (43), with expansion coefficients calculated numerically from Eq. (31) and the quasiclassical wave functions Eq. (49)- (51). Both approaches agree in the whole region except the neighbourhood of the classical turning points, where the quasi-classical approximation breaks down. The eigenvalues ε ,m for the quasi-classical wave functions were obtained from Bohr-Sommerfeld's quantization rule, Eq. (52), which works remarkably well, with relative errors smaller than 10 −3 % (see Figs. 6 and 7).
Finally, we remark that the case m = 0 should be treated with care. In the quasi-classical approximation, for m = 0 the potential v(θ) has no classical turning points around θ = 0 and θ = π. This can be traced back to dropping the contribution from − 1 4 sin −2 θ term in the Schrödinger equation. In particular, the quasi-classical approximation that would be required to vanish for θ = 0, where the potential v(θ) is finite, would have wrong  52)). Dotted blue lines (coinciding with the previous two) represent solutions given by Eqs. (53) and (54). The inset shows the absolute relative error between solutions given by numerical diagonalization and: blue dots -quasiclassical approximation given by Eq. (52) (smallest relative error), dot-dashed -quasi-classical approximation with omitted shift 1/4, i.e., λ = ε − α (medium relative error), and dotted line -the asymptotic form given by Eq. (55) (largest relative error). The spectrum λ , m ofÛ for m = 100 in units of α calculated numerically for α = 3.65 × 10 4 (dashed red line). Dot-dashed black line (coincides with full numerical dashed-red) is the spectrum obtained within quasi-classical approximation (see Eq. (52)). Dotted blue lines represent solutions given by Eqs. (53) and (54), below λ < α and above λ > α, respectively. The inset shows the absolute relative error between solutions given by numerical diagonalization and: red dots -quasi-classical approximation given by Eq. (52) (smallest relative error), red dot-dashed -quasi-classical approximation with ommited shift 1/4, i.e., λ = ε − α (medium relative error), and red dotted line -the asymptotic form given by Eq. (55) (largest relative error). phase in the classically allowed regime. Since for α 1 the contribution to collision rates comes typically from several partial waves with different values of m, here, we omit the quasi-classical analysis for m = 0 and, when needed, refer to numerical calculations.

VIII. SOLUTION OF THE RADIAL PART
We turn now to analysis of the radial part of the Schrödinger equation For the standard scattering problems, where the interaction potentials decays faster than r −2 , the centrifugal potential for all the partial waves except s-wave is always repulsive. However, this is not the case of ion-dipole potential, when eigenvalues λ ,m can be as well positive as negative. Negative values lead to attractive potential, which completely changes properties of the wave functions at short distances. Radial equation (56) can be solved with the Bessel functions of the first kind: where κ = λ , m + 1/4 depends on the indices and m. For large distances kr 1 we have whereas for small kr 1 we have The spectrum of the angular part of the wave function splits into two branches, which have different consequences in the radial part. We shift the eigenvalues as in the previous section, inroducing ε , m = λ , m + α + 1 4 .
It is instructive to analyse properties of the radial solutions in the limit of large angular momenta, which is equivalent to → ∞. Making use of the asymptotic formula (55), we obtain The second term can be neglected for α, and in this case the radial solutions are identical as for standard scattering problem with short-range potentials Different situation occurs for ε , m < α. In that case κ ,m is purely imaginary, so κ , m = i|κ , m |. For small kr 1 we have an oscillating solutions, R ± (r) ∝ (kr) −1/2 exp(±i|κ ,m | log(kr)). For large kr 1 we have the following asymptotic behaviour:

IX. REACTIVE COLLISIONS IN THE UNIVERSAL REGIME
In this section we employ the solutions of the radial equation to analyze reactive collisions in the universal regime [50,51]. To this end, we adopt the approach within QDT, which is based on the parametrization of the wave function at small distances, where with a good approximation, it does not depend neither on the energy nor on the angular momentum of the collision. Specifically, for the solution with imaginary κ, i.e., when λ ,m +1/4 < 0, the short distance wave function oscillates for r → 0, and The parts of the wave function with + and -describe outgoing and incoming probability currents, respectively. Within QDT we parametrize In general both the short-range phase φ and the parameter y, which describes the probability of short-distance reaction, can depend on and m. Below, we analyze the simplest possible case, i.e., we assume that y = 1. In this universal limit, the scattering properties do not depend on the short-distance phase φ. The scattering matrix is given by and is valid for λ ,m + 1/4 < 0.
In the other regime, if λ ,m + 1/4 > 0, κ is real and positive, and then the wave function at short distances takes the form Physically meaningful solution is obtained when assuming d 2 = 0, which, after straightforward calculation, leads to the following S matrix Taking into account (61), in the limit of large we obtain Now, we can proceed to evaluate the reactive collision rate given by Eq. (26). In general, partial waves for λ ,m > −1/4 do not contribute to the K re , since |S ,m | = 1 (cf. Eq. (68)). Physically, this corresponds to the scattering on the repulsive potential, where the particles do not approach the core region r = 0, where the reaction takes place. Therefore, the only contribution to K re comes from partial waves with λ ,m < −1/4. If α 1 we can get an estimate of the reactive collisions rate, since most of the contributing κ are large, and according to Eq. (66), S ,m ≈ 0. The reactive rate is then given by where * m is the maximum for which λ ,m < 0, and m * > 0 is the largest m for which negative λ ,m do exist. In the first approximation, as can be inferred from Eq. (36), m * = α/2 and * m = (m * + |m|)/2. This leads to Since the number of states, for which the particles can collide and react, is proportional to α, the reactive rate is also proportional to α. We note that the reactive rate depends on the energy as 1/k ∼ 1/ √ E. This behaviour agrees with the prediction for power-law potentials V (r) = −C n /r n [87] where in the limit of n → 2 we obtain substituting for probability of reaction P re = 1 and g = 1 for distinguishable particles. Alternatively, the formula from Eq. (71), can be derived based on classical considerations. To be specific, one should solve classical equations of motion and count only the contribution to the reactive rate from the trajectories that fall to the scattering center. By calculating contribution from such trajectories, one can reproduce exactly the formula (71).
In general, the universal collisional rate will have the following form where F is a universal function depending only on the dipol-ion interaction strength α. In Fig. 10 we plot the function F (α) that we calculate numerically. The function exhibits steps, when new states enter below the threshold λ ,m < −1/4. Since a finite value of α is needed to generate such a state, function F is nonzero for α α cr = 1.279. The α cr is the maximal value of α for which radial Schrödinger equation has solutions that are finite at r → 0 [68]. For large values of α the function F approaches the classical limit given by πα/4. In particular, for α 200, the error is smaller than 2%.

X. CONCLUSIONS
In this work, we have investigated reactive scattering in the ion-dipole potential. First, we have introduced modified spherical harmonics, describing angular wave functions in the presence of anisotropic ion-dipole interaction, and their corresponding eigenvalues. For large values of α parameter we have derived a number of analytical approximations based on expansion at small angles, continuous-quantum number approximation and the semi-classical method. We have shown that introduction of modified spherical harmonics, requires some modification of the formulas for elastic and reactive total cross sections and collision rates. We have investigated properties of the radial solutions, which are given in terms of Bessel functions, analysing two different regimes, corresponding basically to attractive and repulsive longrange interaction potentials. Finally, we have calculated the collision rates for reactive scattering in the universal regime, where the short-range reaction probability is equal to unity. In such a case it is not necessary to introduce any other phase parameters describing the shortrange behaviour of the wave function.
It is worth to emphasize that due to the 1/r 2 dependence, the ion-dipole potential exhibits a very peculiar features. Firstly, for such a potential one cannot define any kind of characteristic length scale or the energy scale, as can be done for other power-law potentials [87]. Therefore, the only characteristic parameter, that can be associated with this potential, is a dimensionless α parameter. Secondly, the local de Broglie wavelength is λ(r)/(2π) ∼ r/ √ α, therefore for large α the condition for the quasi-classical approximation: λ (r)/(2π) 1 [88], is fulfilled at all distances. In this case, there is no quantum reflection process at the intermediate distances, as happens for most power-law potentials at low energies [87]. Hence, the relative amplitude of the incoming and outgoing flux will be the same at short and large distances for all collision energies. This means that for α 1, reactive scattering can be very accurately described in the quasi-classical approximation. By summing all the contributions from different partial waves, it turns out that the total reactive cross section is identical to the cross section calculated in the framework of the classical physics.

XI. ACKNOWLEDGEMENTS
This is the Schrödinger equation for a shifted onedimensional harmonic oscillator in which the position is given by l and the wave function is c(l). The solution that converges for large l is given by where D ν (x) is the parabolic cylinder function with ν = −1/2 + (1 + 4α + 4λ)/(4 √ 2 √ α). The value of λ and thus of ν is determined by the boundary conditions for c(l) for small values of l where β m,l deviates from 1/2.
For large x, D ν (x) falls off exponentially, so c(l) becomes negligible for l √ α. Also, from the obtained solution it is clear that c(l) changes smoothly when incrementing l by a unit. These observations are the starting point for the following discussion. Below, we expand the equation (31) in the small parameter = α −1/4 . As in the main text, we drop here the indices and m, rename λ , m to λ and write c (m) , l as c(l). We change the variables x = l and c(l) =c( l) ≡c(x). We notice that for lowest levels λ ≈ −α = −1/ 4 . Thus we write λ =λ/ 4 , whereλ is of the order of unity. Now, we rewrite (31) in the following form: We expand the equation in the powers of , which leads to Eq. (40).

Appendix B: Quasi-classical analysis of the angular part
In this appendix we derive formulas for the eigenvalues that are obtained from the Bohr-Sommerfeld's quantizations rule, see Eq. (52), for small |m|. We show that for small we recover the spectrum given by Eq. (36). We also derive a closed formula for large √ α that is independent of m, Eq. (55). If |m| is not too large, the point θ 0 where α(1 − cos θ) is equal to m 2 / sin 2 θ, is much less than 1. If this is the case, then approximately θ 0 ≈ (2m 2 /α) 1/4 1. The following discussion is valid if m 2 α.
We first assume that 0 ε < 2 . Here, the position of the right turning point is mainly determined by the α(1 − cos θ) term in v(θ), whereas m 2 / sin 2 θ term gives only a small correction which we neglect, i.e., θ R = arccos(1 − ε/α). The position of the left turning point is mainly determined by the m 2 / sin 2 θ term in v(θ). Approximating sin θ ≈ θ we find θ L ≈ (m 2 /ε) 1/2 . We see that if ε |m| α/2, the left turning point is separated from the point where the two parts of the potential are of the same order, i.e., θ L θ 0 . We can use this separation, to effectively evaluate the dependence on m 2 from the integral. To see this, we notice, that the Born-Sommerfeld's quantization rules given by equation (52) can be rewritten in the following form: Because of the mentioned separation, we choose θ L such that θ L θ L θ 0 . In this regime in the first integral on the left-hand side we can neglect α(1−cos θ) and approximate m 2 / sin 2 θ with m 2 /θ 2 . Assuming that θ 2 L m 2 /ε, we can evaluate the first integral to θ L √ ε−|m|π/2. Now, in the second integral we can neglect the part m 2 / sin 2 θ of the potential, and write Here, the second integral cancels the θ L √ ε term from the first integral in (B1). Finally, for 0 ε < 2 we obtain with θ R = arccos(1 − ε/α).
As an application of formula (53), we evaluate the lowlying eigenvalues. For small θ , we expand 1 − cos θ ≈ 1 2 θ 2 in the integrand. The right turning point θ R = 2ε/α. Evaluation of the integral (53) is now straightforward. As a result we obtain the spectrum (neglecting the constant shift −1/4) that is the same as the one given in equation (36).
For the eigenvalues ε 2α we have to take care of the right turning point appropriately. Using the same reasoning to the right region that is accessible to the particle, as we described in the paragraph above, we arrive at the following formula: We see that the eigenvalues here depend only on the quantum number . We emphasize that these results are correct, if |m| is not too large, i.e., m 2 α. Note, that in the special case α = 0, we obtain ε = ( + 1/2) 2 , which gives exactly the correct value ( + 1) if we reintroduce the −1/4 shift to the eigenvalues of the angular part.
(C3) The integral is therefore where we used fact that Θ , m is real. Now we make first approximation, that the integration spans over the region that is classically accessible, neglecting the regions where the function decays exponentially. Therefore the left turning point is θ L = max(θ L , π − θ R ) and the right is θ R = min(θ R , π − θ L ).
We may note that always θ L < π/2, and so if θ R < π/2, the integral is negligible and we have I = 0. At least for small values of m, the condition θ R = π/2 is reached for ε ≈ α. Therefore, the integral is significantly nonzero only when ε > α. In this regime, θ L = π − θ R and θ R = θ R .