A Propagation Method

# Ultra-cold Collisions between Spin-Orbit-Coupled Dipoles: General Formalism and Universality

## Abstract

A theoretical study of the low-energy scattering properties of two aligned identical bosonic and fermionic dipoles in the presence of isotropic spin-orbit coupling (SOC) is presented. A general treatment of particles with arbitrary (pseudo-) spin is given in the framework of multi-channel scattering. At ultracold temperatures and away from shape resonances or closed-channel dominated resonances, the cross-section can be well described within the Born approximation to within corrections due to the -wave scattering. We compare our findings with numerical calculations and find excellent agreement.

## I Introduction

Many novel behaviors and phases in quantum few- and many-body systems can be understood from the competition between kinetic and interaction energies. The extraordinary degree of control in systems of ultracold quantum gases, therefore, provides versatile platforms to study these quantum phenomena Dalfovo et al. (1999); Leggett (2001); Bloch et al. (2008); Giorgini et al. (2008). For instance, short-range interactions (or more precisely, -wave scattering lengths) between atoms can be tuned to virtually arbitrary values via magnetic Feshbach resonances Chin et al. (2010), allowing us to access the unitary regime O’Hara et al. (2002); Makotyn et al. (2014); Hu et al. (2016) and study Efimov physics Kraemer et al. (2006); Berninger et al. (2011); Wang et al. (2012a). On the other hand, long-range interactions, especially the anisotropic dipole-dipole interactions, can significantly change the excitation spectrum Santos et al. (2003) and the stability diagrams of Bose-Einstein Condensation (BEC) Santos et al. (2000); Griesmaier et al. (2005); Koch et al. (2008), which has attracted intense interest in gases of ultracold heteronuclear ground state molecules Ni et al. (2008); Rvachov et al. (2017), and magnetic dipolar atoms such as Cr Griesmaier et al. (2005), Dy Lu et al. (2010, 2011), and Er Aikawa et al. (2012). More recently, manipulating kinetic energies and corresponding dispersion relationship has been realized via the innovative synthetic spin-orbit coupling (SOC) technique, i.e. coupling a particle’s canonical momentum with its (pseudo-) spin degrees of freedom Zhang et al. (2014); Zhai (2015). The realization of SOC provides an important ingredient for studying the fractional quantum Hall effect, topological insulators Hasan and Kane (2010); Qi and Zhang (2011), and has been a fundamental advancement in ultracold quantum gases in recent years Dalibard et al. (2011); Goldman et al. (2014).

There are currently several experimental techniques to realize SOC in cold-atom systems, such as lattice shaking Struck et al. (2012) and Raman coupling Lin et al. (2011). In particular, the Raman laser scheme has been applied to achieve one-dimensional SOC (an equal mixture of Rashba and Dresselhaus spin-orbit coupling) Lin et al. (2011); Cheuk et al. (2012); Wang et al. (2012b); Zhang et al. (2012a); Qu et al. (2013); Khamehchi et al. (2017) and two-dimensional SOC Huang et al. (2016); Meng et al. (2016); Wu et al. (2016) in ultracold gases of alkali atoms. However, Raman coupling for alkali atoms usually also comes along with atomic heating due to spontaneous emission. This heating leads to the loss of quantum degeneracy and trap population, which is a major challenge to study many-body quantum phenomena that manifest at a lower temperature and longer timescales. For atoms with higher ground state orbital angular momentum, spontaneous emission can be eliminated while still producing large Raman coupling Cui et al. (2013), making the open-shell lanthanide atoms Dy and Er suitable candidates. These elements also possess large magnetic dipole moment, allowing studies of the interplay between SOC and long-range dipole-dipole interactions that do not exist in alkali atomic gases, but requiring a more sophisticated theoretical model. Experimentally, SOC has been recently achieved by Ref. Burdick et al. (2016) in Dy, which allowed for the realization of a long-lived SOC degenerate dipolar Fermi gas. The bosonic system of SOC dipolar gas has also been theoretically investigated previously in BEC of Cr Deng et al. (2012).

