General Mathematical Formulation of Scattering Processes in Atom-Diatomic Collisions in the RmatReact Methodology

TL;DR

This paper extends the R-matrix method to model complex atom-diatomic collision processes, including reactive scattering, by developing a comprehensive mathematical framework suitable for ultracold collision scenarios.

Contribution

It introduces a general mathematical formulation for reactive and non-reactive scattering within the R-matrix framework, addressing the challenges of finite inner regions and complex coordinate sets.

Findings

Formulated the R-matrix approach for reactive scattering processes.

Addressed the mathematical challenges of finite domain integrals.

Extended the methodology to include charge exchange and photo-processes.

Abstract

Accurately modelling cold and ultracold reactive collisions occuring over deep potential wells, such as \ce{D+ + H2 -> H+ + HD}, requires the development of new theoretical and computational methodologies. One potentially useful framework is the R-matrix method adopted widely for electron-molecule collisions which has more recently been applied to non-reactive heavy particle collisions such as Ar-Ar. The existing treatment of non-reactive elastic and inelastic scattering needs to be substantially extended to enable modelling of reactive collisions: this is the subject of this paper. Herein, we develop the general mathematical formulation for non-reactive elastic and inelastic scattering, photo-association, photo-dissociation, charge exchange and reactive scattering using the R-matrix method. Of particular note is that the inner region, of central importance to calculable R-matrix methodologies, must be finite in all scattering coordinates rather than a single scattering coordinate as for non-reactive scattering. % The choice of coordinate set and basis function is these cases becomes more complexThis introduces substantial challenges to the basis sets utilised in practical calculations as integrals over finite domains are often much more challenging than over infinite domains for this problem.

Figures

Tables

Table 1. Table 1: Processes that can be studied using the RmatReact code , using H 2 D + system as an example.

Process	Example
Elastic collisions	D+ + H${}_2{}^{}(v,J)$ $\to$ D+ + H${}_2{}^{}(v,J)$
Inelastic collisions	D+ + H${}_2{}^{}(v^{\prime \prime },J^{\prime \prime })$ $\to$ D+ + H${}_2{}^{}(v^{\prime },J^{\prime })$
Photodissociation	H2D+ + $h\nu$ $\to$ HD + H+ or H2 + D+
Photoassociation	H2 + D+ $\to$ H2D+ + $h\nu$
Charge Exchange	D+ + H2 $\to$ D + H${}_{}{}^{+}{}_2{}^{}$
Reactive scattering	D+ + H2 $\to$ HD + H+

Equations

$$\Psi (E)=\sum _i{A}_i(E){\psi }_i$$

$${\psi }_{\text{total}}=\sum _{M\Omega J}{\psi }_J^{vib}{D}_{M\Omega }^{J\ast }(\alpha ,\beta ,\gamma ).$$

$${\mathcal{L}}^A=\frac{1}{2}\delta (r_2^A-a_0^A)\frac{d}{d{r}_2^A}$$

$$(T+V+{\mathcal{L}}^A+{\mathcal{L}}^B+{\mathcal{L}}^C){\psi }_i=E_i{\psi }_i$$

$${\psi }_{\text{total}}=\sum _c^{\text{channels}}\frac{1}{r_2}{F}_c(r_2){\Phi }_c$$

$$\left(-\frac{{\hslash }^2}{2{\mu }_2}\frac{{\partial }^2}{\partial r_2^2}+\frac{{\hslash }^2}{2{\mu }_2}\frac{l(l+1)}{r_2^2}+(e_{nj}-E)\right)F_c(r_2)=$$

$${\psi }_{\text{total}}=\sum _{c_X}\frac{1}{r_2}{F}_{c_X}(r_2){\Phi }_{c_X}+\sum _{c_A}\frac{1}{r_2}{F}_{c_A}(r_2){\Phi }_{c_A}.$$

$${\psi }_{\text{total}}=\sum _{c_A}\frac{1}{r_2^A}{F}_{c_A}(r_2^A){\Phi }_{c_A}+\sum _{c_B}\frac{1}{r_2^B}{F}_{c_B}(r_2^B){\Phi }_{c_B}+\sum _{c_C}\frac{1}{r_2^C}{F}_{c_C}(r_2^C){\Phi }_{c_C}$$

$$|\Psi (E)⟩=\sum _{i=1}^Z|{\psi }_i⟩\frac{1}{E_i-E}⟨{\psi }_i⟩\mathcal{L}\Psi (E)$$

$$F_{c}(a_{0})=\sum_{c}R_{cc^{\prime}}(E)a_{0}\frac{dF_{c}}{dr_{2}}\Bigg{|}_{r_{2}=a_{0}}$$

$$R_{c{c}^{\prime }}(E)=\frac{1}{2a_0}\sum _{i=1}^Z\frac{{\omega }_{c,i}{\omega }_{c^{\prime },i}}{E_i-E}$$

$${\omega }_{c,i}=⟨\frac{{\varphi }_c}{r_2}⟩{{\psi }_i}_{r_2=a_0}^{\prime }$$

$${\omega }_{njl,i}^{JMϵ}(a_0)$$

$$=$$

$${\mu }^{\text{trans-pd}}(h\nu )$$

$$=$$

$${\mu }^{\text{trans-pa}}(h\nu )=\sum _i{A}_i(h\nu -D_0+E_0)⟨{\psi }_f⟩\mu {\psi }_i,$$

$$|\Psi (E)⟩=\sum _{i=1}^n|{\psi }_i⟩\frac{1}{E_i-E}⟨{\psi }_i⟩{\mathcal{L}}_A+{\mathcal{L}}_B+{\mathcal{L}}_C\Psi (E),$$

$$F_{c_{G}}^{G}(a_{0}^{G})=\sum_{c^{\prime}_{A}}R_{c_{G}c^{\prime}_{A}}(E)a_{0}^{A}\frac{dF_{c^{\prime}_{A}}}{dr_{2}^{A}}\Bigg{|}_{r_{2}^{A}=a_{0}^{A}}+\sum_{c^{\prime}_{B}}R_{c_{G}c^{\prime}_{B}}(E)a_{0}^{B}\frac{dF_{c^{\prime}_{B}}}{dr_{2}^{B}}\Bigg{|}_{r_{2}^{B}=a_{0}^{B}}+\sum_{c^{\prime}_{C}}R_{c_{G}c^{\prime}_{C}}(E)a_{0}^{C}\frac{dF_{c^{\prime}_{C}}}{dr_{2}^{C}}\Bigg{|}_{r_{2}^{C}=a_{0}^{C}}.$$

$$R_{c_G{c}_{G^{\prime }}^{\prime }}(E)=\frac{1}{2a_0^{G^{\prime }}}\sum _{i=1}^Z\frac{{\omega }_{c_G,i}{\omega }_{c_{G^{\prime }}^{\prime }i}}{E_i-E}$$

