Characterizing quasibound states and scattering resonances

Characterizing quasibound states from coupled-channel scattering calculations can be a laborious task, involving extensive manual iteration and ﬁtting. We present an automated procedure, based on the phase shift or S -matrix eigenphase sum, that reliably converges on a quasibound state (or scattering resonance) from some distance away. It may be used for both single-channel and multichannel scattering. It produces the energy and width of the state and the phase shift of the background scattering, and hence the lifetime of the state. It also allows extraction of partial widths for decay to individual open channels. We demonstrate the method on a very narrow state in the Van der Waals complex Ar–H 2 , which decays only by vibrational predissociation, and on near-threshold states of 85 Rb 2 , whose lifetime varies over four orders of magnitude as a function of magnetic ﬁeld.


I. INTRODUCTION
Scattering resonances are important in many areas of physics and chemistry. These include nuclear physics [1], electron scattering from atoms [2] and molecules [3], the spectroscopy of Van der Waals complexes [4], and chemical reaction dynamics [5]. They have manifestations in both spectroscopy and collisions, and describe both the decay properties of a quasibound state and the resonant scattering that occurs at energies close to the state.
Zero-energy Feshbach resonances are particularly important in ultracold atomic and molecular physics [6]. It is often possible to tune a near-threshold state across a scattering threshold with an external (usually magnetic) field. In this case, the resonance properties are usually observed as a function of field at essentially constant near-zero energy.
Magnetically tunable Feshbach resonances have been used both to control atomic interactions via the scattering length and to create diatomic molecules by magnetoassociation [7].
When there is a single open channel, the scattering phase shift δ(E ) increases by π across a resonance. Around an isolated narrow resonance, it follows the Breit-Wigner form as a function of energy, where E res is the resonance energy, is its width, and δ bg is the background phase shift. The lifetime of the corresponding quasibound state is τ =h/ . In multichannel scattering, the same behavior is shown by the eigenphase sum [4], which is the sum of the phases of the eigenvalues of the scattering S matrix, which all have magnitude 1. Individual S-matrix elements, however, show more complicated behavior. Scattering resonances are laborious to locate and characterize in scattering calculations. It is usually necessary to carry out calculations on a grid of energies across the resonance and fit the resulting phase shifts or eigenphase sums to Eq. (1) [4]. The width of the resonance is not usually evident until quite late in the process, so many manual iterations are often needed. In addition, the phase shift can usually be calculated only modulo π , and it is easy to miss a narrow resonance entirely if the grid used is too coarse. As an extreme example of this, some vibrational predissociation resonances in the Van der Waals complexes Ar-H 2 and Ar-D 2 are less than 10 −9 cm −1 wide [8], and occur at energies thousands of cm −1 above some of their dissociation channels. Locating such resonances in order to characterize them can be very challenging.
The purpose of this paper is to describe an automated procedure for converging on and characterizing resonances, using coupled-channel calculations of the phase shift or S matrix. Once the location is approximately known, the procedure can often characterize the resonance position, width and background phase shift with calculations at fewer than 10 energies, without the need for manual intervention. We will describe our method in terms of the phase shift, but it applies equally to the eigenphase sum.

