Low-temperature chemistry using the R-matrix method
LLow-temperature chemistry using the R-matrixmethod
Jonathan Tennyson, ∗ a Laura K. McKemmish, a and Tom Rivlin a Techniques for producing cold and ultracold molecules are enabling the study of chemical reac-tions and scattering at the quantum scattering limit, with only a few partial waves contributing tothe incident channel, leading to the observation and even full control of state-to-state collisionsin this regime. A new R-matrix formalism is presented for tackling problems involving low- andultra-low energy collisions. This general formalism is particularly appropriate for slow collisionsoccurring on potential energy surfaces with deep wells. The many resonance states make suchsystems hard to treat theoretically but offer the best prospects for novel physics: resonances arealready being widely used to control diatomic systems and should provide the route to steeringultracold reactions. Our R-matrix-based formalism builds on the progress made in variational cal-culations of molecular spectra by using these methods to provide wavefunctions for the wholesystem at short internuclear distances, (a regime known as the inner region). These wavefunc-tions are used to construct collision energy-dependent R-matrices which can then be propagatedto give cross sections at each collision energy. The method is formulated for ultracold collisionsystems with differing numbers of atoms.
To paraphrase the recent review by Stuhl et al. , a quiet revolu-tion is occurring at the border between atomic physics and exper-imental quantum chemistry. There has been a rapid developmentof techniques for producing cold and even ultracold moleculesthrough techniques such as photoassociation of ultracold alkaliatoms, buffer-gas cooling, Stark deceleration, evaporative cool-ing and laser cooling . This progress is now enabling the ex-perimental study of chemical reactions and scattering at the quan-tum scattering limit with only a few partial waves contributing tothe incident channel (e.g. Quemener and Julienne ). Moreover,the ability to perform these experiments with non-thermal dis-tributions comprising specific states enables the observation andeven full control of state-to-state collision rates in this regime.This is perhaps the most elementary study possible of scatteringand reaction dynamics. The trapping and subsequent study ofchemical reactions involving cold or ultracold chemically im-portant molecules, such as OH and CaH, has opened a wholerange of possibilities that can be explored in chemical and quan-tum mechanical control and exploitation. These experimentalbreakthroughs demand equally transformative theoretical meth-ods for treating ultra cold reactions; these are, for many cases,still lacking. a Department of Physics and Astronomy, University College London, London WC1E 6BT,UK ∗ E-mail: [email protected]
One important feature of ultracold reactions is the pronouncedstructures present in the cross sections due to temporary forma-tion of long-lived quasi-bound states of the compound system,known as resonances. Resonances are ubiquitous in ultracold col-lisions and offer the best opportunity for quantum control and steering: they are already used to steer the formation of ultra-cold diatomic molecules: see, for example, Malony et al. Fur-thermore, elastic and inelastic cross sections can dramaticallychange near resonances, which directly influences the effective-ness of sympathetic cooling and trap losses. These resonances canbe manipulated using magnetic and electric fields. Studying thestructures of resonances in ultracold systems has yielded interest-ing physics, such as chaos in Dyspronium atoms, universalscaling laws/ characteristic behaviour and, when three ormore bodies are involved, Efimov resonances.
There are al-ready a number of examples of novel many-body state physics such as Bose-Einstein condensates (BECs), Efimov trimers, aswell as experiments investigating the crossover region betweenthe superfluidity of bosons in BECs and the Cooper pairing offermions in Bardeen-Cooper-Schrieffer (BCS) theory.
The resonance structure of systems which form over deepwells in their potentials which support many bound states islikely to be particularly rich and thus offer the greatest poten-tial for transformative science. These deep wells also offer themost opportunity for deviations from previously identified uni-versal characteristics. Here we propose a formalism explicitly de-signed to study such systems.1 a r X i v : . [ phy s i c s . c h e m - ph ] S e p rom a theoretical perspective, gas-phase, quantum reactivescattering at room/high temperature is well studied. Time-dependent methods have proved to be particularly powerful forthese problems . However time-dependent methods struggle atultra-low collision energies because of the long collision times in-volved; they are particularly poor at treating resonances. Thereare time-independent methods available which have been usedto treat low-energy collisions. The general physics can often beelucidated using simplified model theories .More molecule-specific theories include, in particular, proce-dures which use hyperspherical coordinates and basis set meth-ods. These theories have been developed and applied to low-energy collisions; see Honvault et al , and Pradhan et al. , forexample. These methods have been used successfully to treat anumber of slow atom-diatom collision problems and are the clos-est in spirit to what is proposed here. However, the hyperspher-ical methods generally involve transforming the problem into aseries of adiabatic potentials for which solutions are then found.While this approach has proved numerically successful, it is notphysically motivated and ultimately involves approximations con-cerning the couplings between the curves which are hard to over-come.The idea behind the new proposed R-matrix method for heavyparticle scattering is the division of space into two regions: an(energy-independent) inner region where most of the physicstakes place, plus an outer region where the interactions are sim-ple. In this inner region, solutions can be obtained by adaptingstandard bound-state programs. The R-matrix is then constructedon the boundary between these regions. Energy-dependent solu-tions to the scattering problem can then be obtained rapidly bypropagating the R-matrix. First principles, or calculable, R-matrixmethods have proved outstandingly successful for studies of lightparticle collisions and are being increasingly used in nuclearphysics . However, such methods have yet to be systematicallyapplied to heavy particle collisions. R-matrix methods were ex-tensively used to study chemical reaction in the 1980s but, apartfrom proof-of-principle studies , these calculations simply used(outer-region) R-matrix propagation over the entire coordinaterange . The proposal here is fundamentally different and ismuch closer in spirit to the methods successfully used by manygroups to study electron collisions.Our R-matrix based formalism builds on the progress made invariational calculations of molecular spectra which are now be-ing used to obtain solutions up to and beyond dissociation forstrongly bound systems such as water and H + . Boththese systems support about a thousand bound vibrational statesand many hundreds of thousands of bound rotation-vibrationstates for which solutions are also being found . These vari-ational calculations provide wavefunctions for the whole systemat short internuclear distances. Indeed, resonances for water and H + have already been studied using these approaches anda complex absorbing potential.There are now a variety of variational nuclear motion methodsand related computer programs available for solving these prob-lems. Here we focus on the codes used within our group: specif-ically the new code D UO , designed for open-shell, coupled-state diatomic problems , DVR3D for three-atom problems and itsfour-atom relative WAVR4 , as well as the general polyatomiccode TROVE . Our group has significant experience with pro-ducing spectroscopic accuracy potential energy surfaces that aregenerally assumed to be essential for quantitative predictions ofultracold collision physics . Hutson presents an interest-ing counter-argument, demonstrating that if there are significantcouplings to inelastic channels, then the sensitivity of the finalcross-section to the details of the potential energy surface is re-duced as the peaks in the cross-section due to the poles producedby the resonances are suppressed.In this paper, we present our proposed methodology and illus-trate it for the case of atom-atom collisions, utilising the new D UO program to obtain the inner region wavefunctions. This sim-ple system allows a proof-of-principle demonstration of our pro-posed methodology. Furthermore, the availability of relevant the-oretical and experimental results will allow thorough bench-marking of our methodology. In particular, we are interested inexplaining the surprisingly high measured cross section for thequenching of metastable, excited argon atoms by ultracold ar-gon . Within the Born-Oppenheimer (BO) approximation, the solutionof the reactive scattering problem divides into two steps: con-struction of the global potential energy surfaces using standardquantum chemistry methodologies, and solution of the nuclearmotion problem on these surfaces to produce scattering cross-sections and other properties of the reaction. We will assume herean appropriate potential energy surface is already available andfocus on the second part of this problem. The desired ‘solution’for scattering problems is the probability of different processes ata given collision energy. Note that generally, the actual wavefunc-tions solving the relevant time-independent Schrödinger equationare not necessary; instead, the observable information is embed-ded in quantities like the phase shifts, scattering S matrix, the K and T matrices and the cross sections.At large separation between the colliding species, the full scat-tering problems can be represented in terms of partial waves. Thedistinguishing characteristic of cold and ultracold scattering prob-lems is that only a small number of these partial wave compo-nents are needed to obtain a very good approximation to the fullanswer. At short separation between the colliding species, a fewpartial waves are no longer sufficient to describe the physics, par-ticularly when the two species interact strongly, i.e. collide overa potential with a deep well. Instead of trying to use a largenumber of partial waves, we propose using an approach whichtreats these two regions separately using methods that are opti-mal for each region. Specifically, we utilise the powerful varia-tional nuclear motion programs discussed earlier to find collisionenergy-independent solutions to the inner region problem, ψ k ,using a single diagonalisation. These energy-independent inner-region wavefunctions are used to construct the so-called R-matrixat the boundary r = a which is given in standard formulations R i , j ( E , a ) = a ∑ k ω k , i ω k , j E − E k , (1)where i and j are the asymptotic channels, and k runs over theinner region solutions, and the ψ k functions have energy E k andamplitude on the R-matrix boundary ω k . The coordinate r is aradial coordinate which asymptotically goes to dissociation prod-ucts. Inner region solutions can be obtained explicitly in terms ofthis coordinate by, for example, working in Jacobi coordinates, orthe amplitudes can be obtained by use of a projection operator onthe boundary.Once the R-matrix has been constructed at r = a , the energy-dependent, but computationally simpler, outer region problemis solved to give K-matrices, from which scattering observables,such as cross sections and resonance parameters, can be deter-mined. Due to the computational simplicity of this propagation,this outer-region propagation can be performed on a fine gridof collision energies, essential to elucidate resonance structure.Note that the R-matrix propagation actually becomes simpler atcolder temperatures because the number of asymptotic channelsdecreases significantly. Figure 1 gives a schematic representationof this solution strategy. Below we develop the theory needed to solve a simple two atomcollision problem on a single potential energy curve. Such a the-ory might apply to ultracold Ar – Ar collisions. Note that whilemuch of this theory is standard, it is often given in atomic units(assuming electron scattering) , such that the reduced massterms, which are important for heavy particle collisions, are miss-ing.Treating the inner region as a diatomic system, we can writea molecular ro-vibrational Hamiltonian operator in the followingway: ˆ H J = − ¯ h µ d dr + ¯ h J ( J + ) µ r + ∑ i ≥ i (cid:48) V ii (cid:48) ( r ) , (2)where µ is the reduced mass of the system of two masses m and m : µ = m m m + m , (3) r is the internuclear separation, J is the total angular momen-tum of the molecule, and V ii (cid:48) ( r ) is an element of the matrix ofpotentials associated with the atomic channels, including the off-diagonal channel coupling elements. These couplings can alsoarise from effects such as spin-orbit interactions which can be rep-resented using a generalisation of rotational operator . At thisstage we are interested in both bound and continuum solution tothis problem.Within the R-matrix method a (hyper-)radius a is chosenwhere the R-matrix is defined and inner regions solutions are ob-tained. There is some flexibility over the choice of a , althoughour plan is for the inner region to contain regions where the po-tential well is deep. However we note that the R-matrix methodhas proved highly successful at finding diffuse, long-range boundstates which extend outside the inner region and such statesare expected near the dissociation limit of polyatomic systems . Solving the Schrödinger equation with the Hamiltonian definedin eq. (2) within a finite region requires the introduction of asurface term, L , known as a Bloch term , to retain Hermiticity.The expression for this term is: L = ¯ h µ δ ( r − a ) ddr , (4)where δ ( r − a ) is the standard Dirac delta function. To solve themolecular problem with the surface term, we introduce a set offunctions { χ Jn ( r ) } . These functions are obtained as eigensolutionsto the matrix problem ( χ Jn | ( ˆ H J + L ) | χ Jn (cid:48) ) = E Jn δ nn (cid:48) (5)where, as is conventional , rounded Dirac brackets have beenused to show that integration in the radial coordinate, r , onlyruns over the finite volume of the inner region, from to a . Theeigenvalues, E Jn , of this equation are usually referred to R-matrixpoles and their associated eigenfunctions are defined using χ Jn ( r ) = ∑ i ∑ j c Ji jn φ Ji j ( r ) , (6)where { φ Ji j ( r ) } is some basis, and the coefficients c Ji jn are deter-mined by the requirement that eq. (5) is diagonal. Since the J isa conserved quantum number, we may label all solutions with it.Final results require the summation over J , but at low energiessuch sums should converge rapidly.The indices i and j in eq. 6 run over the channels, and the basisfunctions within each channel. To isolate the contribution from asingle channel, one can sum over only the basis functions withinthat channel, j , by defining w in ( r ) = ∑ j c i jn φ i j ( r ) . (7)From this, elements of the R-matrix, R J , can be defined on theboundary using the heavy particle generalisation of eq. (1) R Jii (cid:48) ( E , a ) = ¯ h µ a ∑ n w Jin ( a ) w Ji (cid:48) n ( a ) E Jn − E , (8)where w Jin ( a ) is called the surface amplitude (since it is evalu-ated at the boundary), E is the scattering energy of interest, andthe sum is over all n , i.e. over all eigensolutions of eq. (5). Wenote that it is also possible to reformulate the problem to use areduced set of inner region solutions. Note that a single set ofinner region solutions are used to construct the R-matrix at r = a for all scattering energies, meaning that the inner region problemonly needs to be solved once, independent of how many ener-gies the final solutions are needed for. This is particularly usefulfor obtaining high-resolution plots of resonances as a function ofscattering energy.From the scattering energy, E , the scattering wave number, k ,can be defined as k = √ µ E ¯ h . (9)A similar definition exists, and can be obtained from the eigenen-ergies E Jn , for the wave numbers k Jn . These can be written into a3 = a Duo, DVR3D, etc R − m a t r i x bound a r y Variational nuclear motion code R−matrixpropagationPotential energy surface E , k Ψ k Inner region: r < a cross sectionsK−matricesparametersResonance
Long−range potentialRepeated at each energyEnergy−independent
Asymptotic expansion
Outer region: r > a r = a p Fig. 1
Schematic division of space illustrating the use of the R-matrix method. diagonal matrix k J .Defining the outer region wavefunctions for a given targetchannel at some point, r = a , as F Ji ( a ) , the R-matrix representsthe relationship between these functions and their derivatives: F Ji ( a ) = a ∑ i (cid:48) R Jii (cid:48) ( E , a ) dF Ji (cid:48) ( r ) dr (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) r = a , (10)where the sum runs over all channels.Propagating the R-matrix to large r allows the scattering prob-lem to be solved without the explicit need to evaluate the wave-function which, particularly in the presence of closed channels,can be a difficult task numerically and computationally.There are a number of means of propagating the R-matrix, in-cluding those due to Baluja, Burke and Morgan , due to Lightand Walker , and the software FARM (the flexible asymptoticR-matrix package) . As discussed below, we favour the useof the Light-Walker propagator. Furthermore, there are severalways of obtaining the asymptotic wavefunctions, F Ji ( r ) . Inthis work the asymptotic expansion of Burke and Schey is used.Generally speaking, asymptotic expansions follow the form F Ji ( r ) = ∑ i (cid:48) (cid:32) s Jii (cid:48) ( k Ji (cid:48) r ) + ∑ i (cid:48)(cid:48) c Jii (cid:48)(cid:48) ( k Ji (cid:48) r ) K Ji (cid:48) i (cid:48)(cid:48) ( E ) (cid:33) , (11)where both sums are over all channels, K Ji (cid:48) i (cid:48)(cid:48) is an element of the K-matrix, K J ( E ) , and s Jii (cid:48) and c Jii (cid:48)(cid:48) are elements of the matrices s J ( k J r ) and c J ( k J r ) respectively. These matrices are matrices of ‘sine-like’ and ‘cosine-like’ functions respectively, which are differentfor different channels. The Burke-Schey asymptotic expansionspecifies the form of these functions, and is discussed in detail inthe next section.The propagated R-matrix is then combined with the asymp-totic expansion to construct the K-matrix, which has the followingform: K J ( E ) = − s J ( k J r ) − r R J ( E , r ) ˙s J ( k J r ) c J ( k J r ) − r R J ( E , r ) ˙c J ( k J r ) , (12)where ˙s J ( k J r ) and ˙c J ( k J r ) are the derivatives with respect to r of s J ( k J r ) and c J ( k J r ) respectively, and r is evaluated at some large value, denoted a p .From the K-matrix, the S- and T-matrices are defined in thefollowing ways: S J = + i K J − i K J , (13) T J = S J − . (14)Note that while the definition of the S-matrix is general, theprecise definition of the T-matrix depends on the conventionadopted.The eigenphase for each channel, δ Ji ( E ) , is given by the inversetangent of K J ( E ) ’s eigenvalues: δ Ji ( E ) = arctan ( K Ji ( E )) , (15)where K Ji ( E ) is the i th eigenvalue of the K-matrix, associated withchannel i . This, in turn, gives the eigenphase sum for a givensymmetry ( J ): δ J ( E ) = ∑ i δ Ji ( E ) . (16)The total cross section at a given energy, σ tot ( E ) , can be ob-tained in a number of ways, including from the eigenphase sums: σ tot ( E ) = π k J max ∑ J = ( J + ) sin ( δ J ( E )) , (17)where J max is the maximum number of angular momentum val-ues (partial waves) considered. For the ultracold temperaturesbeing considered here, this can be a very small number, possiblya single channel. For multi-channel collisions, the cross sectionfor a transition from channel i to channel i (cid:48) is σ i (cid:48) i ( E ) = π k i J max ∑ J = ( J + ) | T Ji (cid:48) i | . (18) The computational procedure for solving the above equations isessentially made of three steps: (a) the inner region, (b) theboundary and (c) the outer region and asymptopia. The final part,step (c), can be written in a fairly general fashion, which should4ater for a variety of different systems. Therefore our aim in writ-ing the code which constructs the R-matrix on the boundary, step(b), is to make it rather general to allow for the incorporation ofa variety of inner region nuclear motion codes. So far, in practice,we have only used the diatomic code D UO . D UO is designed tocompute spectra for open shell diatomic molecules and allows forexplicit inclusion of coupled potential energy curves through theinclusion of spin-orbit and other coupling terms. D UO is designedto read in potential energy and coupling curves in a variety offormats, including simply as a grid of points. Here we have usedthe Ar – Ar potential of Patkowski et al , defined on a grid. D UO constructs a basis using a Hund’s case (a) representation, whichis then used to obtain a full variational solution of the problem.Further details and discussion can be found elsewhere .The recently published version of D UO is designed only totreat bound rovibronic states. The first task is therefore to ex-tend this to give wavefunctions for the discretised continuum inthe inner region. This is done by constructing a set of functions { φ Ji j } . These functions are intended to be the set of square inte-grable, linearly-independent basis functions which are completeover the [ r min , a ] range up to some appropriate maximum energywhich enter into eq (6). As these functions are used to providethe amplitude of the inner region function { χ Jn } on the boundaryat r = a , one has to be careful how these functions behave at thispoint.In practice, the basis functions are generated in a two step pro-cedure. An initial basis set, { ψ Ji j } , is generated by solving themolecular problem associated with the ro-vibrational Hamilto-nian, ˆ H J , of eq. (2). In solving this problem, an artificial wallis placed in the potential at some distance r wall ( r wall > a ). Testshave shown that use of a wall provides a good representation ofresonance states contained inside it, see Fig. 2. Provided the wallis placed far enough out, the { ψ Ji j } are effectively computationalapproximations of the eigenfunctions of ˆ H J , with each basis func-tion index j belonging to a channel i . Generally speaking, placingthe wall such that a was approximately ≈ of the way to r wall was found to be appropriate .Although the { ψ Ji j } basis is constructed by integrating over thefull inner region [ r min , r wall ] , the rest of the R-matrix constructionmethod involves integrating over the smaller [ r min , a ] region. Inthis region, the { ψ Ji j } are not eigenfunctions of the Hamiltonian,so H J , the matrix of the Hamiltonian ˆ H J in this basis, will not bediagonal when defined over this range. The non-diagonal Hamil-tonian matrix is constructed by using a forward finite differencenumerical differentiation method of order four to evaluate thekinetic term, and an implementation of Simpson’s rule for the nu-merical integration up to a .The wall cannot be placed at a because, by construction, thebasis functions have zero amplitude at the wall. At this point thereare two possible approaches, both of which have been tested byus.The earlier proof-of-principle study of an R-matrix approachto reaction dynamics by Bocchetta and Gerratt simply diag-onalised a generalised version of the eigenvalue eq. (5) by in-cluding the overlap matrix on the left-hand side. Alternatively,the { ψ Ji j } functions can be re-orthonormalised over the [ r min , a ] range, for which we found symmetric or Löwdin orthonormalisa-tion to be the most suitable. Eventually we decided to utilisea generalised eigenvalue scheme, constructing the { χ Jn } functionsdirectly out of the { ψ Ji j } . In both methods we constructed a matrixfor the Bloch operator in the original basis, L , using the aforemen-tioned forwards finite difference method of order four to computethe derivative at a .The generalised eigenvalue problem we arrived at was ( H J + L ) χ Jn = E Jn S J χ Jn , (19)where S J is the overlap, or Gramian matrix, whose elements aremade of all the possible inner products between the different ψ Ji j basis functions. Equation (19) is then solved using the LAPACK routine dsygv to obtain the { χ Jn } and E Jn required to construct theR-matrix.For each value of J , our new R-matrix code reads from a ver-sion of D UO , adapted to implement the potential wall, all of theeigenenergies and eigenfunctions of the molecular system (thenumber of which is user-specified in D UO ), the minimums V ii min of the potentials associated with each channel V ii ( r ) , the wall po-sition r wall , the masses of the atoms m and m , the range overwhich the eigenfunctions are defined and orthonormalised, r min and r max , the step size of the integration, ∆ r , and the zero pointenergy (zpe).This information is used to constuct the R-matrix on the bound-ary, as outlined above, and this is then propagated outwards usingthe Light-Walker propagation method to a point a p . The Light-Walker propagator takes the form of an iteration equation for theR-matrix between the values a and a p , by dividing the region intosub-regions with boundaries a s . The propagator is constructed inthe following way : we diagonalise the matrix V J ( r ) = V J ( r ) − (cid:16) E J (cid:17) + E I , (20)where I is the identity matrix, E is the scattering energy, E J isthe diagonal matrix of eigenenergies (not to be confused withthe scattering energy), and V J ( r ) is the (in general) non-diagonalmatrix of potentials for each channel, including channel cou-pling elements (defined properly below in eq. (29) – note the J -dependence). We call the version of this matrix which has beenevaluated at a s and diagonalised (cid:0) v Js (cid:1) , and the matrix which di-agonalises it we call O Js : (cid:16) O Js (cid:17) T V J O Js = (cid:16) v Js (cid:17) . (21)This allows us to define the real, diagonal matrix λ Js in thefollowing way: (cid:16) λ Js (cid:17) = µ ¯ h (cid:18) E I − (cid:16) v Js (cid:17) (cid:19) . (22)Next we define elements of the diagonal matrix G Js , made up5 Energy (cm-1) r ( Å ) r W a l l a Fig. 2
Solutions of the inner region problem, χ n , for Ar showing both bound and continuum functions. of the following Green’s functions: G Jis ( r , r (cid:48) ) = − λ Jis sin ( λ Jis δ a s ) × (cid:40) cos ( λ Jis ( r (cid:48) − a s )) cos ( λ Jis ( r − a s − )) a s − ≤ r ≤ r (cid:48) cos ( λ Jis ( r − a s )) cos ( λ Jis ( r (cid:48) − a s − )) r (cid:48) ≤ r ≤ a s , (23)where δ a s = a s − a s − . Then, defining G Js ( r , r (cid:48) ) as G Js ( r , r (cid:48) ) = O Js G J (cid:16) O Js (cid:17) T , (24)we can write down the expression for the propagation equation: a s R Js = G Js ( a s , a s ) − G Js ( a s , a s − ) (cid:16) G Js ( a s − , a s − ) + a s − R Js − (cid:17) − G Js ( a s − , a s ) , (25)The size of each step in the iteration is variable, and dependenton the size of the last step. It obeys its own iteration equation,dependent on the derivative of the long-range potential used, inthe following way : δ a s + = β N N ∑ i = (cid:16) λ Ji , s (cid:17) − (cid:16) λ Ji , s − (cid:17) δ a s − / , (26)where i counts over the channels, N is the number of channelsand, β is a control parameter which allows you to specify howmany steps should be taken. β is currently taken to be . . Thevariable step size ensures that for different potentials, the appro-priate number of steps will be used to balance computation speedagainst accuracy. It also means that in the multi-channel case,channels which contribute different amounts can be treated dif-ferently. The initial step size is taken to be . of the distance from a to a p .Next we introduce the Burke-Schey expansion, a specifc versionof Eq. (11). In the Burke-Schey expansion, the matrices s J and c J have the following forms: s Jii (cid:48) = A Jii (cid:48) · sin ( k Ji (cid:48) r ) , c Jii (cid:48) = B Jii (cid:48) · cos ( k Ji (cid:48) r ) , (27)where A Jii (cid:48) = p max ∑ p = α Jpii (cid:48) r − p , B Jii (cid:48) = p max ∑ p = β Jpii (cid:48) r − p , (28)and the alpha and beta coefficients are derived from recurrencerelations. For Ar , the test system being studied, the long-rangeAr potential is written as V Jii (cid:48) ( r ) = V ii (cid:48) ( r ) + ¯ h J ( J + ) µ r = λ max ∑ λ = a λ ii (cid:48) r − λ + ¯ h J ( J + ) µ r , (29)and the diagonal coefficients a λ ii are obtained from Patkowskiand Murdachaew . Note in general λ max can vary for differentchannels. The α Jpii (cid:48) and β Jpii (cid:48) coefficients are then obtained fromthe following interdependent recurrence relations : (cid:18)(cid:16) k Ji (cid:17) − (cid:16) k Ji (cid:48) (cid:17) (cid:19) α Jpii (cid:48) + (( p − )( p − ) − J ( J + )) α Jp − , ii (cid:48) + k Ji (cid:48) ( p − ) β Jp − , ii (cid:48) = N ∑ i (cid:48)(cid:48) = λ max ∑ λ = a ii (cid:48)(cid:48) λ α Jp − λ − , i (cid:48)(cid:48) i (cid:48) , (30)and (cid:18)(cid:16) k Ji (cid:17) − (cid:16) k Ji (cid:48) (cid:17) (cid:19) β Jpii (cid:48) + (( p − )( p − ) − J ( J + )) β Jp − , ii (cid:48) − k Ji (cid:48) ( p − ) α Jp − , ii (cid:48) = N ∑ i (cid:48)(cid:48) = λ max ∑ λ = a Jii (cid:48)(cid:48) λ β Jp − λ − , i (cid:48)(cid:48) i (cid:48) ,. (31)6here N is the number of channels and λ max is the largest valueof λ (with larger values increasing both accuracy and computa-tion time). The derivatives of the s J and c J matrices also generaterelated recurrence relations, which can be derived by differenti-ating their power expansions.Finally, the coefficients obtained from the recurrence relationsare used to construct the asymptotic expansion. This expansion iscombined with the R-matrix to form the K-matrix using eq. (12).From this K-matrix the eigenphases are then obtained, and fromthe eigenphases the cross sections are obtained. As an initial test of the inner region codes we intend to use as in-puts to our new R-matrix code, we have looked for the so-calledshape resonances trapped behind the centrifugal barrier in therotationally excited Ar problem. Tests were performed for therotational state J = , which is sufficiently excited for the Ar system to not support any truly bound states. For this value of J , however, we would expect some quasibound states to existbehind a potential barrier, and for the Argon dimer potential ofPatkowski and Murdachaew , the peak of the centrifugal barrieris at . − . Calculations were performed by inputting thisAr potential on a grid both with Le Roy’s diatomic code LEVEL and with D UO .The D UO results were obtained by computing bound and con-tinuum energy levels both with and without a wall using a sta-bilisation procedure . In the case of a wall energy levels wereobtained in D UO up to v = at J = with the wall placed atvarious locations. A plot of various energy levels against wall lo-cation was then constructed, and the places on that plot wherethe energy levels appeared to lie on a horizontal line (i.e. wherethe energy did not vary with wall location) were used determinethe energies of the shape resonances. This is because the contin-uum energies follow a particle-in-a-box type energy distribution,and so are dependent on the size of the ‘box’. This is not the casefor actual shape resonances.A similar method was used to obtain the results without a wall,only instead of plotting energy against wall position, a plot ofenergy against the position of the end of the grid was used. Again,energy levels which did not vary with grid size were taken to beresonances, as opposed to a continuum state.As Table 1 shows, all three methods give two resonance energylevels. Furthermore, all three methods agreed to an accuracy bet-ter than . − . The fact that in both D UO cases the hori-zontally aligned energies agreed with the LEVEL results was anencouraging indicator that resonance energies had indeed beenfound. Table 1 Ar shape resonance energies obtained with three differentinner region solution methods at J = . All energies are in units of cm − . LEVEL D UO N Wall No wall1 7.7126 7.7126 7.71262 24.6178 24.6178 24.6178Figure 2 shows a further test of the inner region codes used as inputs for the R-matrix method, this time in the form of Ar innerregion wavefunctions obtained for the inner region problem us-ing the Patkowski and Murdachaew potential curve in D UO . Thewall was placed at 16 Å and the inner region boundary, a at15.045 Å. We note that all bound state wavefunctions, exceptthe highest one, are completely confined well inside our innerregion boundary. The highest state’s non-zero amplitude on theboundary is suspect and probably due to residual numerical is-sues. This means that the first 8 states make no contribution tothe R-matrix on the boundary and can be dropped from consider-ation in the scattering region. For polyatomic systems, droppingthe truly bound states from consideration should lead to substan-tial computational savings and could lead to the use of methodswhich do not compute wavefunctions for these states in the firstplace. Our aim is to use the R-matrix formulation of scattering theory tostudy many-particle problems and, by extension, chemical reac-tions. As discussed above, the methodology should be particularlyappropriate for ultra-low energy chemical reactions. Figure 3 il-lustrates how this should work for the prototypical reaction H + D + → HD + H + . (32)This exothermic reaction is likely to display significant resonanceeffects at very low energies. We note that Fig. 3 implies there is achange in coordinates, as the two asymptotes are most naturallyrepresented in different sets of Jacobi coordinates. How preciselythis is best achieved has to be determined, although one possibil-ity would be to solve the inner region problem in one coordinatesystem, here the higher symmetry H – D + Jacobi coordinates, forexample, and then use a projection operator on the boundary toconstruct the R-matrix on the boundary for the other coordinates,in this case HD – H + .In practice, our proposed methodology is by no means limitedto reactive scattering. Table 2 illustrates a number of the possibil-ities, again using the H D + system as an example. +2 eg H D HD(v",J") + H Products
R−matrix boundary scatteringReactants Inner Region (In)elasticDVR3D J ac ob i c oo r d i n a t e s : r( HD ) , R ( HD − H ) , H (v’,J’) + D
Outer Region eg H (v,J) + D
Jacobi coordinates: r(H ), R(H − D), θ θ Fig. 3
Schematic illustrating of the use of the R-matrix method. able 2 Processes that could be studied using a generalised R-matrixcode: the H D + system is simply used as an example, and not allpossible processes or products are listed. Process ExampleReactive scattering D + + H → HD + H + Photodissociation H D + + h ν → HD + H + or H + D + Photoassociation H + D + → H D + + h ν Charge Exchange D + + H → D + H + Elastic collisions D + + H ( v , J ) → D + + H ( v , J ) Inelastic collisions D + + H ( v (cid:48)(cid:48) , J (cid:48)(cid:48) ) → D + + H ( v (cid:48) , J (cid:48) ) Predissociation Not important for H D + The division of space into two regions raises a number of inter-esting possibilities. So far our studies on Ar – Ar collisions havesimply used the same potential energy curves in the inner andouter regions. However, the potential could be divided into tworegions: an inner region potential which captures the full com-plexity of the reaction, a complex intermediate potential, and along-range, outer region potential which can be represented usingknown multipolar forms for the dissociation fragment. Clearly,the two forms should match on the boundary. Standard quantumchemistry methodologies can be used to produce these potentials.Another possibility is to include very weak effects only in theouter region. For example, first principles studies of molecularspectra routinely neglect hyperfine effects. However, these areimportant at ultra-low collision energies. In our proposed methodone could re-couple hyperfine-free inner region wavefunctions onthe boundary so that the outer region problem fully incorporatesthese effects. This approach has been successfully used to treatspin-orbit effects in electron – light atom collision problems formany years . A similar approach could indeed be used to in-clude spin-orbit effects, which are usually totally quenched instrongly bound closed-shell systems, but become important whendissociation occurs to open shell species, such as what happens inwater . Use of the outer region in this fashion offers very signif-icant simplification compared to treating the full problem at allinternuclear separations.Experimental study of ultracold molecules is significantly en-hanced through the use of electric and magnetic fields to tuneresonances and thus increase production rates of the ultracoldmolecules ( ). Similarly, weak field effects could potentiallybe included in a similar fashion in the outer region. We note,however, that Zeeman effects have also recently been included inD UO . . The possibility of studying cold and ultracold collisions processes,and in particular chemical reactions, is one of the most interestingdevelopments of this century. These experiments are stimulatingwhole new areas of scientific investigation, e.g. in quantum con-trol, cold collisions, cold chemistry, accurate measurement, testsof fundamental physics, and more. Thus far most ultracold chem-istry studies have been on alkali metal dimers. Looking to thefuture, the next major stride will involve reactions of chemicallysignificant species and many atoms. Particularly important fornovel aspects of ultracold physics will be the exploitation of res- onances: long-lived quasi-bound states of the compound system.The development of theoretical approaches, for example the onedescribed in this paper, are essential for predicting, interpreting,and modelling this new physics. The R-matrix approach is de-signed to predict this interesting quantum behaviour and simu-late and support experimental studies in a rigorous and flexiblemanner, both theoretically and computationally.Our aim is to construct the harness code which links the innerand outer region segments. Initially this will be an atom-atomcode used for testing numerical and algorithmic aspects of theprocedure; some of these results are presented here. This workwill be used to guide the developments for larger collision sys-tems. The atom-atom code will also be used to study ultra-lowenergy collisions between systems being studied experimentally,starting with the Ar – Ar system mentioned above; this will allowus to explore the treatment of problems with coupled potentialsand magnetic fields, and extend our work to other nuclear motionmethods.
Acknowledgments
This project has received funding from the European Union’sHorizon 2020 research and innovation programme under theMarie Sklodowska-Curie grant agreement No 701962 and fromthe EPSRC.
References
Annu. Rev. Phys. Chem. ,2014, , 501–518.2 B. K. Stuhl, M. T. Hummon, M. Yeo, G. Quemener, J. L. Bohnand J. Ye, Nature , 2012, , 396–400.3 E. S. Shuman, J. F. Barry and D. DeMille,
Nature , 2010, ,820–823.4 V. Zhelyazkova, A. Cournol, T. E. Wall, A. Matsushima, J. J.Hudson, E. A. Hinds, M. R. Tarbutt and B. E. Sauer,
Phys. Rev.A , 2014, , 053416.5 G. Quemener and P. S. Julienne, Chem. Rev. , 2012, ,4949–5011.6 J. D. Weinstein, R. deCarvalho, T. Guillet, B. Friedrich andJ. M. Doyle,
Nature , 1998, , 148–150.7 V. Singh, K. S. Hardman, N. Tariq, M.-J. Lu, A. Ellis, M. J. Mor-rison and J. D. Weinstein,
Phys. Rev. Lett. , 2012, , 203201.8 D. S. Jin and J. Ye,
Chem. Rev. , 2012, , 4801–4802.9 C. Chin, R. Grimm, P. Julienne and E. Tiesinga,
Rev. Mod.Phys. , 2010, , 1225.10 M. Mayle, G. Quemener, B. P. Ruzic and J. L. Bohn, Phys. Rev.A , 2013, , 012709.11 S. Ospelkaus, K.-K. Ni, D. Wang, M. H. G. De Miranda,B. Neyenhuis, G. Quéméner, P. S. Julienne, J. L. Bohn, D. S.Jin and J. Ye, Science , 2010, , 853–857.12 P. K. Molony, P. D. Gregory, Z. Ji, B. Lu, M. P. Köppinger, C. R.Le Sueur, C. L. Blackley, J. M. Hutson and S. L. Cornish,
Phys.Rev. Lett. , 2014, , 255301.13 J. L. Bohn,
Phys. Rev. A , 2001, , 052714.14 J. M. Hutson, M. Beyene and M. L. González-Martínez, Phys.Rev. Lett. , 2009, , 163201.85 T. Köhler, K. Góral and P. S. Julienne,
Rev. Mod. Phys. , 2006, , 1311.16 T. Maier, I. Ferrier-Barbut, H. Kadau, M. Schmitt, M. Wenzel,C. Wink, T. Pfau, K. Jachymski and P. S. Julienne, Phys. Rev.A , 2015, , 060702(R).17 P. S. Julienne, Nature , 2014, , 440–441.18 J. P. D’Incao, H. Suno and B. D. Esry,
Phys. Rev. Lett. , 2004, , 123201.19 T. Lompe, T. B. Ottenstein, F. Serwane, A. N. Wenz, G. Zürnand S. Jochim, Science , 2010, , 940–944.20 F. Ferlaino, A. Zenesini, M. Berninger, B. Huang, H.-C. Nägerland R. Grimm,
Few-Body Systems , 2011, , 113–133.21 Y. Wang and B. D. Esry, Phys. Rev. Lett. , 2009, , 133201.22 J. P. D’Incao and B. D. Esry,
Phys. Rev. A , 2006, , 030703.23 T. Kraemer, M. Mark, P. Waldburger, J. G. Danzl, C. Chin,B. Engeser, A. D. Lange, K. Pilch, A. Jaakkola, H.-C. Nägerl et al. , Nature , 2006, , 315–318.24 S. Knoop, F. Ferlaino, M. Mark, M. Berninger, H. Schöbel, H.-C. Nägerl and R. Grimm,
Nat. Phys. , 2009, , 227–230.25 J. von Stecher, J. P. D’Incao and C. H. Greene, Nat. Phys. ,2009, , 417–421.26 F. Ferlaino and R. Grimm, Physics , 2010, , 9.27 M. Berninger, A. Zenesini, B. Huang, W. Harm, H.-C. Nägerl,F. Ferlaino, R. Grimm, P. S. Julienne and J. M. Hutson, Phys.Rev. Lett. , 2011, , 120401.28 I. Bloch, J. Dalibard and W. Zwerger,
Rev. Mod. Phys. , 2008, , 885.29 M. Greiner, C. A. Regal and D. S. Jin, Nature , 2003, , 537–540.30 S. Jochim, M. Bartenstein, A. Altmeyer, G. Hendl, S. Riedl,C. Chin, J. Hecker Denschlag and R. Grimm,
Science , 2003, , 2101–2103.31 M. Bartenstein, A. Altmeyer, S. Riedl, S. Jochim, C. Chin, J. H.Denschlag and R. Grimm,
Phys. Rev. Lett. , 2004, , 120401.32 T. Bourdel, L. Khaykovich, J. Cubizolles, J. Zhang, F. Chevy,M. Teichmann, L. Tarruell, S. Kokkelmans and C. Salomon, Phys. Rev. Lett. , 2004, , 050401.33 P. Pellegrini, M. Gacesa and R. Côté, Phys. Rev. Lett. , 2008, , 053201.34 S. C. Althorpe and D. C. Clary,
Ann. Rev. Phys. Chem. , 2003, , 493–529.35 R. T. Pack and G. A. Parker, J. Chem. Phys. , 1987, , 3888–3921.36 J. M. Launay and M. Le Dourneuf, Chem. Phys. Lett. , 1989, , 178–188.37 P. Honvault, M. Jorfi, T. González-Lezana, A. Faure and L. Pa-gani,
Phys. Rev. Lett. , 2011, , 023201.38 G. B. Pradhan, N. Balakrishnan and B. K. Kendrick,
J. Phys. B:At. Mol. Opt. Phys. , 2014, , 135202.39 P. G. Burke, R-Matrix Theory of Atomic Collisions: Applicationto Atomic, Molecular and Optical Processes , Springer, 2011.40 J. Tennyson,
Phys. Rep. , 2010, , 29–76.41 P. Descouvemont and D. Baye,
Rep. Prog. Phys. , 2010, , 036301.42 C. J. Bocchetta and J. Gerratt, J. Chem. Phys. , 1985, , 1351–1362.43 R. B. Walker and J. C. Light, Ann. Rev. Phys. Chem. , 1980, ,401–433.44 H. Y. Mussa and J. Tennyson, J. Chem. Phys. , 1998, ,10885–10892.45 G. H. Li and H. Guo,
J. Mol. Struct. (Theochem) , 2001, ,90–97.46 A. G. Császár, E. Mátyus, L. Lodi, N. F. Zobov, S. V. Shirin,O. L. Polyansky and J. Tennyson,
J. Quant. Spectrosc. Radiat.Transf. , 2010, , 1043–1064.47 N. F. Zobov, S. V. Shirin, L. Lodi, B. C. Silva, J. Tennyson,A. G. Császár and O. L. Polyansky,
Chem. Phys. Lett. , 2011, , 48–51.48 J. R. Henderson and J. Tennyson,
Chem. Phys. Lett. , 1990, , 133–138.49 J. R. Henderson, J. Tennyson and B. T. Sutcliffe,
J. Chem.Phys. , 1993, , 7191–7203.50 M. J. Bramley, J. W. Tromp, T. Carrington and G. C. Corey, J. Chem. Phys. , 1994, , 6175–6194.51 J. J. Munro, J. Ramanlal and J. Tennyson,
New J. Phys , 2005, , 196.52 T. Szidarovszky, A. G. Csaszar and G. Czako, Phys. Chem.Chem. Phys. , 2010, , 8373–8386.53 S. Miller and J. Tennyson, Chem. Phys. Lett. , 1988, , 117–120.54 R. Jaquet and T. Carrington, Jr.,
J. Phys. Chem. A , 2013, ,9493–9500.55 A. A. Kyuberis, O. L. Polyansky, L. Lodi, J. Tennyson, R. I.Ovsyannikov and N. Zobov,
Mon. Not. R. Astron. Soc. , 2016.56 T. Szidarovszky and A. G. Csaszar,
Mol. Phys. , 2013, ,2131–2146.57 B. C. Silva, P. Barletta, J. J. Munro and J. Tennyson,
J. Chem.Phys. , 2008, , 244312.58 S. N. Yurchenko, L. Lodi, J. Tennyson and A. V. Stolyarov,
Comput. Phys. Commun. , 2016, , 262–275.59 J. Tennyson, M. A. Kostin, P. Barletta, G. J. Harris, O. L.Polyansky, J. Ramanlal and N. F. Zobov,
Comput. Phys. Com-mun. , 2004, , 85–116.60 I. N. Kozin, M. M. Law, J. Tennyson and J. M. Hutson,
Comput.Phys. Commun. , 2004, , 117–131.61 S. N. Yurchenko, W. Thiel and P. Jensen,
J. Mol. Spectrosc. ,2007, , 126–140.62 A. Yachmenev and S. N. Yurchenko,
J. Chem. Phys. , 2015, , 014105.63 O. L. Polyansky and J. Tennyson,
J. Chem. Phys. , 1999, ,5056–5064.64 S. V. Shirin, O. L. Polyansky, N. F. Zobov, R. I. Ovsyannikov,A. G. Császár and J. Tennyson,
J. Mol. Spectrosc. , 2006, ,216–223.65 O. L. Polyansky, N. F. Zobov, I. I. Mizus, L. Lodi, S. N.Yurchenko, J. Tennyson, A. G. Császár and O. V. Boyarkin,
Phil. Trans. Royal Soc. London A , 2012, , 2728–2748.96 M. Pavanello, L. Adamowicz, A. Alijah, N. F. Zobov, I. I.Mizus, O. L. Polyansky, J. Tennyson, T. Szidarovszky and A. G.Császár,
J. Chem. Phys. , 2012, , 184303.67 J. M. Hutson,
New J. Phys , 2007, , 152.68 P. Barletta, J. Tennyson and P. F. Barker, New J. Phys , 2009, , 055029.69 P. D. Edmunds and P. F. Barker, Phys. Rev. Lett. , 2014, ,183001.70 J. Tennyson, L. Lodi, L. K. McKemmish and S. N. Yurchenko,
J. Phys. B: At. Mol. Opt. Phys. , 2016, , 102001.71 D. A. Little and J. Tennyson, J. Phys. B: At. Mol. Opt. Phys. ,2013, , 145102.72 C. Bloch, Nucl. Phys. , 1957, , 503.73 J. Tennyson, J. Phys. B: At. Mol. Opt. Phys. , 2004, , 1061–1071.74 K. L. Baluja, P. G. Burke and L. A. Morgan, Computer Phys.Comm. , 1982, , 299–307.75 L. A. Morgan, Computer Phys. Comm. , 1984, , 419–422.76 J. C. Light and R. B. Walker, J. Chem. Phys. , 1976, , 4272–4282.77 V. M. Burke and C. J. Noble, Computer Phys. Comm. , 1995, , 471–500.78 P. G. Burke, C. J. Noble, A. G. Sunderland and V. M. Burke, Physica Scripta , 2002,
T100 , 55–63.79 P. G. Burke and H. M. Schey,
Phys. Rev. , 1962, , 147.80 M. Gailitis,
J. Phys. B: At. Mol. Opt. Phys. , 1976, , 843.81 K. Patkowski, G. Murdachaew, C. M. Fou and K. Szalewicz, Mol. Phys. , 2005, , 2031–2045.82 A. T. Patrascu, C. Hill, J. Tennyson and S. N. Yurchenko,
J.Chem. Phys. , 2014, , 144312.83 P.-O. Löwdin,
J. Chem. Phys. , 1950, , 365–375.84 E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel,J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling,A. McKenney and D. Sorensen, LAPACK Users’ Guide , Societyfor Industrial and Applied Mathematics, Philadelphia, PA, 3rdedn, 1999.85 E. B. Stechel, R. B. Walker and J. C. Light,
J. Chem. Phys. ,1978, , 3518–3531.86 R. J. Le Roy, LEVEL 8.0 A Computer Program for Solving theRadial Schrödinger Equation for Bound and Quasibound Levels ,University of Waterloo Chemical Physics Research Report CP-663, http://leroy.uwaterloo.ca/programs/, 2007.87 H. S. Hazi, A. U.and Taylor,
Phys. Rev. A , 1970, , 1109.88 H. E. Saraph, Computer Phys. Comm. , 1978, , 247–258.89 O. V. Boyarkin, M. A. Koshelev, O. Aseev, P. Maksyutenko, T. R.Rizzo, N. F. Zobov, L. Lodi, J. Tennyson and O. L. Polyansky, Chem. Phys. Lett. , 2013, , 14–20.90 G. Quéméner and J. L. Bohn,