$${\omega }_{c_G,i}=⟨\frac{{\varphi }_{c_G}}{r_2^G}⟩{{\psi }_i}_{r_2^G=a_0^G}^{\prime }$$

$$U_{c_G,c_{G^{\prime }}^{\prime }}(r_2)=\sum _{\lambda =1}^{{\lambda }^{\text{max}}}{a}_{c_G,c_{G^{\prime }}^{\prime },\lambda }(r_2^G{)}^{-\lambda -1}{\delta }_{G,G^{\prime }}$$

$$\{M_m(r_1)L_n(r_2)P_{j,k}(\theta )\}$$

$$\{M_{mA}(r_1^A)L_{nA}(r_2^A)P_{jA,kA}({\theta }^A),M_{mB}(r_1^B)L_{nB}(r_2^B)P_{jB,kB}({\theta }^B),M_{mC}(r_1^C)L_{nC}(r_2^C)P_{jC,kC}({\theta }^C)\}$$

$${\int }_0^{\pi }{\int }_0^{\infty }{\int }_0^{a_0^A}F(r_1^A,r_2^A,{\theta }^A)H(a_0^B-r_2^B)d{r}_1^{A}d{r}_2^{A}d{\theta }^A$$

$${\mathcal{Y}}_{j\bar{\Omega }}^{JM}(\widehat{{\mathbf{r}}_2},\widehat{{\mathbf{r}}_1})=\sqrt{\frac{(2J+1)}{4\pi }}{Y}_{j\bar{\Omega }}(\theta ,\gamma )D_{M\Omega }^J(\alpha ,\beta ,0),$$

$${\mathcal{Y}}_{j\bar{\Omega }}^{JM}(\widehat{{\mathbf{r}}_2},\widehat{{\mathbf{r}}_1})=P_{j\bar{\Omega }}(\theta )D_{M\Omega }^{J\ast }(\alpha ,\beta ,\gamma ).$$

$${\mathcal{Y}}_{j\bar{\Omega }}^{JMϵ}(\widehat{{\mathbf{r}}_2},\widehat{{\mathbf{r}}_1})=\frac{1}{\sqrt{2(1+{\delta }_{\bar{\Omega }0})}}\left({\mathcal{Y}}_{j\bar{\Omega }}^{JM}(\widehat{{\mathbf{r}}_2},\widehat{{\mathbf{r}}_1})+ϵ(-1{)}^J{\mathcal{Y}}_{j(-\bar{\Omega })}^{JM}(\widehat{{\mathbf{r}}_2},\widehat{{\mathbf{r}}_1})\right)$$

$${}^{\text{BF}}{\Phi }_{nj\bar{\Omega }}^{JMϵ}(\widehat{{\mathbf{r}}_2},{\mathbf{r}}_1)={\chi }_n(r_1){}^{\text{BF}}{\mathcal{Y}}_{j\bar{\Omega }}^{JMϵ}(\widehat{{\mathbf{r}}_2},\widehat{{\mathbf{r}}_1})$$

$${\Psi }_{\alpha }^{JMϵ}=\sum _{nj\bar{\Omega }}\frac{1}{r_2}{}^{\text{BF}}{F}_{i,nj\bar{\Omega }}^{JMϵ}(r_2){}^{\text{BF}}{\Phi }_{nj\bar{\Omega }}^{JMϵ}(\widehat{{\mathbf{r}}_2},{\mathbf{r}}_1)$$

Full text

††journal: rsta

1 Introduction

The rapid development of techniques for producing cold and even ultracold molecules over the last decade is now enabling the study of chemical reactions and scattering at the quantum scattering limit with only a few partial waves contributing to the incident channel. Moreover, the ability to perform these experiments with nonthermal distributions comprising specific states enables the observation and even full control of state-to-state collision rates in this regime. This is perhaps the most elementary study possible of scattering and reaction dynamics [1]. These experiments are driving the development of new theory to address the new physics encountered in these ultra-slow collisions. Reactions involved charged species are of special interest, in part because of the extra experimental control possible for charge particles [2].

Due to the strength of the \ceH3+ system interaction, the reaction of H+ with H2 can be particularly expected to show quantum behaviour at ultralow collision energies (temperatures); this has meant that the H${}_{3}^{+}$ system has become a benchmark system for the study of ultracold reactions [3, 4, 5, 6, 7, 8, 9]. Processes of interest are described in Table 1 using the \ceH2D+ system as an example. Other near-dissociation properties of H${}_{3}^{+}$ that merit study include the near-infrared photodissociation spectrum which was extensively chararacterised by Carrington, McNab and co-workers [10, 11, 12, 13, 14], but which remains poorly understood [15]. The long-range H+ – H2 potential has been shown to support diffuse states which were called asymptotic vibrational states [16]. These states have some similar characteristics to halo states found in diatomic systems [17] but unlike the diatomic systems are likely to be present in significant numbers although the actual density and structure of these states remains to be determined. Finally, formation of H${}_{3}^{+}$ by radiative association could be important in diffuse environments such as the early Universe. There are low temperature measurements of this rate [18] but no theory and no studies at very low temperatures.

Procedures based on the use of hyperspherical coordinates [19, 20, 21] have been developed to solve the close-coupling equations, and applied to treat ultra-low energy reactive scattering in D++H2 see Lara et. al. [9] and references therein.

RmatReact is promising new methodology that addresses this problem by solving an initial energy-independent problem that encapsulates most of the complicated physics in an inner region. These inner region solutions are then used to describe the scattering-dependent problem in the simpler outer region. This approach can be employed study not just the reactive scattering process but, in principle, all the processes in Table 1. The RmatReact methodology described here is a spiritual successor to the extremely successful electron-molecule collision codes such as UKRMol [22, 23], and indeed our implementation reuses part of this code base. Both methodologies separate space into three regions based on the distance between the two scattering partners: the inner, outer and asymptotic region. The inner region contains the region with significant non-multipole, and indeed non-local, interaction between the two scattering species; the Schrodinger equation solved here does not consider the scattering energy. The outer region contains the region where the interaction between the two scattering species is significant but simple in form (usually a multipole expansions); the outer region equations depend on the scattering energy and are solved once for every different scattering energy under consideration using the inner region solutions and a simple 1D propagation process. The asymptotic region is defined as the region in which the interaction between the scattering species is much less than the scattering energy and thus negligible; at this point, we can calculate properties of the overall scattering interaction such as cross-sections.

The mathematical formulation of our RmatReact methodology for non-reactive single and multi-channel scattering was previously presented by Tennyson et. al. [24]. Two initial applications of this new RmatReact calculable R-matrix method for heavy-particle scattering to atom-atom collisions over the Morse potential and Ar-Ar scattering were presented by Rivlin et. al. [25, 26].

