Low temperature scattering with the R-matrix method: argon-argon scattering

TL;DR

This paper introduces RmatReact, a new computational code for high-accuracy elastic atom-atom scattering calculations, applied here to argon-argon interactions relevant for ultracold chemistry, and compares results with existing data.

Contribution

The paper presents the first practical application of RmatReact for elastic scattering, demonstrating its effectiveness in calculating scattering observables for ultracold atoms.

Findings

Good agreement with previous results for scattering observables.

Resonances are slightly lower in energy and narrower than previous calculations.

Discrepancies in potential energy curves affect scattering length predictions.

Abstract

Results for elastic atom-atom scattering are obtained as a first practical application of RmatReact, a new code for generating high-accuracy scattering observables from potential energy curves. RmatReact has been created in response to new experimental methods which have paved the way for the routine production of ultracold atoms and molecules, and hence the experimental study of chemical reactions involving only a small number of partial waves. Elastic scattering between argon atoms is studied here. There is an unresolved discrepancy between different argon-argon potential energy curves which give different numbers of vibrational bound states and different scattering lengths for the argon-argon dimer. Depending on the number of bound states, the scattering length is either large and positive or large and negative. Scattering observables, specifically the scattering length, effective range, and partial and total cross-sections, are computed at low collision energies and compared to previous results. In general, good agreement is obtained, although our full scattering treatment yields resonances which are slightly lower in energy and narrower than previous determinations using the same potential energy curve.

Tables

Table 1. Table 1: The five PECs studied in this work, with their minima and equilibrium distances. The number of J = 0 𝐽 0 J=0 bound states derived in this work is the same as in all cited references (see Section III )). Note here the differing numbers of calculated bound states ( N bound subscript 𝑁 bound N_{\rm{bound}} ) between the methods.

Label	Authors	Citation	$N_{\mathrm{bound}}$	$V_{\min }$ / ${\mathrm{cm}}^{-1}$	$r_{\min }$ / Å
PM	Patkowski et al.	Patkowski et al. (2005)	9	-99.269	3.7673
Aziz	Aziz	Aziz (1993)	8	-99.554	3.7570
TT	Tang et al.	Tang and Toennies (2003)	8	-99.751	3.7565
PS	Patkowski et al.	Patkowski and Szalewicz (2010)	9	-99.351	3.7624
MD	Myatt et al.	Myatt et al. (2018)	8	-99.490	3.7660

Table 2. Table 2: Positions ( E res subscript 𝐸 res E_{\rm{res}} ) and widths ( Γ res subscript Γ res \Gamma_{\rm{res}} ) with standard errors of three quasibound states extracted from the supplementary data of Myatt et al. (MD) Myatt et al. ( 2018 ) , compared to the three shape resonances produced in this work by fitting eigenphases to Equation ( 12 ) (with background resonance paramters A 0 subscript 𝐴 0 A_{0} and A 1 subscript 𝐴 1 A_{1} ). The widths extracted from Myatt et al. Myatt et al. ( 2018 ) have been multiplied by two to match the convention employed in this paper.

$v$	$J$	$E_{\mathrm{res}}$ (MD) /${\mathrm{cm}}^{-1}$	${\Gamma }_{\mathrm{res}}$ (MD) /${\mathrm{cm}}^{-1}$	$E_{\mathrm{res}}$ (this work) /${\mathrm{cm}}^{-1}$	${\Gamma }_{\mathrm{res}}$ (this work) /${\mathrm{cm}}^{-1}$	$A_0$	$A_1$ /(${\mathrm{cm}}^{-1}$)-1
6	9	0.129	$0.660\times {10}^{-6}$	0.1287	$0.663\times {10}^{-6}$	0	0
6	10	0.448	0.00330	0.4486(2)	0.00247(46)	-0.107	0.757
7	5	0.071	0.00605	0.06993(5)	0.004841(5)	0.00213	-0.997

Table 3. Table 3: Scattering lengths ( A 𝐴 A ) and effective ranges ( r eff subscript 𝑟 eff r_{\rm{eff}} ) generated using four potentials compared to previous values. For the first three potentials, Aziz ( 1993 ); Patkowski et al. ( 2005 ); Tang and Toennies ( 2003 ) , the scattering lengths and effective ranges cited are from Barletta et al. Barletta et al. ( 2010 ) . For the fourth potential, the potential and scattering length are from the same source: Myatt et al. Myatt et al. ( 2018 ) .

Potential	$A$/Å (literature)	$A$/Å (this work)	$r_{\mathrm{eff}}$ /Å (literature)	$r_{\mathrm{eff}}$ /Å (this work)
Aziz Aziz (1993)	-505.6	-647.1	35.94	35.53
PM Patkowski et al. (2005)	1285	844.0	33.87	33.53
TT Tang and Toennies (2003)	-60.79	-62.50	50.12	49.20
MD Myatt et al. (2018)	-714	-709.3	–	35.41

Equations

$$R^J(k,a_0)=\frac{1}{a_0}\sum _{n=1}^N\frac{{\left(w_n^J(a_0)\right)}^2}{{\left(k_n^J\right)}^2-k^2},$$

$$k=\sqrt{\frac{2\mu E}{{\hslash }^2}},$$

$$\left({\hat{H}}^J+\mathcal{L}(a_0)\right)w_n^J(r)=E_n^J{w}_n^J(r),$$