The existence of long-range dipole-dipole interactions in these systems is expected to give an interesting interplay with the SOC, leading to intriguing new quantum phases. Previous theoretical studies on the interplay between short-range two-body interactions and SOC in Fermi gases have explored novel superfluid states in the BEC-BCS crossover Gong et al. (2011); Hu et al. (2011); Yu and Zhai (2011); Vyasanakere et al. (2011); Han and Sá de Melo (2012). For a BEC with SOC, new quantum phases, such as a stripes phase, have been predicted for a certain range of the Raman coupling strengths determined by the inter- and intraspecies scattering lengths Wang et al. (2010); Ho and Zhang (2011); Zhang et al. (2012b); Li et al. (2013). The anisotropic and long-range dipole-dipole interaction can be regarded as an additional degree of freedom, which might lead to new physics, but also brings new challenges in theoretical studies. To construct a concrete theoretical model for SOC dipolar quantum gases, the low-energy scattering between two dipoles in the presence of SOC needs to be understood first, which is the main topic of our study here.

Our theoretical formalism is inspired by several previous studies on ultracold collisions between two non-dipole particles in the presence of the three-dimensional (3D) isotropic SOC (which is a 3D analog of Rashba SOC) Cui (2012); Zhang et al. (2012c); Duan et al. (2013); Wang and Greene (2015); Guan and Blume (2016, 2017). While the 3D isotropic SOC has not yet been realized experimentally in cold-atom systems, proposals have been made that are based on adding more laser fields. Recently, the realization of 2D isotropic SOC using this scheme has been reported Sun et al. (2017). The laser scheme to realize 3D isotropic SOC is ideally suited for lanthanide atoms, where there is less atomic heating. On the other hand, the 3D isotropic SOC is more closely related to the cases in condensed-matter physics due to the high symmetry Wan et al. (2011); Burkov and Balents (2011). This symmetry also allows for a fully analytical treatment of low-energy scattering in the presence of SOC.

The low-energy scattering is strongly affected by the asymptotic behavior of atoms at large distances. In our system, the SOC persists even for atoms with infinite separation, which changes the threshold energies, and modifies the dispersion relation Dalibard et al. (2011); on the other hand, a dipole-dipole interaction also dominates potential energies at large distances. The competition between SOC and dipole-dipole interaction, therefore, gives rise to special threshold behaviors. The details of our formalism to study this problem is outlined below in Sec. II. An analysis by applying the first-order Born approximation is given in Sec. III. A comparison with numerical results for spin- dipolar fermions is given in Sec. IV. A summary of our study is given in Sec. V.

## Ii Formalism

In our model, each dipole is treated as a point particle with mass . The interaction potential between two dipoles aligned to the -axis and seperated by a large distance is therefore given by . Here is the polar angle of in spherical coordinates, and denotes the dipole moment, where is the magnetic dipole moment and is the vacuum permeability. The characteristic length scale of the dipole potential is given by the dipole length , where is the two-body reduced mass. Correspondingly, a natural energy scale can be defined by the dipole energy . To mimic the experimentally available control of short-range interactions by using methods such as Feshbach resonances, we model the short-range potential by a simplistic hard-wall potential, i.e. for , and for . This specific chosen form of the short-range potential, however, does not limit the generality of our study of ultracold collisions especially near potential resonances, which will be elaborated on later.

Similar to Refs. Cui (2012); Zhang et al. (2012c); Duan et al. (2013); Wang and Greene (2015); Guan and Blume (2016, 2017), we focus on the scattering in the center-of-mass frame. With the presence of 3D isotropic SOC, the Hamiltonian in relative coordinates is given by,

 H=→p22μ+kso2μ→p⋅(→s1−→s2)+V(→r), (1)

where and are the spin operators for atom 1 and atom 2, is the relative momentum operator, and is the strength of SOC in the units of inverse length. The energy scale for SOC can therefore be defined by the recoil energy .