Here, we extend the mathematical formulation of the RmatReact methodology to triatomics and consider for the first time photo-association, photo-dissociation, charge exchange and reactive scattering.

2 Mathematical Formulation

2.1 Overview of general RmatReact methodology

Our methodology models scattering of atoms or molecules A and B with center of mass separation distance $R_{AB}$ using the following steps:

Variationally solve the (3$N$-3)D Schrodinger equation of the joint system AB in a finite region at zero scattering energy (including Bloch terms to ensure the Hermicity of the Hamiltonian) where $N$ is the total number of nuclei in the system and $R_{AB}$ less than some box size $a_{0}$: this produces a discrete number, $Z$, of inner-region energies, $E_{i}$, and wavefunctions, $\psi_{i}$. 2. 2.

Map the (3N-3)D problem onto a 1D Hamiltonian in $R_{AB}$, with a reduced potential, $U(R_{AB})$. Simultaneously, each inner-region wavefunctions, $\psi_{i}$, can be mapped onto outer region channels, $\phi_{c}$ producing the surface amplitudes, $\omega_{c,i}$. 3. 3.

Using the energy-independent solutions, for each scattering energy under consideration, construct the scattering-energy-dependent R-matrix (to be defined below) at the boundary, thereby using the first set of solutions as an effective basis for describing the desired second set. 4. 4.

For each scattering energy, propagate this R-matrix to an asymptotic distance, at which point scattering observables such as cross-section can be evaluated using simple formulae.

The Schrodinger equation solved in Step (i) is generally different for each system size. The bound state nuclear motion problem has been extensively considered in the context of high resolution spectroscopy studies [27]. Introducing a finite region boundary into the problem does modify the Schrodinger equation somewhat and considerably change the nature of the solution, particularly in modified boundary conditions at $R_{AB}=a_{0}$ and through discretisation of the solutions above dissociation.

The Schrodinger equation solved in Step (iii) has the same general form for all system sizes (and indeed is the same as the equations used in R-matrix theory for electron-atom and electron-molecule collisions aside from a reduced mass factor); a single program can hence be used in this outer region, with only numerical considerations (e.g. step-size) changing between systems. However, the form of the reduced potential, $U(R_{AB})$, changes with system size and is often non-trivial, particularly if coordinate transformations are involved.

2.2 Developing the mathematical description of scattering processes using the RmatReact methodology

2.2.1 Establishing the problem

Much of the mathematics and analysis in this paper does not rely on the system being triatomic, and none relies on the system being \ceH3+ or its isotopologues, though this will be used extensively as a illustrative and useful example of the methodologies discussed.

Non-reactive elastic and inelastic scattering

A triatomic non-reactive scattering problem, e.g. \ceH+ + H2 -> H+ + H2, is most effectively solved using Jacobi coordinates, with $r_{1}$ as the diatomic bond distance, $r_{2}$ as the distance between the centre of mass of the diatomic and the scattering atom, and $\theta$ as the angle between these two vectors.

Photo-association and photo-dissociation

Non-reactive scattering coordinates are most appropriate here. These processes are half-collision processes involving scattering energy equal to the photon energy, $h\nu$, minus the difference in energy between the dissociation energy $D_{0}$ and the energy of the initial state, $E_{0}$, i.e. E = $h\nu-D_{0}+E_{0}$.

Charge Exchange

Treating charge exchange requires at least two electronic potential energy surfaces, called X and A here for simplicity. It is simplest to represent both potentials in a single coordinate system, the same as that used for non-reactive scattering.

Reactive scattering

Consider a reactive scattering system which includes scattering channels with products A+BC, B+AC and C+AB with coordinates A, B and C respectively. Jacobi coordinates can be defined for each scattering coordinate when considering the triatomic reactive scattering case. Simplification to just two channels (i.e. modelling just a single reaction) is straightforward, and the extension to more than three channels logical.

The mathematics described here builds on that presented in Chapter 7 of Burke [28]; this has been successfully used for the study of the positron-atom and positronium-ion reactive collision problem [29].

2.2.2 Nuclear Motion Schrodinger Equation: Solving the Inner Region Hamiltonian at energy-independent

The first stage of the RmatReact methodology is to find the wavefunctions and energies of the combined system in a finite inner region with zero scattering-energy. As these solutions will formally form a complete basis set in this finite region, these solutions can be used as a basis to describe the solutions to the scattering-energy-dependent problem in this inner region, i.e.

[TABLE]

where $E$ is the scattering energy, $\Psi(E)$ is the inner region solution to the scattering-energy-dependent Schrodinger equation, $i$ count the solutions to the energy-independent inner region problem $\psi_{i}$ and $A_{i}$ are expansion coefficients.

In traditional quantum chemistry treatments, the full Schrodinger equation can be simplified by ignoring translational wavefunction and separating electronic, vibrational and rotational wavefunction. The separation of the electronic component is an approximation, known as the Born-Oppenheimer approximation [30], while the separation identification of the vibrational and electronic components is not an approximation provided the Corolois term is included (as we do). The electronic component is considered in electronic structure packages to produce potential energy curves. In nuclear motion packages when treating a single electronic state, the total wavefunction can be represented as a sum of products between rotational and vibrational wavefunctions, i.e. $\psi_{\textrm{total}}=\sum_{i}\psi_{i}^{\textrm{rot}}\psi_{i}^{\textrm{vib}}$ (coefficients of the summation are absorbed into the vibrational wavefunction typically). The rotational wavefunction is a function of Euler angles $\alpha$, $\beta$, $\gamma$ and is quantised in terms of $J$ (the total angular momentum of the triatomic system), $M$ (the projected total angular momentum of the system onto the space-fixed $z$ axis), and $\Omega$ (the projected total angular momentum of the system onto the body-fixed $z$ axis) and described using Wigner $D$-functions [31], $D_{M\Omega}^{J*}$. Note in the absence of an external field, $M$ does not affect the energy of a molecular system and can be dropped from consideratio. Using this ansatz, the full Schrodinger equation is simplified to a set of $J$-dependent $3N-6$ dimensional Schrodinger equations that are typically given in internal vibrational coordinates. The total wavefunction thus becomes

[TABLE]

This basis set, or appropriately symmetrised versions of it, are used in variational nuclear motion programs, such as DVR3D [32] for triatomic systems, to yield the vibrational wavefunctions, $\psi_{J}^{\textrm{vib}}$. For triatomics, using Jacobi coordinates gives $\psi_{J}^{\textrm{vib}}(r_{1},r_{2},\theta)$.