$$\mathcal{L}(a_0)=\delta \left(r-a_0\right)\frac{d}{dr}.$$

$$F^{J}(k,a)=aR^{J}(k,a)\frac{dF^{J}(k,r)}{dr}\biggr{\rvert}_{r=a},$$

$$F_i^J(k,r)\underset{r\to \infty }{\sim }{s}_i^J(kr)+K^J(k)c_i^J(kr).$$

$$K^J(k)=\frac{R^J(k,a_p)k{a}_p{s}_J^{\prime }(k{a}_p)-s_J(k{a}_p)}{c_J(k{a}_p)-R^J(E,a_p)k{a}_p{c}_J^{\prime }(k{a}_p)},$$

$$s_J(kr)$$

$$c_J(kr)$$

$$\delta (k)=\arctan K(k).$$

$${\sigma }_{\mathrm{tot}}(k)=\sum _{J=J_{\min }}^{J_{\max }}\frac{4\pi }{k^2}(2J+1)\sin {({\delta }^J(k))}^2,$$

$$k\cot \delta (k)=-\frac{1}{A}+\frac{1}{2}{r}_{\mathrm{eff}}{k}^2,$$

$${\delta }^J(E)=A_0+A_{1}E+\arctan \frac{{\Gamma }_{\mathrm{res}}}{E-E_{\mathrm{res}}},$$

Full text

I Introduction

Laser cooling, Stark deceleration, buffer-gas cooling, and evaporative cooling are among a wide variety of cooling techniques developed in recent decades, which have allowed for precise control over individual molecules, especially diatomic molecules Stuhl et al. (2014, 2012); Shuman et al. (2010); Zhelyazkova et al. (2014). As such, a variety of experiments at millikelvin and microkelvin temperatures (so-called ‘ultracold’ temperatures) have now become routine, including experiments that probe collisions between small species in unprecedented levels of detail Molony et al. (2014); Dawid et al. (2018); Bohn (2001); Köhler et al. (2006); Heazlewood and Softley (2015); Bell and Softley (2009); Bell et al. (2009). This has led to the possibility of the fine-tuning of state-to-state reaction dynamics for reactions involving only a small number of partial waves Quemener and Julienne (2012); Ospelkaus et al. (2010). To quote Stuhl et al. Stuhl et al. (2014), this is “perhaps the most elementary study possible of scattering and reaction dynamics”. These experiments have led to the discovery of a variety of intriguing quantum phenomena, including shape resonances, Feshbach resonances Chin et al. (2010); Köhler et al. (2006); Pellegrini et al. (2008); Hutson et al. (2009), universal scaling laws D’Incao et al. (2004); Lompe et al. (2010), and Efimov trimers Lompe et al. (2010); Ferlaino et al. (2011); Wang and Esry (2009); D’Incao and Esry (2006).

Ultracold experiments have also revealed that, as with the well-known, near-dissociation $\rm{H}_{3}^{+}$ spectrum Carrington et al. (1982); Henderson and Tennyson (1996); Kemp et al. (2000); Tennyson et al. (2017), ultracold atomic collisions can have an overwhelmingly large density of resonances in scattering energy Mayle et al. (2013). Resonance states also offer the best opportunity for quantum control and steering: they are already being used to steer the formation of ultracold diatomic molecules Molony et al. (2014).

In response to these developments, this paper demonstrates a novel algorithm for the simulation of collisions between atoms, with the intention of extending the methods to collisions involving larger systems. This algorithm, known as RmatReact, is based on the computable R-matrix-based method widely applied to electron-atom and electron-molecule collisions Tennyson (2010); Burke (2011), which here has been adapted to the atom-atom case. The method is designed to study reactive and non-reactive, and elastic and inelastic collisions occurring over deep wells. With the exception of a single proof-of-principle study by Bocchetta and Gerratt Bocchetta and Gerratt (1985), this method has not been applied to so-called heavy particle scattering before.

The R-matrix method, being time-independent, is well-suited to studying the narrow, short-lived resonances considered here. In contrast to existing methods using only R-matrix propagation for heavy particle collisions Light and Walker (1976); Walker and Light (1980), the R-matrix method employed in this work makes full use of the partitioning of space into inner, outer, and asymptotic regions. This is in order to leverage the efficiency of variational nuclear motion programs at solving the short-range (inner region) problem and the R-matrix method for generating high resolution plots of scattering observables such as the cross-section. The R-matrix method is very similar in spirit to the multichannel quantum defect theory (MQDT) which has been extensively used to study ultracold atom-atom Burke et al. (1998); Raoult and Mies (2004); Gao et al. (2005a) and atom-molecule collisions Croft et al. (2011, 2012). Both methods consider the problem in two regions and the treatment of the outer region can be very similar. However, while MQDT approximates a full solution of the close-coupling equations by using quantum defects which only have a weak dependence on the collision energy, the R-matrix method aims to provide an exact solution to the close-coupling problem based on an inner region with no energy dependence.

The RmatReact algorithm developed in this paper has been discussed in two previous papers. In Tennyson et al. Tennyson et al. (2016) we presented a preliminary formalism, though the method has evolved since then. In Rivlin et al. Rivlin et al. (2018) we provided a demonstration of the method with comparisons to analytic Morse potentials.