Following the same spirit as Refs. Duan et al. (2013); Wang and Greene (2015); Guan and Blume (2016, 2017), we solve the relative Schrödinger equation formally as a multichannel problem, i.e. using channel functions (basis) of , all degrees of freedom except for , to expand the ’th independent solution as,

 Ψτ(→r)=∑νΦν(Ω)Fντ(r)r. (2)

The channel functions adopted here are the tensor spherical harmonics that are simultaneous eigenstates of whose eigenvalues are collectively represented by . Here, is the (relative) orbital angular momentum operator, is the total spin operator, is the total angular momentum, and is the projection to the -axis in the laboratory frame. The tensor spherical harmonics are defined as,

 Φν(Ω)≡⟨θ,ϕ|(ℓs)jmj⟩=iℓ∑mℓ,msCjmjℓmℓ;smsYℓmℓ(θ,ϕ)χ(s,ms), (3)

where are the Clebsch-Gordan coefficients, are the usual spherical harmonics, and denote the spin states. The phase term is introduced to make the matrix elements of the Hamiltonian all real, which will be convenient for carrying out numerical propogation later. The matrix elements for the first two terms in Eq. (1) (except for an additional phase associated with the term) have been derived previously in Refs. Duan et al. (2013); Wang and Greene (2015):

 ⟨(ℓ′,s′)j′m′j∣∣→p22μ∣∣(ℓ,s)jmj⟩≡ℏ22μ(−Iν′νd2dr2+B(2)ν′ν)=(−ℏ22μd2dr2+ℏ2ℓ(ℓ+1)2μr2)δjj′δmjm′jδℓℓ′δss′, (4)

and

 ⟨(ℓ′,s′)j′m′j∣∣kso2μ→p⋅(→s1−→s2)∣∣(ℓ,s)jmj⟩≡ℏ22μ(Aν′νddr+B(1)ν′ν1r) Unknown environment '% ×[−(−1)s√s1(s1+1)(2s1+1){s1s2ss′1s1}+(−1)s′√s2(s2+1)(2s2+1){s2s1ss′1s′2}] Missing or unrecognized delimiter for \left (5)

which are all real. Here, the curly bracket denotes the 6j symbol. Since the isotropic SOC preserves total angular momentum, different ’s are not coupled by these two terms. However, the anisotropic dipole-dipole interaction will couple different ’s, and only is still a good quantum number due to the azimutual symmetry. The matrix elements for dipole-dipole interaction are then given by,

 ⟨(ℓ′,s′)j′m′j∣∣Vd(→r)∣∣(ℓ,s)jmj⟩≡ℏ22μB(3)ν′ν1r3=−iℓ−ℓ′2d2r3(−1)ℓ′−ℓ+s+jδs′sΠℓjℓ′Cj′m′jjmj;20{ℓsjj′2ℓ′}(ℓ′2ℓ000), (6)

where . These matrix elements are also real despite the factor, since the 3-j symbol at the end of Eq. (6) ensure that . In addition, the Clebsh-Gordan coefficient shows that is a good quantum number and only channels with can be coupled (specially, if , only channels with are coupled.) In principle, one needs to include channels with all possible angular momenta for an exact calculation, however, only a finite number of channels, , are needed in practice to obtain a converged result for the scattering cross-sections, where is sufficiently large (about for our chosen parameters) con ().

In terms of these matrix elements, the Schrödinger equation in matrix form can then be written as,

 (−I–d2dr2+A––ddr+B(r)–––––−k2I–)F(r)–––––=0, (7)

where the underlined variables, , denote matrices with matrix elements , is the identity matrix, is the incident energy and . The logarithmic derivative matrix can be obtained by propogating from to a sufficiently large distance (about for our chosen parameters). The details of the propagation method is elaborated in Appendix A. The -matrix and -matrix are therefore given by,

 K––=(L––g–−g′––)−1(L––f––−f′––), (8)