The solutions to the traditional problem are bound state normalisable wavefunctions. In describing scattering, however, we need to include non-bound solutions corresponding to energies above the dissociation energy. Therefore, we move from an infinite region to a finite region, i.e. form a finite inner region, thereby discretising the continuum solutions. We ultimately desire the wavefunction solutions to the scattering Schrodinger equation at a large number of specific low scattering energies. An effective basis set to describe this large number of solutions can be formed by solving the single problem at zero scattering energy, as long as some of these solutions have non-zero value at the boundary between the inner and outer region (the R-matrix boundary). This is the first task of any R-matrix approach.

Defining the inner region is a key component of the RmatReact methodology; slightly different concerns are necessary for each type of scattering process.

Non-reactive elastic and inelastic scattering

We formulate the inner region by introducing a finite domain constraint in the scattering coordinate, i.e. $0\leq r_{2}\leq a_{0}$ (utilising Jacobi coordinates), where $a_{0}$ is known as the R-matrix boundary and defined as the scattering coordinate beyond which the two scattering systems interact via multipoles only to within the desired error, i.e. the potential beyond $a_{0}$ can be reduced from a full $(3N-6)$D potential to an effective 1D potential in the scattering coordinate $r_{2}$. In the inner region, we need at least one basis function which has a non-zero value at the R-matrix boundary $r_{2}=a_{0}$ in order to describe the scattering wavefunction. We have found Lobatto shape functions [33, 34, 35] to be a suitable choice of basis functions for this coordinate [36, 26].

Photo-association, photo-dissociation

The inner region for photo-association and photo-dissociation will be defined as for non-reactive scattering with the additional caveat that the treatment of these processes [37] becomes more complicated [38] if the quantum state of the combined system has significant magnitude beyond the R-matrix boundary (this case will not be considered here). In photo-dissociation, one would expect the reactant state to be well-bound and this will thus generally not be an issue. This is more likely to arise in photo-association when the product system (e.g. \ceH2D+) may be weakly bound in asymptotic halo states whose wavefunction is extended.

Charge-exchange reactions

The inner region R-matrix boundary $a_{0}$ should be defined such that it contains all non-negligible coupling between the two electronic states.

Reactive scattering

We need to expand this definition from a inner region defined by a single restricted domain to one that is constrained in all relevant scattering coordinates, i.e. $0\leq r_{2}^{A}\leq a_{0}^{A}$, $0\leq r_{2}^{B}\leq a_{0}^{B}$ and $0\leq r_{2}^{C}\leq a_{0}^{C}$. Again, the three $a_{0}$ values should be chosen such that the interaction between $A$, $B$ and $C$ can be modelled as a function of the scattering coordinate alone. To ensure the Hermiticity of the Schrodinger equation in this restricted region, we construct Bloch operators [39] at each of these boundaries like

[TABLE]

We then want to find the solutions to the Schrodinger equation with these Bloch terms, i.e. find $Z$ solutions $E_{i}$, $\Psi_{i}$ to

[TABLE]

where $T$ is the kinetic energy operator and $V$ is the potential energy operator for the system of interest. This is most easily done by using existing nuclear motion programs, e.g. DVR3D for triatomic systems, in which the Hamiltonian is modified to incorporate the Bloch terms and the basis functions are modified for this finite region. These inner region solutions $i$ will include both bound states of the molecule and a finite number of discretised continuum states. The strongly bound states will have energies and wavefunctions indistinguishable from the infinite region problem, but more weakly bound states will be influenced by the finite region constraints and be modified. As the box size increases, the differences for these weakly bound states will be smaller. The number of solutions will formally be equal to the number of basis functions utilised; however, we will only need to consider solutions with energies close to the scattering energy (often just above the dissociation energy of the reactant channel).

The Bloch terms and finite region introduces substantial requirements for the basis set, which must now ideally consist of basis functions with finite domain in three non-orthogonal scattering coordinates $r_{2}^{A}$, $r_{2}^{B}$ and $r_{2}^{C}$, with tractable resulting integrals. Furthermore, for each scattering boundary, at least one basis function must have a non-zero value but zero derivative; the zero derivative boundary condition introduced by necessity by the Bloch operator is a non-trivial and not often understood condition that has become obvious in our considerations of wavefunctions for RmatReact but which has been obscured in electron-molecule collision problems due to their far inferior basis sets: this fact is demonstrated in a simple system in the Appendix. These multiple boundary conditions are unusual and non-trivial constraints on the design of the basis set that are not yet fully understood, and will be the most challenging part of utilising the RmatReact methodology for reactive collisions. We will thus defer its consideration to Section 3. We should, however, be reassured by the fact that this type of approach has been successfully utilised in light particle reactive collisions, e.g. where positron-atom reactants react to positronium-ion products [28, 29].

2.2.3 Scattering Theory: Describing the Outer Region using Channels and Reduced Radial Functions

The close-coupling equations [40] also simplifies the full Schrodinger equation in order to progress; however, instead of separating based on rotational and vibrational wavefunction, the wavefunction in the scattering coordinate (the reduced radial function) is separated from the other components of the wavefunction, which are described as channel in standard close-coupling treatments.

Non-reactive scattering

The wavefunction in terms of channels and the reduced radial function, $F$ as

[TABLE]

where the definition of $F$ is by convention, $c$ goes over all channel functions and $\Phi_{c}$ are the channel functions with all necessary coordinates orthogonal to $r_{2}$. The product of the quantum states of the two isolated systems and their relative angular motion is typically used define the channels. For example, consider the triatomic system \ceA + BC. The scattering coordinate is $r_{2}^{A}$, the diatomic vibrational wavefunction is $\chi_{n}(r_{1})$, the rotational angular momentum of the diatomic is $j$ and the relative angular momentum between the atom and diatomic is $l$, with the channel defined by $nlj$ and having energy equal to the diatomic vibrational and rotational energy. Mathematically, we can use our knowledge of the energies and wavefunctions of the asymptotic reactant and products that define the channel based on solving the bound state $(N-1)$-dimensional problem for which there are well-developed program solutions in the nuclear motion community; in the case of atom-diatomic scattering, this means using codes such as level [41] or Duo [42]. The form Equation 5 can in principle be used to describe the wavefunction in any region of space. However, for small $r_{2}^{G}$, the number of channels will be high especially if the combined system is strongly bound (e.g. \ceH3+, \ceH2O), i.e. the potential energy surface is deep compared to the vibrational spacings. This is one key reason why the inner region problem is solved separately using traditional nuclear motion techniques in the full dimensional space rather than using the equations in this section. In the outer region, however, for non-reactive scattering problem, only a small number of channels are important for describing systems, especially for ultracold collsions where often only a single rotational levels is populated.

The Schrodinger equation to be solved in the outer region for non-reactive scattering is given by the close-coupling [40] expansion:

[TABLE]