In this paper, numerical results from this new algorithm are presented for the elastic scattering of argon atoms off other argon atoms at ultracold temperatures, ranging from sub-$\rm{\mu K}$ temperatures up to approximately 1 K ($=0.695$ $\rm{cm}^{-1}$, where $\rm{cm}^{-1}$ is used as a unit of energy). Much of the work in this paper is dedicated to re-creating existing results, and confirming known pieces of physics, in order to demonstrate the efficacy of the new methods developed as part of RmatReact.

Several different ground state $\rm{Ar}_{2}$ potential energy curves (PECs) Patkowski et al. (2005); Aziz (1993); Aziz and Chen (1977); Song et al. (2010); Tang and Toennies (2003); Patkowski and Szalewicz (2010); Myatt et al. (2018) are examined in order to simulate this scattering. These PECs are listed in Table 1. Despite having a shallow PEC formed from van der Waals forces, similar to other noble gas dimers, and in contrast to the more deep-well systems that the method was designed to study Tennyson et al. (2016), the argon-argon system has been chosen as a test system for the algorithm. This is in order to compare against existing experimental and computational results. The large number of high-accuracy PECs available for $\rm{Ar}_{2}$ make it a good candidate for testing the RmatReact method. Experiments have also been performed on cold ground state argon atoms Edmunds and Barker (2014).

Barletta et al. Barletta et al. (2010) studied low-energy Ar-Ar collisions in support of experimental studies using Ar for sympathetic cooling Barletta et al. (2009); they assessed four PECs for the $\rm{Ar}_{2}$ system. Of these four, three are also assessed in this work (PM, Aziz, TT, see Table 1); the fourth PEC of Slavíček et al. 03SlKaPa.Ar2 is not considered here, but two additional ones are. All five PECs studied here superficially appear very similar. However, as Table 1 shows, they do have slight differences which have significant impacts on their low-energy scattering properties. The PM and PS PECs were generated ab initio, whilst AZ, TT and MD used experimental results in their fit.

Barletta et al. Barletta et al. (2010) generated scattering lengths and effective ranges for the PECs they studied. These values, especially the scattering lengths, diverge significantly from each other. A further PEC, and the associated scattering length prediction (computed with the method of Meshkov et al. 11MeStLe.Rmat), from Myatt et al. Myatt et al. (2018), is also recreated here.

The issue of the highly varying scattering lengths appears to be closely linked to a long-standing debate over the number of vibrational ($J=0$) bound states belonging to the $\rm{Ar}_{2}$ system. Some PECs appear to support only eight bound states, while others appear to support a ninth bound state. If this state exists, it has a binding energy on the order of magnitude of 1 $\rm{\mu K}$, which is approximately 0.7 \rm{\mu}$$\rm{cm}^{-1}, or 86 picoelectronvolts, and thus would be difficult to detect. Nevertheless, the value of the scattering length of a particular system is highly dependent on the position of the highest bound state Burke (2011). Consequently, whether or not this state exists has important implications for the physics of the scattering.

Sahraeian et al. 18SaHa.Ar2 study two $\rm{Ar}_{2}$ PECs; those which we have labelled PM and PS. They claim to have detected the ninth bound state in both cases. This result for the PM PEC is in agreement with Barletta et al. Barletta et al. (2010).

The RmatReact method is described in Section II. The results presented in Section III include predictions of the scattering length and effective range, and partial and total cross-sections for a variety of partial waves, including the detection and characterisation of three ultracold shape resonances. These results are compared against literature results. In Section IV, along with the conclusions, some allusions to intended future works with this algorithm are presented.

II The RmatReact Method

R-matrix theory has existed in various forms since its invention by Eugene Wigner in the 1940s 46Wigner.Rmat; 47WiEi.Rmat. The underlying principle behind the R-matrix method is the partitioning of space into an inner region, an outer region, and an asymptotic region along the reaction coordinate $r$ Burke (2011). The radius of the boundary between the inner and outer regions is often designated $a_{0}$.

Since the reaction is assumed to be spherically symmetric, it can be modelled as taking place over one dimension, here represented by the internuclear distance, $r$. The angular dependence of the scattering observables is accounted for by splitting the overall three-dimensional scattering wavefunction into one-dimensional partial waves and summing over these waves. Each partial wave is labelled by a different value of $J$, which for the system studied here is the total angular momentum of the system.

In the inner region, the reactants are treated as a bound system. For two atoms scattering off each other, this means that the inner region consists of a one-dimensional diatomic PEC. The RmatReact method solves the one-dimensional, time-independent Schrödinger equation for this system over a range of $r$ values from a minimum of $r_{\rm{min}}$ to a maximum of $a_{0}$. Because the Schrödinger equation is being solved over a finite region instead of over all space, an extra surface term must be added to the equation to account for the surface term in the integration. This is known as the Bloch term 69Robson.Rmat. Note that this method differs from some R-matrix implementations, where a Buttle correction 67Buttle.Rmat; Burke (2011) is used to account for this issue.

In solving the time-independent Schrödinger equation with the Bloch term, the method diagonalises the inner region to produce finite-region rovibronic eigenenergies and eigenfunctions of the diatomic system, which are needed to construct the R-matrix on the boundary $a_{0}$ Burke (2011). The calculations in this region are independent of scattering energy, and so can be performed once for a given symmetry and for all scattering energies, hence greatly reducing the computational expense of the method.