II. THEORY
The Breit-Wigner functional form, Eq. (1), has three unknowns, and thus we need the value of δ(E ) at a minimum of three energies to characterize the resonance. We first define the dimensionless quantity a(E ) = tan δ(E ). This has a pole near the resonance position and Eq. (1) can be rewritten as where δ bg = arctan a bg , The quantities a bg , , and E res do not have immediate physical meaning, but they put Eq. (1) into a form (2) that is more convenient for numerical solution.
Our numerical procedure is to evaluate the phase shift (or eigenphase sum) δ(E ) from coupled-channel scattering calculations at three energies E 1 , E 2 , and E 3 to obtain the corresponding values a 1 , a 2 and a 3 . Solving three simultaneous equations allows us to extract a bg , , and E res . Defining we obtain and finally Equations (4) to (7) and (3) allow us to estimate the position, width and background phase shift of the resonance from three points in its vicinity. However, the estimates become numerically unstable if any two of the points are too close together. It is thus not satisfactory to converge on the resonance simply by generating a sequence of points that approach closer and closer to E res . Instead, we aim to generate a final set of three points, one near E res and two others at distances away from it comparable to the resonance width. In the present paper, we converge the central point upon E res . One of the outer points is placed about t lo from E res (with tolerance ±ξ t lo ), and the other is about t hi from it (with tolerance ±ξ t hi ). The three points can be regarded as allowing characterization of E res , , and δ bg , respectively. The logic we use to select which point to discard and where to place the next point is shown in Fig. 1. We terminate the iteration when the estimated value of E res is within a small amount of the closest of the three points and the other two points satisfy the criteria above. The algorithm shown in Fig. 1 is similar to the one we presented previously for converging on a zero-energy Feshbach resonance in the scattering length as a function of external field [9]. However, it allows better control of the placement of the final points and improves the sequence of points close to convergence. It is beneficial in combination with the procedures of Ref. [9] as well as for the present purpose.
The spacings t lo and t hi are signed, with |t lo | |t hi |. The values t lo = −0.1, t hi = 1.0 and ξ = 0.25 are usually appropriate for characterization of isolated resonances, and are used in the present work. However, for special purposes, it may be necessary to use different choices of the points. Larger values of t lo and t hi are sometimes needed for very narrow resonances, to reduce the effects of numerical noise, and smaller values may be needed for very wide resonances, to reduce variation in δ bg across the range. For overlapping resonances, it may be desirable to place both outer points on the same side of E res .
We need three energies in the vicinity of the resonance to start the procedure. We choose to use equally spaced points separated by a small amount δE , based on a physical estimate of the resonance width. The initial estimate of the resonance position could come from various sources, including approximate calculations or experimental measurements; we usually use the program BOUND [10,11], which solves the coupledchannel problem subject to bound-state boundary conditions.
The algorithms described here make the approximation that δ bg (E ) is constant across the range of points. This approximation improves as the convergence proceeds and the range of points becomes smaller. Nevertheless, it is the limiting factor that determines the distance from which convergence can be achieved. At least one of the initial points must give a phase shift that is affected by the resonance by more than the variation of δ bg (E ) across the range of the points. For very narrow resonances, computational noise in the phase shift can also affect convergence.

A. Vibrational predissociation of Ar-H 2
The first example we use to demonstrate this method is vibrational predissociation of the Van der Waals complex Ar-H 2 . The state with H 2 vibrational quantum number v = 1, rotation j = 0 and total angular momentum J = 0 lies about 4100 cm −1 above the ground state and can predissociate to form H 2 , v = 0, j 8. We use the interaction potential of Le Roy and Carley [12], and perform close-coupled scattering calculations using the MOLSCAT package [11,13] to evaluate the S matrix and its eigenphase sum as a function of energy. We use a space-fixed basis set that includes all functions with j 10 for v = 0 and with j 8 for v = 1 and solve the coupled equations using the symplectic log-derivative propagator of Manolopoulos and Gray [14] with the six-step fifth-order symplectic integrator of McLachlan and Atela [15].
This state has previously been characterized in coupledchannel calculations [8], but the purpose of this example is to show the efficiency of the present method. We locate the resonance that describes this state at E = 4139.075 cm −1 with respect to the threshold with v = 0, j = 0. We find a width = 2.02 × 10 −8 cm −1 , in agreement with the result reported in Ref. [8]. The decay is very slow because the interaction potential includes only low-order anisotropy and the nearest open channel, with j = 8, is only indirectly coupled to the bound state. Open channels with lower j are more directly coupled to the bound state, but decay to them releases more kinetic energy and this reduces their contributions. Table I shows the sequence of points E n generated by our procedure and how the successive estimates converge on the resonance parameters. Further details of the logic used at each step are given in the program output in Ref. [16]. To demonstrate the power of the method, we have purposefully picked an initial estimate that is further from the resonance than the best estimates from bound-state calculations. Our

013291-2
Input energies E1, E2, E3 and corresponding δ1, δ2, δ3 Evaluate Eres, Γ, δbg from Eqs. (4) to (7) and (3) Calculate dn = En−Eres |Γ| Sort 3 points and relabel min, mid, max such that Discard point corresponding to dmax  We have also characterized the same state using the improved interaction potential of Le Roy and Hutson [17]. In this case, we find the resonance at E = 4138.884 cm −1 , with width = 9.4 × 10 −10 cm −1 , which is a factor of about 20 narrower than for the potential of Le Roy and Carley [12]. The potential coupling terms responsible for vibrational predissociation are anisotropic terms off-diagonal in the H 2 vibrational quantum number, and these are dominated by the coefficient V λk (R) with λ = 2 and k = 1 in the potential expansion of Refs. [12,17]. Figure 3 of Ref. [17] shows that this coefficient is significantly weaker for the potential of Ref. [17] than for that of Ref. [12], and this is responsible for the decreased width.