where $\frac{\hbar^{2}}{2\mu_{2}}$ is the reduced mass along the Jacobu ‘scattering’ coordinate, $E$ is the scattering energy, $e_{nj}$ is the energy of the $n,j$ state of the diatomic and $U_{c,c^{\prime}}(r_{2})$ is the reduced potential $\braket{\phi_{c}}{\Delta V}{\phi_{c}^{\prime}}$ with $\Delta V$ equal to the difference between the total potential energy operator for the triatomic system and the diatomic potential energy operator; thus $U_{c,c^{\prime}}\rightarrow 0$ as $r_{2}$ increases. The values of $a_{0}$ should be determined such that $U_{c,c^{\prime}}(r_{2})$ is well represented by a multipole expansion in $r_{2}$.

Photo-association, photo-dissociation

In photo-association and photo-dissociation, only the reactant and product respectively consist of separated species (e.g. \ceH2 and \ceD+) that are described by channels; the combined system (e.g. \ceH2D+) will be described by the quantum numbers of the combined system.

Charge-exchange

The full wavefunction in the outer region will be described by channels in both electronic states, $c^{X}$, and $c_{A}$, i.e.

[TABLE]

The $c^{X}$ and $c_{A}$ channels are two sets of uncoupled channels. We can extend Equation 6 to sum over all $c=\{c^{X},c_{A}\}$ as long as $U_{c^{X},c_{A}}<<E$ for all $c^{X},c_{A}$, i.e. the interaction between channels on the different electronic surfaces are negligible compared to the scattering energy .

Reactive scattering

In this more complicated case, we need to be able to deal with multiple scattering coordinates (i.e. $r_{2}^{A},r_{2}^{B},r_{2}^{C}$), and thus channels in multiple coordinates, $\Phi_{c_{A}}$, $\Phi_{c_{B}}$ and $\Phi_{c_{C}}$, and multiple reduced radial functions $F_{c_{A}}(r_{2}^{A})$, $F_{c_{B}}(r_{2}^{B})$ and $F_{c_{C}}(r_{2}^{C})$. For generic coordinates, we will use ‘G’ as our notation. We can thus write the total wavefunction as

[TABLE]

where the coordinates of each of the channels are orthogonal to their associated $r_{2}^{C}$ coordinate.

The number of channels in the outer region is dependent on the difference in energy between the reactant and product; modelling more exothermic reactive scattering processes will necessitate a much larger number of channels with a consequent considerable increase in the calculation time.

To find $F_{c_{G}}(r_{2}^{G})$ beyond the R-matrix boundary $a_{0}^{G}$, we reduce the full-dimensional Schrodinger to three one-dimension scattering Schrodinger equation in each set of coordinates $G=A,B,C$ equal to the Equation 6. Note that by construction of the multiple R-matrix boundaries, the reduced potential $U$ connecting channels that scatter in different coordinates should be negligible compared to the scattering energy.

2.2.4 Energy-dependent Scattering using RmatReact methodology

Non-reactive scattering

The mathematics here has been discussed by Tennyson et. al. [24] and is summarised here to make clear the differences in treatment necessary between non-reactive scattering (the simplest kind of process) and the other processes, particularly reactive scattering.

We have discussed earlier the fact that we solve the inner region energy-independent problem in order to provide an efficient basis for describing the non-zero-scattering energy problem at the boundary, i.e. Equation 1. Mathematically, this approach utilises resolution of the identity (i.e. $1=\sum_{i}\ket{i}\bra{i}$) and provides a spectral representation of the Green’s function to find

[TABLE]

where $\mathcal{L}$ is the Bloch term given by $\mathcal{L}=\frac{\hbar^{2}}{2\mu_{2}}\delta(r_{2}-a_{0})\frac{d}{dr_{2}}$, $\Psi(E)$ is the scattering-energy-dependent wavefunction and $E$ is again the scattering energy. The summation runs over all inner region solutions, $i$, which is formally infinite but in practice finite due to the representation of the inner region solutions in a basis set of size $Z$. As this equation has the desired solution $\Psi(E)$ on both sides, this representation cannot be directly utilised in a computational solution.

Instead, mathematical transformations described in Burke [28] and Tennyson et. al. [24] yield an expression linking the reduced radial functions $F$ and the R-matrix $R$ at the R-matrix boundary $r_{2}=a_{0}$ as

[TABLE]

with the R-matrix defined by

[TABLE]

where the summation $i$ is over all inner region solutions, $c$, $c^{\prime}$ go over all channels and $\omega_{c,i}$ are the surface amplitudes defined by

[TABLE]

where the prime in the Braket notation indicates the integral is over all coordinates except $r_{2}$.

Note that the coordinate systems for the inner and outer region, and quantum numbers for the inner region solutions and outer region channels need to be carefully considered and can have a substantial effect on the calculation time and accuracy. In particular, we highlight that inner-region codes are generally based on body-fixed coordinates, whereas outer-region propagation codes will generally use space-fixed coordinates to allow identification of the asymptotic channels with the states of the fragmented systems. Thus, a frame transformation is required to convert between these; this is discussed in Appendix A (which also defines the below quantum numbers and notation). The key integral that must be calculated is the surface amplitude between inner region solution $i$ and outer region channel $njl$, denoted as $\omega_{njl,i}^{JM\epsilon}(a_{0})$, where the $JM\epsilon$ denote the quantum numbers for the combined triatomic system. As made clear by the notation, the outer region channel function is usually defined in space-fixed (SF) coordinates as $\phi_{njl}^{JM\epsilon}$. The inner region solution, ${}^{BF}\psi_{i}^{JM\epsilon}$, is defined in body-fixed coordinates such that it can be expanded as ${}^{BF}\psi_{i}^{JM\epsilon}=\sum_{n^{\prime}j^{\prime}\bar{\Omega^{\prime}}}F_{n^{\prime}j^{\prime}\bar{\Omega^{\prime}}}^{JM\epsilon}(r_{2})^{BF}\phi_{n^{\prime}j^{\prime}\bar{\Omega^{\prime}}}^{JM\epsilon}$. Thus the surface amplitude can be calculated using

[TABLE]

where we use the fact that the channels in body-fixed coordinates are defined such that they are orthonormal, i.e. $\Braket{{}^{BF}\phi_{n^{\prime}j^{\prime}\bar{\Omega}^{\prime}}^{JM\epsilon}}{{}^{BF}\phi_{nj\bar{\Omega}}^{JM\epsilon}}^{\prime}_{r_{2}=a_{0}}=\delta_{nn^{\prime}}\delta_{jj^{\prime}}\delta_{\bar{\Omega}\bar{\Omega^{\prime}}}$.

Photo-association and photo-dissociation

To describe photo-association and photo-dissociation processes using the R-matrix approach, we need to follow the approach to atomic photo-ionisation developed by Burke and Taylor [37] which links the initial wavefunction with the final wavefunction in the inner regions through the electric dipole.