Figure 1 shows the PEC of Myatt et al. (MD) Myatt et al. (2018). The eight bound states predicted by Myatt et al. are located below dissociation at their appropriate eigenenergies (with two states close to dissociation at $-1.58$ $\rm{cm}^{-1}$ and $-0.307$ $\rm{cm}^{-1}$). The lowest of the continuum states used by the RmatReact algorithm in the R-matrix calculation are also shown. At the scale shown in Figure 1, the MD PEC is not distinguishable from the other PECs studied in this work, and the eigenfunctions are very similar.

Slight numerical instability can be seen in the continuum states in the vicinity of $a_{0}$ in Figure 1. This is in part due to the smaller number of grid points used in the calculation to produce Figure 1 (as opposed to the larger number of grid points used in the results section), and is not sufficient to significantly impact the results. In practice, it is only the very last point in the grid which is important for the R-matrix calculations.

In the outer region, the reactants are treated as being unbound. However, they still interact over a long-range PEC. In this work only PECs which are polynomial in $r^{-1}$ at large values of $r$ are considered. The boundary between the outer and asymptotic regions is denoted here as $a_{p}$. In the outer region, the RmatReact method uses R-matrix propagation techniques, such as those due to Light and Walker Light and Walker (1976) or Baluja et al. 82BaBuMo.Rmat to extract the value of the R-matrix function at $a_{p}$ from its value at $a_{0}$. Although this propagation method does depend on scattering energy, it is considerably less expensive than the inner region calculation.

In the asymptotic region, the potential is assumed to be zero. Here, the scattering observables are calculated by determining the K-matrix at a given energy Burke (2011). Eventually the RmatReact method will utilise an asymptotic expansion, such as those developed by Burke and Schey 62BuSc.Rmat or Gailitis 76Ga.Rmat, for this calculation. However, at this point a simpler K-matrix formulation is used (as described below).

Low-energy resonances often have very narrow widths when plotted as a function of scattering energy. There are also often many resonances close together. As such, it is important to determine the scattering observables on a fine grid of energies. The R-matrix method’s inner-outer region separation is ideal for this task.

The details of the R-matrix method for the single-channel (elastic scattering) case have been discussed extensively in our previous RmatReact papers Tennyson et al. (2016); Rivlin et al. (2018), although there are differences from the version of the method in Tennyson et al. (2016). The following is an abridged explanation derived from Burke’s R-Matrix Theory of Atomic Collisions Burke (2011).

II.1 R-matrix theory

The R-matrix is a quantity with two equivalent definitions. In the single-channel case, the first definition is:

[TABLE]

where $k$ is the scattering wavenumber associated with the scattering energy $E$ via the equation

[TABLE]

$R^{J}(k,a_{0})$ is the R-matrix for a certain partial wave $J$, and $\mu$ is the reduced mass of the system.

$k_{n}^{J}$ are the wavenumbers associated with the rovibronic eigenenergies $E_{n}^{J}$ of the diatomic system in the inner region (labelled by quantum numbers $J$ and $n$, and following a similar relationship to Equation (2)), and $w_{n}^{J}(a_{0})$ are known as boundary amplitudes. $E_{n}^{J}$ and $w_{n}^{J}(a_{0})$ are the eigenvalues and eigenfunctions (evaluated at $a_{0}$) respectively of the time-independent Schrödinger equation with the Bloch term:

[TABLE]

where $\hat{H}^{J}$ is the Hamiltonian for the system for a given $J$, which includes kinetic and potential operator components, and the Bloch term, $\mathcal{L}(a_{0})$, is given by

[TABLE]

Burke Burke (2011) also provides a detailed explanation of how to derive the first definition of the R-matrix given in Equation (1) from the Schrödinger equation with the Bloch term.

The eigenenergies and eigenfunctions of Equation (3) are not restricted to bound states (in fact bound states tend to contribute little to the R-matrix sum). The numerical diagonalisation method used in the inner region creates a discretised continuum of $N-N_{\rm{bound}}$ above-dissociation states, see Figure 1, which all contribute to the R-matrix sum.

Note that this definition is entirely dependent on parameters which appear in the inner region and arise from the bound diatomic problem. In contrast, the second definition can be written as

[TABLE]

where $F^{J}(k,a)$ is the wavefunction associated with a particular partial wave $J$, which is evaluated at a particular point $a$. In the multichannel case, $F_{i}^{J}(k,a)$ is associated with a particular atomic channel $i$. Hence it is known as a channel function. $R^{J}(k,a)$ is the R-matrix, as in Equation (1). This definition is based on quantities that exist in the outer region. From this definition one can see that the R-matrix can be thought of as a form of ‘log-derivative’ of the channel function. This definition of the channel function is consistent with the definition used in MQDT 01Gao.Rmat; 05GaTiWi.Rmat, which also has similar definitions for the K-matrix and S-matrix (introduced below).

As a result of the equivalence of the two definitions of the R-matrix given by Equation (1) and Equation (5), information about the inner region bound problem can be used to obtain information about the scattering channel functions in the outer region. From these channel functions, scattering observables can be constructed via the K-matrix, $K^{J}(k)$, which is dependent on the asymptotic boundary condition involving the channel functions at arbitrarily large distances:

[TABLE]

Here $s_{i}^{J}(kr)$ and $c_{i}^{J}(kr)$ are ‘sine-like’ and ‘cosine-like’ functions which, in general, can have a variety of forms depending on the specific asymptotic region implementation of the R-matrix method being used (as described below).

II.2 Implementation

The RmatReact method is designed to act as a ‘harness’ between other codes that solve the inner and outer region problems. Ultimately it is intended for the harness to function with a variety of inner and outer region codes, with a “plug and play”mentality in mind.

In this work, the harness is only used with one inner region code: a modified version of the diatomic nuclear motion code Duo jt609. The version of Duo used here has been modified to use a discrete variable representation (DVR) basis 82LiPaLi.Rmat; 85LiHaLi.Rmat based on Lobatto shape functions, and to solve the inner region problem with the additional Bloch term. The Lobatto functions are derived from work by Manolopoulos, Wyatt, and others 93Manolopoulos.Rmat; 88MaWy.Rmat; 94Meyer.Rmat; 89MaWy.Rmat, which explain how to derive expressions for the kinetic and potential components of the Hamiltonian.

This is in contrast to the ‘sinc DVR’ method 92CoMi.Rmat currently implemented in Duo, which enforces a zero boundary condition on its eigenfunctions at the ends of the grid – clearly an unacceptable property for a method which relies on the amplitudes of eigenfunctions at the boundary. The Lobatto DVR method has boundary conditions that set the derivatives of the eigenfunctions at the boundary to zero, but the amplitudes themselves are allowed to take on non-zero values at $a_{0}$. This is also in contrast to the method of Bocchetta and Gerratt Bocchetta and Gerratt (1985), which used non-orthogonality and a grid that extended slightly beyond $a_{0}$ to produce arbitrary boundary conditions at $a_{0}$. This approach was tested in earlier versions of this work (see Tennyson et al. (2016)), but has since been supplanted by the Lobatto DVR methods.

The algorithm for generating Lobatto shape function nodes and weights is derived from Manolopoulos 93Manolopoulos.Rmat, with some modifications. The Duo code with Lobatto functionality used in this work is provided on the Duo GitHub page.

In this work, the outer region is handled in the harness code itself, with an iteration method in space based on the R-matrix propagation methods of Light and Walker Burke (2011); Light and Walker (1976); 76ZvLi.Rmat; 79ScWa.Rmat. In future work, however, this will be replaced with the fast R-matrix propagation code PFARM 02SuNoBu.Rmat, based on the sector diagonalisation method of Baluja et al. 82BaBuMo.Rmat, which also includes the asymptotic expansion of Gailitis 76Ga.Rmat. Preliminary testing with PFARM has demonstrated that it is able to re-create the resonances described in the results section of this paper. The specific implementation of the Light-Walker propagator used in this work can be seen in Equation (3) and Equation (4) of Rivlin et al. (2018).

The asymptotic region is addressed in this work using the following expression for the K-matrix:

[TABLE]

where $E$ is the scattering energy as before, and where $s_{J}(ka_{p})$ and $c_{J}(ka_{p})$ are given by:

[TABLE]

Here $j_{J}(kr)$ and $n_{J}(kr)$ are the spherical Bessel and Neumann functions respectively 14GrRy.Rmat, and the derivatives of $s_{J}(kr)$ and $c_{J}(kr)$ with respect to $r$, at the point $a_{p}$, are defined as $s_{J}^{\prime}(ka_{p})$ and $c_{J}^{\prime}(ka_{p})$. The distance, $a_{p}$, should be chosen such that the potential is sufficiently small by that point that the value of the K-matrix is not affected by the specific choice of $a_{p}$.

II.3 Scattering observables

The four quantities generated by the RmatReact method which are presented in this work are the eigenphase, cross-section, scattering length, and effective range of the argon-argon interaction. All of these can be constructed from the K-matrix, $K(k)$. In the single-channel case, they have simplified forms Burke (2011). The eigenphase, $\delta(k)$, sometimes known as the phase shift, is given by:

[TABLE]

As a result of this definition, the eigenphase (in radians) is the same modulo $\pi$. The eigenphases presented in this work are given in the range $[-\frac{\pi}{2},\frac{\pi}{2}]$. This leads to seeming discontinuities, e.g. in Figure 5 when the eigenphase passes through $|\frac{\pi}{2}|$. These discontinuities are characteristic of resonances which are also present in the eigenphases. Although the eigenphase is technically not a scattering observable itself, it can be used to construct the other three observables used in this work, and is consequently the best variable for the detection of resonances.

The total cross-section for an interaction can be given as the sum over partial waves from a minimum $J$ value $J_{\rm{min}}$ to a maximum $J_{\rm{max}}$, $\sigma_{\rm{tot}}(k)$. It is given by:

[TABLE]

and the cross-section for a given partial wave, $\sigma^{J}(k)$ is merely the summand of Equation (10).

The scattering length, $A$, and the effective range, $r_{\rm{eff}}$, are numbers that characterise the properties of the PEC and the scattering process at low energy. They can be defined in terms of a linear expansion at low energy. If one plots $k\cot{\delta(k)}$ for $J=0$ as a function of $k^{2}$, then for sufficiently low energy the plot should be linear. In this case, $A$ and $r_{\rm{eff}}$ are defined in the following way:

[TABLE]

ignoring higher order terms in $k^{2}$.

Another observable it is possible to detect using the eigenphase is a resonance. A resonance will appear as a feature in a plot of the eigenphase or cross-section as a function of $E$. Furthermore, the energy of the resonance, and its width – the inverse of its lifetime – can be determined by fitting a function to the eigenphase following the form of Breit and Wigner 36BrWi.Rmat; Burke (2011); jt31:

[TABLE]

where $\delta^{J}(E)$ is the eigenphase for partial wave $J$ at scattering energy $E$, $\Gamma_{\rm{res}}$ is the width of the given resonance, and $E_{\rm{res}}$ is the energy of the resonance. Note this definition of $\Gamma_{\rm{res}}$ follows that in standard use in scattering (eg jt31), and differs by factor of two from the definition a full width at half maximum (FWHM).

The non-resonant shape of the eigenphase (the ‘background’ eigenphase) is accounted for by the two terms $A_{0}$ and $A_{1}$, where it is assumed that the width is narrow enough that the background eigenphase can be approximated by a linear function of $E$ over its length. By fitting a generated eigenphase to a function of this form, values for $\Gamma_{\rm{res}}$ and $E_{\rm{res}}$ can be obtained. Note that it is sometimes necessary to replace the final term in Equation (12) with its negative, if it is required by the resonance shape.

III Results

III.1 Bound states

When performing the inner region calculations with Duo in this work, the number of bound states was found to be in agreement with literature values Barletta et al. (2010); 18SaHa.Ar2 for all five PECs studied (see Table 1).

However, there were considerable complications when attempting to detect the ninth bound state in this work for the PECs where it was predicted to exist – the PS and PM potentials. As this state is so weakly bound, it was necessary to extend the inner region calculations out to large distances in order to detect it. This ninth bound state has many similarities to a halo state 18OwSp.Rmat, as seen in Figure 2, which shows the ninth bound state as a function of $r$ for the PM PEC for when $a_{0}=50$ Å.

The results in Table 3 of this work were obtained only by extending $a_{0}$ out to distances of over 35 Å for the PM PEC and over 40 Å for the PS PEC. As such, a very large number of points needed to be used in order to maintain precision. The difficulty in detecting the ninth bound state is underlined by the fact that when the diatomic nuclear motion code LEVEL 16LeRoy.Rmat was used, the ninth bound state was never detected for any of PECs considered here, no matter how far out or how many points the inner region was integrated over. Sahraeian et al. 18SaHa.Ar2 also cited difficulties in detecting this state, which they quote a value of -0.86233\ \mu$$\rm{cm}^{-1} for.

Consequently the actual binding energy of the ninth bound state, for PECs in which it was detected, varied as a function of the $a_{0}$ used in the integration here, up to 100 Å. This is seen in Figure 3, which shows the value of the ninth bound state, $E_{9}$, of the PM and PS PECs as a function of $a_{0}$, for all values of $a_{0}$ under 105 Å for which the state was actually bound. If the calculations are converging on fixed values of $E_{9}$, they are significantly different from the value obtained by Sahraeian et al.

Once a threshold value of $a_{0}$ was reached, every PEC that the literature claimed had nine bound states were consistently found to do so, even if its value changed with $a_{0}$. No ninth bound state was detected in this work for any PEC for which it was claimed that there are only eight bound states, even when using large values of $a_{0}$ over 100 Å. More assessment of the numerical issues faced by the R-matrix method can be found in Rivlin et al. Rivlin et al. (2018).

We note that diffuse bound states whose wavefunctions extend significantly beyond $a_{0}$ can be found rather efficiently within an R-matrix formalism by performing scattering calculations with negative energies jt106; jt560. We plan to implement such a procedure within the RmatReact framework.

III.2 Resonances

The supplementary data provided by Myatt et al. Myatt et al. (2018) (MD) includes the rovibrational eigenenergies of the $\rm{Ar}_{2}$ system obtained using LEVEL 16LeRoy.Rmat. The supplementary data also quotes values for states which lie above the dissociation threshold but below a centrifugal barrier for $J>0$, known as quasibound states. The quasibound states from Myatt et al. Myatt et al. (2018) which have $J$ quantum numbers $J\leq 10$ are quoted in Table 2.

In this work, the quasibound states quoted for the MD potential in Myatt et al. Myatt et al. (2018) were characterised by analysing resonances in the scattering calculation. The diatomic nuclear motion code used in this work, Duo, does not have the capacity to detect quasibound eigenvalues directly (although it is possible to detect them using a stabilisation method with continuum states). However, these quasibound states should correspond to shape resonances, which can be detected in plots of the eigenphase and cross-section.

In order to detect the shape resonances, the RmatReact method was used to generate the eigenphase, and from it the partial cross-sections for all the partial waves with $J\leq 10$. The inner region was calculated using 500 Lobatto grid points between $r_{\rm{min}}=2.5$ Å and $a_{0}=22.5$ Å. The outer region propagation was performed from $a_{0}=22.5$ Å to $a_{p}=45$ Å, with over 1,000 propagation iterations.