B. Lifetimes of 85 Rb 2 Feshbach molecules
The second example is for 85 Rb 2 molecules produced by magnetoassociation at the ( f , m f ) = (2, −2) + (2, −2) threshold [18]. This is not the lowest threshold for 85 Rb atoms in a magnetic field, so the molecules can decay to lower thresholds in which one or both of the atoms has m f > −2. The Feshbach resonance used for magnetoassociation at a magnetic field near 155 G is caused by a state that has a lifetime around 82 μs at energies well below the (2, −2) + (2, −2) threshold. As it approaches threshold, the bound state acquires a large fraction of threshold character. The amplitude of its wavefunction at short range decreases in both the closed channel and the threshold channel. This reduces its width and increases its lifetime τ =h/ close to threshold.
We use the interaction potential of Strauss et al. [19], without retardation corrections. We perform coupled-channel calculations using the MOLSCAT [11,13] and BOUND [10,11] packages. The methods used are as described in Ref. [20]. The Hamiltonian includes atomic hyperfine coupling, electron and nuclear Zeeman terms, and coupling due to the singlet and triplet interaction potentials. The state decays through a combination of dipolar spin-spin and second-order spin-orbit interactions [21]. The wave function is expanded in a fully uncoupled basis set that contains all allowed spin functions for each value of the end-over-end angular momentum L of the colliding pair. We first perform calculations using a restricted basis set that includes only functions with L = 0. This removes channels with M F = m f ,a + m f ,b > −4, including the open channels, so that the quasibound states of interest become truly bound. We locate these bound states using the BOUND package, and use the results as initial estimates of the resonance positions. We then switch to a full basis set, including functions with L = 0 and 2, and carry out scattering calculations using the MOLSCAT package. Table II shows the convergence on the bound state below the (2, −2) + (2, −2) threshold at a magnetic field of B = 155 G, where it is bound by only 140 × h Hz and is less than 1 × h Hz wide. Under these conditions, the initial estimate of the binding energy from bound-state calculations is inaccurate by over 1900 times the width, which is about 13.5 times the binding energy. Nonetheless, our method converges with calculations at only 9 energies. Remarkably, even the estimate of the width from the first three points is within 20% of the converged value.
We repeat this convergence procedure for a range of magnetic fields and obtain the lifetimes shown in Fig. 2. At fields below 155 G, where the quasibound state is extremely narrow, the initial estimate of the energy from bound-state calculations is inadequate. In this region we obtain the initial estimate from linear extrapolation based on previous energies.
Köhler et al. [22] developed a model for the decay of this bound state. They obtained an analytic expression for the lifetime at magnetic fields just above the resonance, where the bound state is dominated by the open-channel component. Figure 2 includes the result of Eq. (9) of Ref. [22]; we use values for the real and imaginary parts of the scattering length from our coupled-channel calculations. The analytic expression of Ref. [22] is accurate close to the resonance, but breaks down far from it. The coupled-channel results correctly approach the lifetime of the underlying quasibound state, but the analytic expression for the lifetime drops unphysically to zero, because it is based on an expression for the closedchannel amplitude that is valid only close to resonance. Köhler et al. performed coupled-channel scattering calculations at nine specific values of the magnetic field. The calculations were too laborious for them to show the results as a curve rather than isolated points. It should be noted that they used an earlier interaction potential than the present work, which produced a significantly shorter lifetime (32 μs) for the underlying state. We have investigated the dependence of this lifetime on the dipolar spin-spin and second-order spinorbit interactions, which have opposite sign for the interaction potential of Ref. [19]. We find that the lifetime is strongly dependent on the degree of cancelation between the two.

C. Energy dependence of background scattering
At each step, the procedure described here obtains estimates of the resonance position and width. In their simplest form, these estimates neglect the energy dependence of the background phase shift or eigenphase sum δ bg (E ). They are in error if δ bg (E ) varies significantly across the current range of points. A background energy dependence may prevent identification of the resonance from far away and, in drastic cases, may prevent convergence on the resonance even from points near to it. However, in many cases, it is possible to obtain an estimate of δ bg (E ). It is then straightforward to redefine a(E ) = tan(δ(E ) − δ bg (E )) and use the rest of the procedure exactly as before.
The estimate of δ bg (E ) may come from a variety of sources. It may be as simple as estimating the background gradient from a plot of δ(E ), or it might include the effects of other resonances whose properties are already known. More sophisticated estimates may also be envisaged, based for example on semiclassical phase integrals for nonresonant open channels. In any case it is not necessary for the estimate to include the constant part of δ bg (E ), which is accounted for by the three-point approach, but only its energy dependence.
In the presence of uncorrected background variation in δ bg (E ), the resonance parameters obtained depend on the positions of the outer two points within the range permitted by the tolerance ξ . Different convergence sequences may result in differently positioned points and produce small apparently random variations in the parameters. Such variations may be unacceptable in some applications, such as calculating numerical derivatives of resonance parameters when performing least-squares fits of potential parameters to experimental data. They may be reduced or eliminated by choosing a small value of the tolerance ξ .