In the case of photo-dissociation, let the initial wavefunction be $\Psi_{0}$, an eigenstate of the energy-independent inner region, i.e. $\psi_{i}$ for some $i$; this will often be the ground or low-lying rovibrational states of the combined system (e.g. \ceH3+). In the inner region, the final wavefunction in the inner region, $\Psi^{\textrm{inner}}(E)$, can be written as a sum of inner region wavefunctions $\psi_{f}$ where $f$ sums over all inner region states as in Equation 1, where the photon energy, $h\nu$, is equal to the scattering energy, $E$, plus the difference in energy between the dissociation energy $D_{0}$ and the energy of the initial state, $E_{0}$, i.e. $h\nu=E+D_{0}-E_{0}$. Then the integral that characterises the strength of the photodissociation process as a function of $h\nu$ is

[TABLE]

Then, we need to have the dipole moment function for the inner region so that $\Braket{\psi_{i}}{\mu}{\psi_{f}}$ can be calculated over the full finite inner region, as well as $A_{f}$, i.e. the coefficients of expansion for the scattering-energy-dependent inner region wavefunction in terms of the energy-independent wavefunction which are generated using non-reactive scattering with half-collision boundary conditions [37].

This approach should be capable of providing a complete model for the challenging H${}_{3}^{+}$ near-dissociation spectrum of Carrington and coworkers [10, 11, 12, 13, 14] including, for example, modeling the different resonance widths observed in their spectra.

The case of photo-association is analogous except that the initial and final states are switched, i.e.

[TABLE]

where this expression is only valid if the final wavefunction $\psi_{f}$ can be assumed to have negligible extent beyond the R-matrix boundary $a_{0}$. Within an R-matrix formulation it is also possible to consider the contribution of dipole transitions arising from the outer region [38], but this is beyond the scope of this paper.

Charge-exchange

Equation 9 applies for the charge exchange process with the channel index $c$ going over both the $c_{X}$ and $c_{A}$ channels.

Reactive scattering

Following the same logic as for the non-scattering case, we obtain

[TABLE]

with the key difference here being the need for multiple R-matrix boundaries and thus multiple Bloch terms. As we need to consider multiple boundaries and their associated coordinates, there is a substantially more complex form for the reduced radial functions $F$ and for the R-matrix, $R$. Following the methodology for a related derivation in Burke [28], we start from Equation 22, project onto channel functions $\phi_{c_{G}}$ (different for each coordinate), evaluate on the boundary, and ultimately obtain (in analogy with Eq. (7.26)–(7.27) of Burke [28]), yielding

[TABLE]

and the R-matrix found by

[TABLE]

where ${c_{G}}$ runs over all channels associated with coordinate $G$, i.e. $A,B,C$ if \ceA + BC, B + AC, C + AB are all considered explicitly, and where $\omega_{c,i}^{G}$ is the surface amplitude defined as

[TABLE]

Thus, with these equations, the R-matrix at a boundary $a_{0}^{G^{\prime}}$ can thus be obtained for any scattering energy if the inner region problem can be solved and if the surface amplitudes can be computed.

There are a number of important things to note about these expressions. First, the R-matrix is not symmetric with respect to its indices, i.e. the denominator is a function only of the second coordinate boundary, $a_{0}^{G^{\prime}}$. Second, the summation in this definition is formally over all inner region solutions to the energy-independent problem; however in practice bound state solutions corresponding to states well below dissociation will have near zero amplitudes at all three boundaries (i.e. $\omega_{ci}^{G}\approx 0,\forall G$) and can be excluded from the summation without errors. It should also be possible to significantly reduce the sum over inner-region solutions with positive scattering to within, say, two orders-of-magnitude of the scattering energy, $E$, using the so-called partitioned R-matrix approach [43, 44] which uses simple formulae based on perturbation theory to correct the R-matrix for the contributions due to higher energy poles. This ability to trim solutions should provide further contribute to the computational efficiency of this method. Finally, the key required integrals that must be evaluated to utilise these expressions are given by Equation 25, in which an inner region solution to the energy-independent problem $\psi_{i}$ is projected onto an outer region channel. How this is done must be considered carefully as it is this point in the calculation that one can introduce either a coordinate change and/or a frame transformation. Thus for example, one might wish to project solutions on outer region channels which have different non-orthogonal scattering coordinates. Issues arising from this and possible choices of inner-region basis sets are discussed in Section 3.

The issue of multiple coordinates to consider does substantially increase the difficulty of this problem. Note, for example, how the surface amplitudes for the inner region solutions need to be evaluated in all scattering coordinates, A, B, C, in order to evaluate $R$.

2.2.5 Propagation from R-matrix boundaries, and Asymptotic Expansion

There are standard procedures available for R-matrix propagation [45, 46, 47] which have been widely and successfully used for both light and heavy particle scattering. Our proposed solution is to use the parallel fast asymptotic R-matrix (PFARM) code [48]. Therefore considerations of how to propagate the R-matrix from the boundary to asymptotic distances, and to use asymptotic expansion and calculate scattering observables, is largely a solved problem. However, we make some notes in this section.

Non-reactive scattering

The propagation algorithm to go from the R-matrix boundary to its asymptotic value, and the way in which these asymptotic R-matrix is used to calculate the K-matrix and other scattering observables, is discussed for the single channel case in [36]. Extensions to the multi-channel case are in progress, based on the use of PFARM in the outer region.

Photo-association and photo-dissociation

Propagation proceeds in the same manner as for non-reactive scattering, except an explicit form for the wavefunction must be calculated (i.e. the $A_{i}$, $A_{f}$), and that half-boundary-conditions must be applied [37].

Charge exchange

Techniques for propagating the R-matrix from the R-matrix boundary to asymptotic regions for sets of uncoupled channels is developed for the case of two sets of channels in Appendix E6 of Burke [28], and can be adopted with minimal changes to this problem. The crucial thing here is that two different electronic states are assumed be non-interacting in the outer region; this assumption provides the criterion for choosing an appropriate R-matrix boundary. Note that the R-matrices associated with the channels associatd with the two different surfaces are not always zero, otherwise the rate of the charge exchange process would be zero. Therefore, the propagation of the entire global R-matrix containing the two set of uncoupled channels needs to be performed at one time [28].

The form of the K matrix expected from a charge exchange calculation is shown diagrammatically in Figure 1. Diagonal elements of this matrix represent elastic scattering processes, off-diagonal elements within a block represent inelastic processes and off-diagonal elements outside the two block diagonals represent charge transfer processes as the molecule moves from the X to A state.

Reactive scattering