and

 S––=(I–+iK––)(I–−iK––)−1, (9)

respectively, where and are the regular and irregular solutions in matrix form.

The regular solutions are obtained by projecting the plane-wave solutions onto the tensor spherical harmonics, Eq. (3). The plane-wave solution with scattering energy can be written as

 ⟨→r|ξ,ζ;+^k⟩=√12+2δξζ×[ei→kξζ⋅→r∣∣ξ,+^k⟩∣∣ζ,−^k⟩+(−)pbe−i→kξζ⋅→r∣∣ζ,−^k⟩∣∣ξ,^k⟩], (10)

where equals for identical bosons (fermions) and is a single-particle state with the projection of spin along the quantization axis being . In the presence of 3D isotropic SOC, the particle can be well described by its helicity state, where the quanization axis is along the direction of its canonical momentum. Here is the canonical momenta with direction and magnitude , where . The expansion gives the matrix elements of as , where is the normalization constant chosen to ensure flux density conservation, and

 uντ=√2ℓ+12j+1Cj,ξ−ζℓ,0;s,ξ−ζCs,ξ−ζs1,ξ;s2,−ζ1+(−)s1+s2−s+ℓ+pb√2+2δξζ. (11)

Correspondingly, the matrix elements of the irregular solutions are given by , where and are the regular and irregular spherical Bessel functions respectively. Hereafter, unless specified otherwise, we use to collectively represent the quantum numbers of the channel function in Eq. (3) and to represent the partial wave of a particular helicity state. Therefore, Eq. (11) can also be regarded as a unitary tranformation between the helicity basis denoted by and the spin singlet/triplet basis in the absence of SOC indicated by quantum number for a particular partial wave of . The explicit values of for spin- fermions and spin- bosons has been previously obtained for in Ref. Duan et al. (2013) and Ref. Wang and Greene (2015) respectively, which can be used to verify our Eq. (11) (after carefully taking care of the factor). The form of regular and irregular solutions guarantee the -matrix to be real and symmetric (and hence the -matrix to be unitary), where the proof is given in Appendix B.

Finally, the cross-section from one partial wave of helicity states to another partial wave of another helicity channel is given by,

 στ′τ=2πk2τ′|Sτ′τ−δτ′τ|2. (12)

## Iii First-order Born Approximation

One of the most important observables in ultracold collisions is the threshold law behavior determined by the competition between SOC, short-, and long-range interactions. In this work, we are mostly interested in magnetic dipolar atoms whose dipole lengths is about nm for different species. Therefore, we are focusing on the case , as the SOC strength is reasonably estimated to be . We remark that the parameter regime of might be achieved in systems of heteronuclear molecules or Rydberg atoms, which is, however, beyond the scope of this paper.

The threshold behaviors for dipolar scattering without SOC have been discussed in Refs. Bohn et al. (2009); Wang and Greene (2012); Bohn and Jin (2014); Zhang and Jie (2014), where the scattering cross-sections of partial waves with are universally determined by the dipole length. The physical explanation is that scattering at low energy can only occur at distances larger than the centrifugal barrier, where the potential is dominated by the dipole-dipole interaction. In addition, the behavior of dipole-dipole interaction is weak at large distances, which allows the application of the perturbative first-order Born approximation. We apply the first-order Born approximation within the multi-channel framework, where the -matrix can be approximated by Wang and Greene (2010); Wang et al. (2011),

 K(Born)τ′τ=π∫∑ν′νf∗ν′τ′(r)2d2~B(3)ν′νr3fντ(r)dr, (13)

where by comparing with Eq. (6). Inserting the expression of the regular solution , we arrive at,

 K(Born)τ′τ=4Dℏ2π2μkτ′kτ√Nτ′Nτ∑ν′νu∗ν′τ′~B(3)ν′νuντΓτ′τℓ′ℓ, (14)