IV. PARTIAL WIDTHS
The product state distribution from decay of a quasibound state is characterized by a set of partial widths i for each open channel i. For an isolated narrow resonance, the partial widths sum to the total width . Across a resonance, each S-matrix element describes a circle in the complex plane [23,24] The partial widths are defined as real quantities, i = |g i | 2 , and the circles in the complex plane have radii √ i i / ; the complex coefficients g i also have phases φ i , which determine the orientations of the circles in the complex plane.
For each S-matrix element, Eq. (8) contains six independent parameters. We can determine the magnitude and phase of the product g i g i , but not g i and g i individually. As each element has real and imaginary parts, we need calculations at three points to determine these six parameters. However, the position E res and total width are common to the behavior of all S-matrix elements. Our usual procedure is therefore to keep these two parameters fixed at the values determined from fitting to the eigenphase sum, as described in Sec. II, and use calculations at only two points in order to characterize the real and imaginary parts of S bg,ii and g i g i . Equation (8) may be written where x = 2(E − E res )/ and D is a matrix with elements D ii = −2g i g i / . We insert calculated S matrices S 1 and S 2 at two energies E 1 and E 2 , and solve the resulting simultaneous equations to give The partial widths are obtained from the diagonal elements, We evaluate the partial widths after convergence on E res and as described in Sec. II. We use S matrices at two of the three energies used in the final iteration of the convergence procedure; these S matrices have already been calculated during the convergence, so we need no additional coupledchannel calculations. For narrow resonances (including the examples covered in this paper), where convergence is limited by numerical noise, it is best to use the points closest to and furthest from the resonance position. For wide resonances, where convergence is limited by a nonconstant background, it is best to use the two points closest to the resonance position. 013291-5 We have applied this procedure to the Ar-H 2 example described above, and the results are shown in Table III. It may be seen that they are in excellent agreement with those of Ref. [8], which were obtained by a much more laborious fitting procedure.
It is also possible to use a three-point approach based on the resonant circle in a single S-matrix element. This is analogous to the "fully complex" procedure of Ref. [9] and allows convergence on E res and as well as S bg,ii and D ii . In multichannel scattering, we have sometimes found that this procedure converges successfully, even for resonances where uncorrected energy dependence in the background prevents convergence based on the eigenphase sum. This occurs because the eigenphase sum contains background contributions from all open channels, whereas a diagonal S-matrix element contains a background contribution from only a single open channel. Nevertheless, it is usually more convenient to locate resonances initially in the eigenphase sum, using an estimated background energy dependence if necessary.

V. CONCLUSIONS
We have developed an automated procedure to converge on a scattering resonance and characterize it in terms of the resonance energy and width and the phase shift (or S-matrix eigenphase sum) of the background scattering. The procedure may be used for both single-channel and multichannel scattering. It eliminates the extensive manual iteration and fitting inherent in previous approaches.
The procedure requires an initial set of three energies in the vicinity of the resonance, one of which must be close enough to the resonance that its phase shift is affected by the resonance by more than the variation of the background phase across the range of the points. The strategy employed aims to generate a set of three points, one very close to the center of the resonance and two others at varying distances on either side. The three points allow characterization of the resonance position, width and background phase. Subsequent processing of the S-matrix elements allows extraction of partial widths for decay to individual open channels.
We have demonstrated the method on two very different test cases. One is a very narrow state in the Van der Waals complex Ar-H 2 , which lies far from any threshold and decays only by vibrational predissociation. The other is a nearthreshold state of 85 Rb 2 , whose lifetime varies from ∼80 μs to over 1 s as a function of magnetic field.
An important use of the procedure will be in least-squares fitting to determine interaction potentials from experimental results [25][26][27][28][29]. Such fits require numerical derivatives of calculated properties with respect to potential parameters, and the calculation of these derivatives is not feasible if manual intervention is required.
We have implemented the procedure in version 2020.0 of the general-purpose quantum scattering package MOLSCAT [11].