Techniques for propagating the R-matrix and reduced radial wavefunction from the R-matrix boundary to asymptotic regions for sets of uncoupled channels is developed for the case of two sets of channels in Burke [28]; extensions to three sets of channels, if necessary for a particular problem, do not introduce fundamental changes to the approach here and will thus not be considered further in this paper. The mathematical framework used to solve the outer region problem of reactively scattering atoms from positrons to form ions and positronium (considered by Burke) is highly analogous to the framework needed to describe reactive scattering of atoms and diatomics. A full computational implementation of this mathematical framework will be the subject of a future paper.

The key required component for this outer region propagation is an expression of the reduced 1D potential energy curves, $U_{c,c^{\prime}}$, defined in Equation 6. A multipole expansion of the form

[TABLE]

is generally sufficient in the outer region, simplifying the propagation procedure. The value of the coefficients, $a_{c,c^{\prime},\lambda}$, can be obtained through appropriate integration of the non-diatomic potential, $\Delta V$, over the channels in all coordinates other than the scattering coordinate of interest, $r_{2}^{G}$. Thus for channels associated with H2 + H+\ceH3+. Note that the reduced potential between channels in different scattering coordinates must be negligible (and set as zero) for this methodology to be implementable.

The form of the K matrix for the full reactive scattering problem involving all three scattering coordinates will be a logical extension from Figure 1, with diagonal elements representing elastic scattering processes, off-diagonal in block elements representing inelastic processes and off-diagonal, off-block elements representing reactive scattering where the product and reactant are described using different Jacobi scattering coordinates.

3 Inner Region Coordinates and Basis Sets for Reactive Scattering

For polyatomic systems the variational nuclear motion programs we plan to use to provide solutions to the inner region problem, namely DVR3D [32], WAVR4 [49] and TROVE [50], all provide the choice of using a variety of different internal coordinates and some control over the basis set employed. So far we assumed the use of Jacobi coordinates and have not defined the basis sets will be used to compute the inner region energy-independent solutions $\psi_{i}$. Basis sets consist of basis functions that are generally products of one-dimension basis functions in each coordinate of interest. In practice, the nuclear motion problems are generally solved on a discrete variable representation (DVR) grid but transformation between polynomial basis functions and a DVR is straightforward [51], so it is sufficient at this stage to simply consider basis functions.

For studies of the bound states of \ceH3+ and its isotopologues including those up to dissociation, Jacobi coordinates have proven extremely successful [52, 53, 54, 55], despite not representing the full symmetry of the problem. In describing these processes in the inner region problem, we use Laguerre polynomials (either Morse-like oscillators [56] or spherical oscillator [57]) for the “diatomic”$r_{1}$ Jacobi coordinates, (associated) Legendre polynomials for the $\theta$ coordinate and Lobatto shape functions for the finite $r_{2}$ Jacobi coordinates. At the R-matrix boundary, we need at least one basis function to have a non-zero value (and formally the true solution has a zero-derivative boundary though this condition is less important); the need to satisfy this boundary condition and the finite domain in this coordinate is the main reason for utilising Lobatto shape functions.

For the reactive scattering problem a number of considerations need to be take into account meaning that there a several possible options for internal coordinates. Explicitly, we consider:

• •

A single set of Jacobi coordinates,

• •

Multiple sets of Jacobi coordinates,

• •

Hyperspherical coordinates,

• •

Radau coordinates.

Each of these coordinate choices leads naturally to a set of basis functions.

When assessing our options, we want to consider a few factors. Ideally, we would like to consider the ingong and outgoing channels on an equal footing. Second (and most importantly), we need to be able to computationally effectively evaluate the overlap integrals and Hamiltonian matrix elements arising from the basis set and coordinates in a doubly- or triply-finite region (depending on how many scattering channels are energetically accessible) as well as evaluating the surface amplitude integrals (i.e. Equation 25).

Single Jacobi coordinate basis set

For a triatomic reactive scattering code, it is potentially possible to use a traditional single set of Jacobi coordinate and define basis functions in this coordinate of the form

[TABLE]

for a set of $m,n,j,k$, where $M_{m}$ are Laguerre polynomials, $L_{n}$ are Lobatto basis sets and $P_{j,k}$ are associated Legendre polynomials.

For this type of single-coordinate basis set runs into the following problems:

The ingoing and outgoing channels are not treated equivalently; 2. 2.

The magnitude of the basis functions at the other R-matrix boundaries might not be sufficiently large to describe scattering properly in that coordinate; if the inner region basis functions cannot describe the wavefunction involved in scattering properly, then the RmatReact methodology will not be able to describe the scattering process properly; 3. 3.

The zero-derivative boundary conditions will not be met in the other coordinate systems (this is a desirable but in practice not necessary condition); 4. 4.

Accurately evaluating integrals with the finite boundary conditions imposed by the other scattering coordinate constraints is extremely complicated, especially as there will usually be no simple relationship between the coordinates.

Multiple Jacobi coordinate basis sets

To address some of the problems associated with using basis sets defined by a single set of Jacobi coordinates, we can introduce sets of basis functions associated with each Jacobi coordinate, i.e.

[TABLE]

for a set of $mA,nA,jA,kA,mB,nB,jB,kB,mC,nC,jC,kC$.

This multiple-coordinate basis sets approach alleviates the first two problems for a single Jacobi coordinate basis set, but the problem of efficient computation of integrals for a non-orthogonal and over-complete basis will be necessary for the methodology to be of practical usefulness; this is not a general feature of variational nuclear motion programs. A related approach was successfully utilised by Day and Truhlar [58], who used multiple Jacobi coordinate basis set to calculate the bound state energy levels of \ceH3+. This approach has the additional advantage of restoring the full symmetry to the treatment of the \ceH3+ problem.

The key difference between our problem and that solved by Day and Truhlar [58] is that our problem has multiple boundary conditions. For example, one of the simpler integrals required is of the form

[TABLE]

where $H(x)$ is the heavisidetheta function, i.e. $H(x)=1$ for $x>0$ and $H(x)=0$ for $x<0$. The best way to approach these integrals is to design the basis functions, $B$, such that $B(r_{1}^{A},r_{2}^{A},\theta^{A})\ll E$ for $r_{2}^{B}>a_{0}^{B}$. This is eminently feasible and will have the additional benefit of reducing linear dependency issues.

Hyperspherical coordinates

Hyperspherical coordinates represent one way to avoid the need for multiple internal region coordinate sets and thus non-orthogonal basis sets. Hyperspherical coordinates have already been used successfully for H+ + H2 reactive scattering [8]. However, hyperspherical coordinates become increasingly inefficient at large scattering distances and thus applications to these sorts of problems have generally require transformation into Jacobi coordinates at some large scattering distance for efficient treatment of the problem [9]. Furthermore, the inner region nuclear motion calculations using hyperspherical coordinates are much less efficient than DVR codes based on orthogonal Jacobi and Radau coordinates [59].

Radau coordinates