where . Since we are in the perturbative regime, i.e., , the cross-section can be approximated by

 σ(Born)τ′τ≈2πk2τ′|2K(Born)τ′τ|2. (15)

We would like to remark here that due to the absence of a centrifugal barrier in the -wave channel, the first-order Born approximation should not apply to terms in the expansion when as the integral is divergent. However, this is not a problem, as for , and hence gives no contribution to . Similar to the argument for the non-SOC case Bohn et al. (2009), the -wave contributions can be included later by supplementing the Born approximation with a short-range contribution that can be determined from the full closed-coupling calculations.

The threshold behavior for partial waves satisfying the first-order Born approximations can be further explored by using the analytical properties of the integral :

 Γτ′τℓ′ℓ=(kτ′kτ)ℓ′πΓ[12(ℓ+ℓ′)]8Γ[12(3+ℓ−ℓ′)]Γ(32+ℓ′)×2F1[12(−1−ℓ+ℓ′),12(ℓ+ℓ′),32+ℓ′,k2τ′k2τ], (16)

for , where is the hypergeometric function. Due to the symmetry of the integral, one can obtain for by simply switching the primed and unprimed indices on the right-hand side. In addtion, we only need to focus on the cases with , since equals zero otherwise. The integral can therefore be further simplified for ,

 Γτ′τℓ′ℓ=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩12ℓ′(ℓ′+1),ℓ′=ℓ,16(ℓ′+1)(ℓ′+2),ℓ′=ℓ−2,16ℓ′(ℓ′−1),ℓ′=ℓ+2. (17)

For , the integral can also be simplified in the limit :

 Unknown environment '%' (18)

where the (double) exclamation marks denote a (double) factorial.

We remark here that, at the regime , , can therefore always be approximated by Eq. (17). Furthermore, if we investigate the case of (and therefore , and ), one can verify that agrees with the -matrix element found in Ref. Bohn et al. (2009); Bohn and Jin (2014) for dipole-dipole scattering without the presence of SOC.

## Iv Example: Two Spin-1/2 Fermionic Dipoles

We apply our analysis to systems of two identical spin-1/2 fermionic dipoles as an example, and focus on the and even channels. In order to simplify notification we use and to represent the helicity and respectively in this section. Furthermore, to avoid double counting, we always choose . Therefore, the three possible two-particle helicity states for the two dipoles are given by and , with canonical momentum , , and , and normalization constant , and . Using Eq. (11), we find that for , only and are involved and coupled to the channel functions and . For higher even ’s, all three possible helicity states are involved and coupled to , and . Therefore, only the partial waves with are coupled to -wave, and the first-order Born approximation can be applied to all other higher partial waves.

Some of the elastic scattering cross-sections and corresponding -matrix elements are shown in Fig. 1. The solid curves are numerical calculations using parameters and , while the dashed curves are the first-order Born approximations found from Eq. (14) and Eq. (15). The first-order Born approximation agrees excellently with the essentially exact calculations at low scattering energy (small ) and high angular momentum partial waves. The almost perfect agreement can be understood by realizing the scattering occurs at larger distances for lower scattering energy and a higher centrifugal barrier, where the first-order Born approximation becomes almost exact.

The first-order Born approximation also agrees well for the essentially numerically exact inelastic process. In Fig. 2, we present some inelastic cross-sections and the corresponding -matrix elements from channels of to , but the state remains in the same helicity, which also show very good agreement. We can write out the explicit form for the elastic and inelastic cross-sections whose incoming and outgoing two-particle helicity states are the same:

 σj,mj→j′m′j(−−)→(−−) =128πD2C2τ′τ(√k2+(kso/2)2+kso/2)2k2+(kso/2)2, (19) σj,mj→j′m′j(−+)→(−+) =128πD2C2τ′τ, (20) σj,mj→j′m′j(++)→(++) =128πD2C2τ′τ(√k2+(kso/2)2−kso/2)2k2+(kso/2)2, (21)