Figure 5 shows the eigenphase and cross-section generated using the MD potential for $J=0$, $J=5$, and $J=10$. Figure 6 shows the eigenphase and cross-section generated using the MD potential for $J=9$. In all these cases, the eigenphase and cross-section were calculated for energies between $E=0.001$ $\rm{cm}^{-1}$ and $E=1$ $\rm{cm}^{-1}$.

The $J=0$ partial wave plots are included in Figure 5 to indicate what a typical eigenphase and cross-section looks like for this system when no resonances are present: in the $J=0$ cross-section plot the cross-section sharply rises at low energies.

Myatt et al. Myatt et al. (2018) predicted (see Table 2) that there should be quasibound states in the $J=5$, $J=9$, and $J=10$ partial waves. These resonances can clearly be seen in our calculations (Figures 5 and 6) where their positions are marked with dashed lines. These three states are the only quasibound states given by Myatt et al. for $J\leq 10$ and the only resonances detected in this work.

For the $J=5$ and $J=10$ resonances, the energy $E_{\rm{res}}$, width $\Gamma_{\rm{res}}$, and $A_{0}$ and $A_{1}$ parameters were fitted to the Breit-Wigner form of Equation (12), using the values quoted by Myatt et al. as the starting point of the fitting procedure. The very narrow $J=9$ resonance could not be fit in this way, and so the energy location of the width was determined by identifying where the eigenphase suddenly went from $\approx\frac{\pi}{2}$ to $\approx-\frac{\pi}{2}$ and identifying the two points either side of this jump; $E_{\rm{res}}$ was taken as the mid-point between them. This energy was then inserted directly into the Breit-Wigner fit.

Figure 4 shows the result of this procedure for the resonance in the $J=10$ partial wave. The fitting was performed using the energy range $E=0.4006$ $\rm{cm}^{-1}$ to $E=0.499501$ $\rm{cm}^{-1}$, using the Levenberg-Marquardt algorithm as implemented in the software Origin (OriginLab, Northhampton, MA).

Table 2 contains the results of this fitting procedure for all three resonances studied in this work (all using the same software and algorithm with appropriate energy ranges). The narrowest resonance is for $J=9$ and there is very good agreement between our results and those quoted by Myatt et al. Myatt et al. (2018). For the other two, broader resonances we find slightly different positions and widths. This is consistent with the full treatment of coupling to the continuum obtained in a scattering calculation: LEVEL, as used by Myatt et al. for their quasibound states, is known to be less well-adapted for characterising broader resonances 16LeRoy.Rmat; jt372. Both the resonance position and width for $J=10$ are also similar to the figures quoted by Čížek et al. 96CiHo.Ar2.

As Figure 6 and Table 2 show, narrow resonances can be hard to detect. The only resonances detected in this work were ones which had been previously predicted and only needed to be corroborated. In the future, a more sophisticated resonance-detecting software such as those by Tennyson and Noble jt31 or Noble et al. 93NoDoBu.Rmat, or possibly a procedure based on the complex analysis of the S-matrix Burke (2011) such as that of Čížek and Horáček 96CiHo.Rmat, will be used to to detect resonances which may otherwise be missed.

Finally, Figure 7 shows the total cross-section generated using the RmatReact method with the MD potential. The quasibound states predicted by Myatt et al. Myatt et al. (2018) are also pictured. This figure gives a good overview of the properties of argon-argon scattering at low energy. It is notable for having many features. Besides the three resonances, there is also more structure to the plot – something that is more prevalent in heavy particle scattering than electron-atom or electron-molecule scattering due to the greater number of partial waves contributing to the scattering process. Furthermore, the cross-section tends to a large value at the lowest energies on the graph. This corroborates the feature seen in Figure 8 towards the lowest energies where the cross section becomes very large.

Thus far we have not considered the consequences of the Pauli principle. 40Ar is a Boson with zero nuclear spin; as a consequence collisions with odd $J$ are forbidden. Figure 7 shows the observable cross section obtained by simply summing partial waves with even $J$. As a consequence the resonances with $J=5$ and $J=9$ disappear and there is a pronounced Ramsauer minimum at about 0.01 $\rm{cm}^{-1}$.

III.3 Low-energy Scattering

In order to analyse low-energy scattering behaviour, the cross-section for $J=0$ was plotted for $E=10^{-8}$ $\rm{cm}^{-1}$ to E = 1 $\rm{cm}^{-1}$ on a log-log axis, see Figure 8. The same numerical parameters were used as in Section III.2. The plot shows that the cross-section tends towards a constant at lower energies, which is predicted by Equation (11).

Figure 9 analyses the region of validity of the low-energy linear fit of Equation (11). It is designed to re-create a plot shared as private communications by the authors of Ref. Barletta et al. (2010).

The solid, red line of Figure 9 represents the eigenphase calculation generated by the RmatReact method, using an R-matrix inner region ranging from $r_{\rm{min}}=2.5$ Å to $a_{0}=82.5$ Å, an integration over 1600 Lobatto grid points, and an R-matrix propagation from $a_{0}=82.5$ Å to $a_{p}=165$ Å, with 1000 propagation iterations. The dashed line represents Equation (11), with the parameters $A$ and $r_{\rm{eff}}$ determined by using a least-squares linear fit of the lower-energy portion of the red line (intercept $=0.00146$ Å*-1*, slope $=18.42$ Å), again using the software Origin. As with Table 3, this Figure is in agreement with results provided in private communications by Barletta et al. Barletta et al. (2010), who also computed the scattering length of the $\rm{Ar}_{2}$ collision based on the potential due to Aziz Aziz (1993). It can be seen from Figure 9 that the plot is only linear at a very low energy.