Problems where there are two reaction channels, such as \ceD + OH -> DO + H, can be efficiently treated using a single set of Radau coordinates (eg [60, 61, 62]). In this case, one would use Lobatto shape functions (or Radau shape functions [63]) for both the $r_{1}$ and $r_{2}$ coordinates (with the $\theta$ coordinate basis functions being Legendre polynomials as usual). The surface amplitude integrals will thus require projection of the basis functions defined in Radau coordinates onto the channel basis functions that are naturally defined in the Jacobi $r_{1}$ and $\theta$ coordinates; such transformations are relatively easily performed in a DVR representation (see Appendix A of Tennyson et al. [32]).

4 Conclusions

The RmatReact methodology is a new theoretical and computational methodology designed to treat ultracold heavy-particle scattering. The mathematical formalism borrows heavily from the highly successful calculable R-matrix methods that have been used extensively to treat electron-atom and electron-molecule collisions. However, there are crucial differences in the approach, particularly in the definition and solution for the inner region problem where the two scattering particles strongly interact. This paper extends for the first time the previously presented [24] mathematical framework for treating atom-atom inelastic and elastic scattering to all atom-diatomic scattering processes relevant for the \ceH3+ system: elastic and inelastic non-reactive scattering, photo-association and photo-dissociation, charge exchange and reactive scattering. The RmatReact methodology has the potential to revolutionise modelling of cold and ultracold heavy particle scattering by exploiting the inherent division of space into two regions: an inner region where the particle interactions are strong and should be treated in their full dimensionality with basis sets and calculation methods designed for molecular systems, i.e. nuclear motion methods, and the outer region where particle interaction is weak but must be considered to a very large interparticle distance due to the small collision energies.

Our detailed consideration of the mathematics required to describe all scattering processes using the new RmatReact methodology shows that the key difficulty is probably in the choice of coordinates and basis sets to describe reactive heavy-particle scattering problems like \ceD+ + H2 -> H+ + HD. This occurs because the inner region becomes finite as defined over every non-orthogonal scattering coordinate. The modifications to this inner region Schrodinger equation require that for each boundary at least one basis function must have a non-zero contribution, raising challenges in defining an appropriate basis set and evaluating the resultant integrals with multiple boundary conditions. To describe reactions with only two relevant reaction channels, e.g. \ceD + OH -> H + OD, Radau coordinates seem to be the most logical path forward. For systems where three reaction channels are of interest, e.g. symmetric H${}^{+}+\ce{H2}$ collisions, use of multiple Jacobi coordinates or hyperspherical coordinates appear to offer the best prospects for success.

Previous studies on the H${}_{3}^{+}$ system for reactive problems, eg D+ + H2($v=0,j=0$), have been performed on potentials with accurate long-range behaviour but a relatively poor representation of the H${}_{3}^{+}$ well region. We have recently developed a global H${}_{3}^{+}$ ground state potential energy surface [64] which joins the highly accurate ab initio spectroscopic potential of Pavanello et al. [65] with the correct treatment of the above dissociation and asymptotic regions due to Velilla et al. [66] to provide an accurate global potential for the \ceH3+ system. This will be used for our future studies on this system.

Acknowledgement

We thank Tom Rivlin and Eryn Spinlove for many helpful discussions during the course of this work.

\aucontribute

LKM carried out the main body of investigation and drafted the manuscript. JT conceived of the study and provided critical insight to develop and interpret the results. All authors read and approved the manuscript. \fundingThis project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 701962.

Appendix: Frame Transformation for Triatomic Systems

Channels are defined by three coordinates: the vibrational quantum number of the diatomic vibration $n$, the rotational quantum number of the diatomic $j$ and a relative angular momentum quantum number that is $l$ in space-fixed coordinates and $\bar{\Omega}$ in body-fixed coordinates. The inner region problem is generally solved in body-fixed coordinates with solutions labelled by $nj\bar{\Omega}$, while the outer region solutions utilise space-fixed coordinates $njl$. It is imperative to be able to convert between these representations. The frame transformation mathematics described here are adapted from Launay (1976) [67] to be appropriate to the present problem.

In both coordinate systems, the $r_{1}$ coordinate is the diatomic vector and the $r_{2}$ coordinate is the vector between the centre-of-mass of the diatomic and the scattering atom. For simplicity, the angular coordinates associated with each of these vectors are often referred to collectively as $\bar{\mathbf{r_{1}}}$ and $\bar{\mathbf{r_{2}}}$; they are defined differently for each coordinate system as described below.

The total angular momentum of the triatomic system is $J$, with a projection of $M$; each solution to the inner region problem, each channel, each full solution to the scattering problem and each reduced radial function solution to the scattering problem are labelled by these quantum numbers.

Body-fixed coordinates

In 2-angle embedding [68], the angular component of the wavefunction is given by

[TABLE]

where $Y_{lm}(\theta,\phi)$ is a spherical oscillator function and $\bar{\Omega}=0,1,....\min(j,J)$ if $\epsilon=(-1)^{J}$ and $\bar{\Omega}=1,....\min(j,J)$ if $\epsilon=(-1)^{J+1}$, while in 3-angle embedding, it is instead given by

[TABLE]

where $P$ is a Legendre polynomial.

The combined eigenfunctions of definite total parity are given by

[TABLE]

The expression for the channels in body-fixed coordinates, labeled by $nj\bar{\Omega}$, are thus

[TABLE]

An inner region solution $\alpha$ wavefunction in terms of channel functions is given by

[TABLE]

A scattering-energy solution $s$ at scattering-energy $E$ in terms of body-fixed channel functions is given by

[TABLE]

Space-fixed coordinates

The combined basis function for the angular coordinates takes the components of the angular momentum functions for the diatomic ($j,m_{j}$) and the atom rel. to the diatomic ($l,m_{l}$) that contribute to state with quantum numbers $J,M$ [67]:

[TABLE]

In space fixed coordinates, the channels are labelled by $njl$, and given by

[TABLE]

A scattering-energy solution $s$ at scattering-energy $E$ in terms of space-fixed channel functions is given by

[TABLE]

Connecting body-fixed and space-fixed coordinates

Define

[TABLE]

Then, since our functions are all real

[TABLE]

Since $\mathcal{Y}_{jl}^{\textrm{SF},JM\epsilon}(\hat{\mathbf{r_{2}}},\hat{\mathbf{r_{1}}})$ for all $l$ form a complete set of eigenfunctions in $\hat{\mathbf{r_{1}}},\hat{\mathbf{r_{2}}}$, we can use resolution of the identity to demonstrate that

[TABLE]

Similarly,

[TABLE]

Equating the body-fixed and space-fixed representations of the full wavefunction,

[TABLE]

Then, by comparing terms within the expression

[TABLE]

Similarly,

[TABLE]