where is a constant with respect to , and is given by Eq. (17). The threshold power law at can therefore be given by , , and .

For inelastic scattering processes with different incoming and outgoing helicity states, simple formulas for cross-sections without involving the hypergeometric function do not exist. However, in the limit , the threshold behavior can still be analyzed by using Eq. (18). For example, in the process , we have in the low limit, which leads to . Therefore, the leading order of the -matrix element is , where is the lowest that can couple to , which equals 1 in this case. The additional factor of in the exponent comes from the factor . The cross-section, therefore, obeys the power law . The same analysis can be applied to other scattering processes, and are summarized in Fig. 3 for cross-sections from channels of to . We can see that the power-law describes the threshold behaviors well.

We have also carried out calculations for different , and observe that the cross-sections for discussed previously are insensitive to , i.e., the cross-sections at low scattering energy are universally determined by the dipole length and This universality also applies better for lower and higher angular momentum due to the better application of the first-order Born approximation. However, for the partial cross-sections in the subspace of , the universality implied by the first-order Born approximation no longer exist due to the absence of a centrifugal barrier. Indeed, we find that the cross-section in the subspace of depends on and can change by orders of magnitude in our numerical calculation, as shown in Fig. 4. One can also see that the cross-sections in the subspace share a lot of similarities with the short-range results presented in Ref. Duan et al. (2013), such as the power-law behaviors, and the identical cross-sections for the same final helicity state regardless of the initial helicity state at very small .

In addition, the cross-sections for final states at low scattering energy goes to a constant and can reach resonance by tuning , similar to the non-SOC situation studied in Ref. Kanjilal and Blume (2008). In particular, near the broad potential resonances found in Ref. Kanjilal and Blume (2008) (or more specificly, away from shape resonances and closed-channel dominated resonances), we have observed another universality, i.e., the cross-section can be universally determined by the dipole length , SOC coupling , and the non-SOC singlet -wave scattering length . This universality implies that all the short-range physics can be absorbed into one single parameter , and the detailed form of the short-range interaction is not important. The underlying physics is the classical suppression of the WKB wave function amplitude by the large potential well at short distances Wang et al. (2012a), where a frame transformation is allowed. Therefore, we do not expect this universality can be applied to shape-resonances or closed-channel dominated resonances, where details of short-range potential becomes important.

As shown in Fig. 5, for a fixed , the elastic cross-sections are calculated near two different resonances and can be fitted as a function of ,

 σj=0(−,−)→(−,−)/D2=σres/D2(D/as−D/ares)2/Γ2res+1, (22)

where , , and are fitting parameters, which only depend on . Figure 5 also shows the effect of on resonances: the resonance shifts further to the positive side and becomes broader for larger SOC coupling. Interestingly, the shift of resonance due to the presence of SOC can be explained by the interplay between the short-range interaction and SOC (see Appendix C for details).

Another interesting feature in ultracold scattering with the presence of SOC is that the particles are preferentially scattered into the lowest helicity states (where the particle’s momentum is antiparallel to its spin direction), regardless of their incidence channel Duan et al. (2013); Wang and Greene (2012). This spontaneous handedness is an analog of an antiferromagnetic phenomena induced by the momentum-dependent magnetic field. The presence of dipole-dipole interaction would not change this spontaneous handedness effect, as can be seen by comparing the ratios of the different scattering cross-sections . Therefore, after sometime, we expect all the particles in our system to be in a “” helicity state, and any rethermalization due to a perturbation should be described by the total cross-section of , where

 σ(Born)(−,−)→(−,−)=∑j,j′,mjσj,mj→j′mj(−,−)→(−,−). (23)

Summing the partial cross-section from Eq. (19) in the limit of gives . Noticing that implies that the total cross-section from the first-order Born approximation equals the average of cross-sections for identical fermions and identical bosons without the presence of SOC Aikawa et al. (2014); Bohn and Jin (2014). We believe this reflects the fact that the total cross-section sums over the singlet and triplet cross-sections, which corresponds to the non-spin identical bosons and fermions respectively. In addition, when the particles are in helicity, they have equal probability to be projected into singlet and triplet states.