A similar low-energy fitting procedure in Origin was performed for all five PECs studied. The values of $A$ and $r_{\rm{eff}}$ were calculated in this work for the four PECs where corresponding literature values could be obtained, the comparison of which can be seen in Table 3.

The effective ranges featured in Table 3 all appear to be in broad agreement. This is to be expected since this quantity is not especially sensitive to fine changes to the quantity of the potential, and is not affected significantly by the number of bound states Burke (2011).

The values obtained for the scattering length are found to be sensitive to the energy range used in the fitting procedure, and so whilst numbers are quoted in Table 3, it should be noted that these numbers are not intended to be definitive. When using the energy range of Figure 9 for the low-energy fit, it is possible to obtain the scattering lengths quoted in Barletta et al. (2010) to within a 5% relative difference. However, when using a much lower energy range for the fit of $k^{2}\approx 10^{-10}$ Å*-2* to $k^{2}\approx 10^{-8}$ Å*-2*, the scattering lengths change significantly. (The effective ranges also change slightly, but are still in agreement.) The values quoted in Table 3 are the ones created using the lower energy range fit. As Figure 9 shows, this lower range is where the expansion of Equation(11) is most appropriate.

The features seen towards the right of Figure 8 correspond to energies where the eigenphase pass through zero. On a log-log plot of the cross-section these crossings manifest as the dips seen in the Figure.

Although the scattering length values diverge from each other very significantly, the RmatReact method was able to qualitatively corroborate each one. The PECs in Table 3 which have a negative scattering length correspond to PECs for which there are eight bound states in literature (see Table 1). The only PEC considered which supports nine bound states, PM Patkowski et al. (2005), has a large, positive scattering length.

This is in line with the observation that the scattering length is strongly affected by the energy of the highest bound state. If the scattering length is plotted as a function of $V_{\rm{min}}$, the minimum of the potential, then there is a pole at points where the number of bound states increments by one, going up to positive infinity in one direction and down to negative infinity in the other Burke (2011). That means that either side of this pole, the scattering length can be very different: any real number is a potentially valid scattering length.

It is known 00FaVaSn.Ar2; Patkowski et al. (2005) that relativistic and nonadiabatic effects can impact potential parameters such as the depth of the potential. The different PECs studied in this work all incorporate these effects to different degrees. Whilst this work attempts to verify the scattering observables produced using these potentials, no attempt is made to assess the quality of each potential relative to the other ones. These effects, along with the other sources of uncertainty related to the PECs, are by far the biggest source of uncertainty and error in the results, and contribute much larger error to the numbers quoted here than numerical errors in the algorithm itself.

No previous values are available for the scattering length and effective range of the PS PEC Patkowski and Szalewicz (2010); the scattering length and effective range were calculated, using the same lower energy range fitting as the results in Table 3. The scattering length was found to be 1669 Å to four significant figures. This is noteworthy because both this work and Sahraeian et al. 18SaHa.Ar2 claim to have detected nine bound states for this system, and so the PS PEC continues the pattern of large, positive scattering lengths for $\rm{Ar}_{2}$ PECs with nine bound states, as seen in Table 3. Finally, the effective range was found to be 33.82 Å, in good agreement with most of the other effective ranges cited in the literature and this work.

IV Conclusions and outlook

In this paper the validity and accuracy of the RmatReact method for the single-channel, diatomic case was demonstrated by comparing results generated using it to other literature results. In doing so, the accuracy of the scattering length, and the positions of the resonances generated by Myatt et al. Myatt et al. (2018) were confirmed. Most of the widths of the resonances generated by Myatt et al. were also confirmed.

This paper corroborated the qualitative features of the highly divergent scattering lengths quoted in Barletta et al. Barletta et al. (2010). This has interesting implications for the study of the low-energy behaviour of the argon-argon scattering interaction – the debate over the scattering length remains unresolved. Novel experimental techniques such as those in 18BeMe.Rmat may help to resolve the dispute over the $\rm{Ar}_{2}$ scattering length and the alleged ninth bound state.

Further study of the single-channel, atom-atom scattering problem is intended. A resonance finder will be useful for detecting any narrow resonances missed by other authors. The S-matrix can be used for this purpose, and also for the equally useful purpose of detecting weakly-bound bound states Burke (2011); jt106.

In resolving the numerical difficulties of adapting pre-existing codes to the ‘harness’ of the RmatReact method, this work paves the way for the study of more complex interactions with the method. Other follow-ups to this work will include a study of a multichannel collision between atoms, and collisions between an atom and a diatom.

Eventually the RmatReact method is intended to evolve into a method that can be applied to even more complex reactants and reactions, to resolve the many outstanding questions in the field of ultracold scattering. A formulation of the method for treating chemical reactions in three particle systems has recently been presented jt758.

Acknowledgements

This project has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie grant agreement No 701962, and from the EPSRC under grants EP/M507970/1 and EP/R029342/1. We would also like to acknowledge the contributions of the late Professor Robert J. Le Roy, in this and many other works we have been involved in.