## V Conclusion

In summary, this paper extends previous theoretical studies of ultracold scattering in the presence of 3D isotropic SOC to a dipolar system. Our formalism is general in the sense that it can be applied to either bosons or fermions with arbitrary spin and the inclusion of any angular momentum partial waves. Similar to the non-SOC cases, the cross-sections involving high angular momentum partial waves can be well described by the first-order Born approximation, and can be determined universally by the dipole length and spin-orbit coupling strength However, the cross-sections that can couple to an -wave channel depend on the short-range physics and can have resonances. Nevertheless, all the short-range physics can be described by one additional parameter, , near a broad potential resonance. We have tested our theory in the example system of spin- dipolar fermions, and find excellent agreement with our numerical calculations.

While this work focuses on the ultracold regime , our formalism can be easily extended to the negative scattering energy following the same approach in Ref. Guan and Blume (2016). In this energy regime, however, the canonical momentum should be understood as given by (which has the same definition of in Ref. Guan and Blume (2016)). This topic, however, will be addressed elsewhere in the future.

###### Acknowledgements.
We would like to thank Hui Hu for his insightful discussions. This research was supported under Australian Research Council’s Future Fellowships funding scheme (project number FT140100003) and Discovery Projects funding scheme (project number DP170104008). The numerical calculations were partly performed using Swinburne new high-performance computing resources (Green II).

## Appendix A Propagation Method

Based on the R-matrix propagation method using discrete variable representations (DVR) basis, we develop a numerically stable method for propagating the logarithmic derivative matrix , where is the solution of Eq. (7). Numerically, we seperate the whole regime into many sectors. For each sector the propagation method allows us to calculate the logarithmic derivative matrix at one end for a given on the other end. One key ingredient in this method is the construction of the DVR basis. Our DVR basis functions are in the form of Langrange polynomial:

 πi(r)=√1~wiN∏j≠ir−rjri−rj, (24)

where and are defined by the Gauss-Lobatto quadrature points and weights correspondingly, with and Rescigno and McCurdy (2000). One can easily verify that a DVR basis satisfies that leads to the DVR approximation, i.e. with a smooth function . In addition, the derivative of DVR basis can be derived analytically.

Under the DVR approximation, Eq. (7) can be written as,

 H––→c(μ)≡(T––+M–––+V––−k2I–)→c(μ)=Λ––→c(μ), (25)

after integrating by parts, where (in vector notation, ) are the expansion coefficients for the matrix elements of ,

 Fνμ(r)=∑jc(μ)νjπj(r). (26)

The matrix elements of other terms in Eq. (25) are given by

 Tμi,νj=δμν∑m~wmπ′i(rm)π′j(rm), (27)
 Mμi,νj=Aμν12∑m~wm[πi(rm)π′j(rm)−π′i(rm)πj(rm)], (28)
 Vμi,νj=Bμν(ri)δij, (29)

which are all symmetric. The surface term is given by

 Λμi,νj={πi(r)[δμνπ′j(r)−12Aμνπj(r)]}|a2a1. (30)

From the form of the surface term, we define a matrix , so that we have,

 F––′(r)=[L––(r)+12A––]F––(r), (31)

which gives,

 ∑jνΛiτ,jνc(μ)jν=∑ν[δiNδjN~wNLμν(a2)c(μ)νN−δi1δj1~w1Lμν(a1)c(μ)ν1]. (32)

Defining the matrix with the elements

 hcc′iτ,jν=H––iτ,jν+δi1δj1~w1Lμν(a1), (33)

where and are a collective index for selected DVR basis indices, i.e. and . The radial equation can therefore be written in a matrix form:

 (hsshsNhNshNN)(→c(μ)s→c(μ)N)=(000L(a2)/~wN)(→c(μ)s→c(μ)N), (34)