---
abstract: |
    In this contribution, we introduce numerical continuation methods and bifurcation theory, techniques which find their roots in the study of dynamical systems, to the problem of tracing the parameter dependence of bound and resonant states of the quantum mechanical Schrödinger equation. We extend previous work on the subject [@Broeckhove2009] to systems of coupled equations.

    Bound and resonant states of the Schrödinger equation can be determined through the poles of the $S$-matrix, a quantity that can be derived from the asymptotic form of the wave function. We introduce a regularization procedure that essentially transforms the $S$-matrix into its inverse and improves its smoothness properties, thus making it amenable to numerical continuation. This allows us to automate the process of tracking bound and resonant states when parameters in the Schrödinger equation are varied. We have applied this approach to a number of model problems with satisfying results.
address: |
    Department of Mathematics and Computer Science, Universiteit Antwerpen,\
    Middelheimlaan 1, B-2020 Antwerpen
author:
- 'Przemysław Kłosiewicz , Jan Broeckhove and Wim Vanroose'
bibliography:
- 'coupled.bib'
title: Numerical Continuation of resonances and bound states in coupled channel Schrödinger equations
---

Introduction {#sec:Introduction}
============

The appearance of resonances is of ever-growing interest in the study of wave phenomena as they are considered among the most important features of systems described by wave equations. They appear in systems that are penetrable by an impacting wave. Such systems allow the interior field to couple to the external domain which leaves a characteristic fingerprint on the far-field pattern of the scattered wave. Many examples appear naturally in acoustic scattering [@fahy2007sound] and in fluid-mechanical structure interaction [@dhia2007resonances]. In all cases the appearance of resonant states has a profound and important influence on the system’s dynamics.

In the context of quantum mechanical scattering, resonant behavior also strongly influences the interactions between microscopic particles, which in turn has its influence on the reactivity of molecules and atoms described by such quantum mechanical models. In molecular systems these resonances can easily turn into bound states if the molecular configuration changes.

In [@Broeckhove2009] we developed a framework for applying numerical continuation techniques in the context of bound states and resonances in spherically symmetric short-range potentials. We have shown that numerical continuation methods, originally developed in the study of dynamical systems, can be applied successfully to track bound and resonant states efficiently in terms of a varying system parameter. Moreover, this technique can be used to reveal subtle and interesting transitions and connections between states automatically.

The present work focuses on the extension of that procedure to coupled channel short-range systems. This extension is a logical step towards automated, efficient and robust methods for the study of interactions in scattering experiments. In all generality, these techniques can be applied to systems of coupled Helmholtz equations with variable wave numbers, as long as the short-range conditions are met.

The outline of the paper is as follows. Section \[sec:Scattering\] sets the coupled channel Schrödinger equations in the context of non spherically symmetric quantum mechanical problems. In section \[sec:Regularization\] a regularization procedure is discussed that allows application of numerical continuation even though the underlying functions that characterize resonances and bound states are numerically and analytically not very well-suited. Section \[sec:NumCont\] provides a brief overview of basic numerical continuation methods and gives some pointers on the available implementations. Finally, in section \[sec:results\] we present several results obtained with our implementation of the discussed methods.

Quantum scattering in coupled channel problems {#sec:Scattering}
==============================================

In this paper we discuss a coupled channel problem that derives from a one particle Schrödinger equation with a non-spherical potential $$\label{eq:3dschrodinger}
    \left(-\frac{1}{2\mu}{\Delta} + V(\mathbold{r},\lambda) - E  \right) \psi(\mathbold{r}) =0,$$ where $\Delta$ is the three-dimensional Laplacian, $V(\mathbold{r},\lambda)$ is a potential with a system parameter $\lambda$ and $E$ is the complex-valued energy of the system. The problem is such that for all $\mathbold{r}$ outside a bounded domain $V(\mathbold{r},\lambda) \approx 0$, i.e. the potential becomes negligibly small. Formally, the limitation to potentials that are negligible outside a certain radius is termed as the restriction to so called short-range potentials: $V(\mathbold{r},\lambda)$ must decay faster than $r^{-3}$ as $r=|\mathbold{r}|\to\infty$ and must be less singular than $r^{-2}$ in the origin $r=|\mathbold{r}|\to0$ [@Taylor2006]. The long-range Coulomb interaction requires a significantly different approach and is not discussed here. The boundary conditions that are appropriate in the short-range case require the solution to be zero at the origin and force the solution to tend to a linear combination of free waves (i.e. solutions of for $V=0$) at infinity [@Taylor2006].

In this section we briefly introduce elements from quantum scattering and partial wave expansion to arrive at the concept of resonant states. In the subsequent sections we then focus on applying numerical continuation to study the dependence of such resonances on the system parameter $\lambda$ in . For this reason, we will faithfully record the $\lambda$ dependence in our notations, even if it is at times somewhat cumbersome.

Equation can be written in spherical coordinates $(r,\theta, \varphi)$ around the center of the system. The differential operator $\Delta$ then splits into angular and radial differential operators and the solution can be expanded as a sum $$\psi(\mathbold{r}) = \psi(r,\theta,\varphi) = \sum_{l=0}^{\infty}\sum_{m=-l}^{l} \psi_{lm}(r) Y_{l}^{m}(\theta, \varphi),
    \label{eq:PartialWave}$$ where $Y_{l}^{m}(\theta,\varphi)$ are the spherical harmonics, eigenfunctions of the angular differential operators in the Laplacian. In physics this decomposition is commonly referred to as the partial wave expansion. The labels $(l,m)$ are intimately connected to the irreducible representations of the rotational symmetry groups $SO(3) \supset SO(2)$ of equation . Each $(l,m)$-component of is known as a partial wave and the system is said to be modeled by multiple partial wave channels. We are now interested in localizing the resonances in such systems.

Single-channel scattering {#sphericalsymmetry}
-------------------------

To introduce some notations and concepts we first briefly look at a spherically symmetric potential, i.e. $V(\mathbold{r},\lambda) = V(r,\lambda)$. In this case, the partial wave channels are decoupled and the radial wave function $\psi_{lm}$ is identical for all $m$. Hence we can drop the index $m$ and equation turns into: $$\label{eq:equationspherical}
    \left(-\frac{1}{2\mu} \frac{d^2}{d r^2} + \frac{l(l+1)}{2\mu r^2} - E \right) \psi_{l}(r,\lambda) + V(r,\lambda) \psi_{l}(r, \lambda) = 0,$$ for each $l$. The boundary condition at $r=0$ specifies $\psi_{l}(r,\lambda) = 0$ and at large $r$ the solution must be a linear combination of the spherical Riccati-Hankel functions (see appendix \[app:spherical\]), the free incoming and outgoing waves: $$\label{eq:sphericalboundary}
    \psi_{l}(r, \lambda) \xrightarrow{r\to\infty} \frac{i}{2} \left( \hat{h}^{-}_{l}(kr) - \hat{h}^{+}_{l}(kr) S_{l}(k, \lambda) \right),$$ The spherical Riccati-Hankel functions are the solutions obtained in the absence of the potential term in , i.e. when the equation reduces to a Helmholtz equation with a constant wave number $k=\sqrt{2\mu E}$. The $S_{l}(k, \lambda)$ is called the $S$-matrix for channel $l$ and it determines the scattering properties associated with potential $V(\mathbold{r},\lambda)$.

Multi-channel scattering
------------------------

For a non-spherical potential the partial wave channels do not decouple. The wave function must be represented by a sum as in equation , which is typically truncated at an $l_\text{max}$. In most cases of physical interest the potential still has axial symmetry. It is well known that as a consequence the channels with different $m$ are decoupled and the solutions can be represented by $$\psi(\mathbold{r}) = \sum_{l=0}^{l_\text{max}} \psi_{lm}(r) Y_{l}^{m}(\theta, \varphi).
    \label{eq:PartialWaveAxial0}$$ Upon substitution in , this generates a separate set of equations for each m. For the purpose of our exposition we may, without loss of generality, take $m=0$ and drop the index $m$ in the notation for the radial wave function $$\psi(\mathbold{r}) = \sum_{l=0}^{l_\text{max}} \psi_{l}(r) Y_{l}^{0}(\theta, \varphi).
    \label{eq:PartialWaveAxial}$$ This significantly simplifies notations further on in the discussion. Lifting the restriction of axial symmetry presents no fundamental difficulty: the number of channels increases and summations over $m$ have to be reintroduced. We will however consider only the axially symmetric case with $m=0$ as it provides a sufficient basis to demonstrate that numerical continuation can be used to track resonant states in multi-channel systems.

Insertion of in , followed by a projection onto each of the partial wave channels (by multiplying with $Y_{l'}^{0}$, integrating over the angular variables and using the orthonormality of the spherical harmonics) leads to a set of coupled radial equations $$\label{eq:coupledchannel}
    \left(-\frac{1}{2\mu} \frac{d^2}{d r^2} + \frac{l(l+1)}{2\mu r^2} - E \right) \psi_{l}(r, \lambda) \ + \ \sum_{l'=0}^{l_\text{max}} V_{ll'}(r,\lambda) \psi_{l'}(r, \lambda) \ = \ 0,$$ for each channel $l$ and where $$V_{ll'}(r,\lambda) = \int \sin\theta \, \text{d}\theta \, \text{d}\varphi \, {Y_{l}^{0}}^{*}(\theta, \varphi) V(r,\theta,\lambda) Y_{l'}^{0}(\theta,\varphi),$$ and where $^{*}$ denotes the complex conjugate. For every partial wave $\psi_{l}$ the boundary condition at $r=0$ is still $\psi_{l}(r, \lambda) = 0$ while for large $r$ one formulates $l_{\text{max}}+1$ independent boundary conditions (one for each $l'$) $$\label{eq:incomingpartial}
    \psi^{l'}_{l}(r, \lambda) \xrightarrow{r\to\infty} \frac{i}{2} \left( \hat{h}_{l}^{-}(kr) \delta_{ll'} - \hat{h}_{l}^{+}(kr) S_{ll'}(k, \lambda) \right),$$ where $\delta_{ll'}$ is the usual Kronecker symbol.

Each boundary condition corresponds to an incoming wave in a different channel $l'$ and gives rise to an independent associated solution with its channel components labeled by $l$. A key difference with the spherical case is that, due to the non-spherical action of the potential, outgoing wave components appear in channels other than the incoming channel, as indicated by the $S$-matrix defined by its elements $S_{ll'}(k, \lambda)$.

We now rewrite the above in terms of matrices by introducing the square matrices $\mathbold{T}$, $\mathbold{V}$ and the identity matrix $\mathbold{I}$ with dimensions equal to the number of channels: $$\left(\mathbold{T}(r)\right)_{ll'} =  \delta_{ll'}\: \left(-\frac{1}{2\mu} \frac{d^2}{d r^2} + \frac{l(l+1)}{2\mu r^2} \right),\qquad \left(\mathbold{V}(r,\lambda)\right)_{ll'} = V_{ll'}(r,\lambda).$$ Eq.  then is equivalent to $$\left( \mathbold{T}(r) - E\mathbold{I} + \mathbold{V}(r, \lambda) \right) \mathbold{\Psi}^{l'}(r, \lambda) = 0,$$ where the wave function is a vector function with dimension equal to the number of channels, corresponding to a solution with (easily rewitten in vector format). In multi-channel calculations it is however customary to introduce the square matrix $\mathbold{\Psi}$ whose columns are the $\mathbold{\Psi}^{l'} $ (corresponding to ) for each of the $l'$. This then leads to the set of coupled channel equations and boundary conditions that determine the full $S$-matrix which is the focus of our study $$\label{eq:matrixeq}
    \left( \mathbold{T}(r) - E\mathbold{I} + \mathbold{V}(r, \lambda)  \right) \mathbold{\Psi}(r, \lambda) = 0,$$ $$\label{eq:matrixboundary1}
    \mathbold{\Psi}(0, \lambda) = 0,$$ $$\label{eq:matrixboundary2}
    \mathbold{\Psi}(r, \lambda) \xrightarrow{r\to\infty} \frac{i}{2}  \left( \hat{\mathbold{h}}^-(kr) - \hat{\mathbold{h}}^+(kr) \mathbold{S}(k,\lambda) \right).$$ The matrices $\hat{\mathbold{h}}^{\pm}(kr)$ are diagonal matrices with the incoming and outgoing spherical Riccati-Hankel functions on the diagonal $$\left(\hat{\mathbold{h}}^{-}(kr)\right)_{ll'} =  \delta_{ll'} \hat{h}_{l}^{-}(kr), \qquad \left(\hat{\mathbold{h}}^{+}(kr)\right)_{ll'} = \delta_{ll'} \hat{h}_{l}^{+}(kr).$$ The multi-channel $S$-matrix $\mathbold{S}$ is the object of interest that we will use to find the resonant states. We can extract the $S$-matrix from the numerical solution of with its boundary conditions by evaluating the Wronskian in a point $r$ outside the range of the potential $$\label{eq:Smatrix_and_wronskian}
    \mathbold{S}(k,\lambda) = \mathcal{W}\left[\mathbold{\Psi}(r,\lambda),\hat{\mathbold{h}}^{+}(kr)\right]^{-1}\,\, \mathcal{W}\left[\mathbold{\Psi}(r,\lambda),\hat{\mathbold{h}}^{-}(kr)\right],$$ where $\mathcal{W}$ denotes the Wronskian of two matrices defined as $$\label{eq:Wronskian}
    \mathcal{W}[\mathbold{A},\mathbold{B}] = \mathbold{A} \left(\frac{d \mathbold{B}}{dr} \right)-\left(\frac{d \mathbold{A}}{dr} \right)\mathbold{B},$$ with $\mathbold{A}$, $\mathbold{B} \in \mathbb{C}^{n\times n}$.

Resonances and regularization {#sec:Regularization}
=============================

The aim of this article is to develop a robust method to trace the parameter dependence of resonances that appear in multichannel problems. We define a resonance or bound state to be the solution of equation for an energy $E\in\mathbb{C}$ for which the $S$-matrix has a pole. This definition is widely used and it is motivated, amongst others, in [@Newton1982; @Taylor2006].

It is natural to locate the resonances numerically by searching for the zeros of $\mathbold{S}(k,\lambda)^{-1}$. However, it is not the best numerical strategy to search for these zeros directly. A key property of the $S$-matrix is the symmetry $\mathbold{S}(k^{*},\lambda) = \mathbold{S}(k,\lambda)^{-1}$, where $k^{*}$ denotes the complex conjugate of $k$. This makes $\mathbold{S}(k,\lambda)$ unitary for $k\in\mathbb{R}$. As a result a pole of the $S$-matrix at $k$ is always accompanied by a zero of $\mathbold{S}(k,\lambda)$ at $k^*$. This makes the search for the zeros of $\mathbold{S}(k,\lambda)^{-1}$ a numerically hard problem, particularly when $k \ll 1$. In that case $\mathbold{S}(k,\lambda)$ has zeros and poles close to each other and it does not satisfy the smoothness conditions required to ensure rapid convergence if, for instance, Newton iterations are used to find the zeros.

For single-channel radial problems, as in section \[sphericalsymmetry\], it was proposed in [@Broeckhove2009] to use a regularized function rather than the $S$-matrix itself to converge and to track down bound states and resonances. This function was defined as, per $l$-channel, $$\label{eq:one_channel_regularized}
    F_{l}(k,\lambda) = \frac{k^{2l+1}}{S_l(k,\lambda)-1},$$ with $S_{l}$ as in . It was shown that for $k\neq 0 $ this function has a zero iff $S_l(k,\lambda)$ has a pole. In addition this function is bounded as long as $\int_0^\infty dr\, \hat{j}_l(kr) V(r,\lambda) \psi_l(k,r) \neq 0$, which does not hold only for very specific potentials.

In this paper we extend this result to the more general class of coupled channel problems described in the previous section. We introduce, in extension of , the $(l_\text{max} + 1)$ by $(l_\text{max} + 1)$ square matrix function $$\label{eq:multi_channel_regularized}
    \mathbold{F}(k,\lambda) = \mathcal{K}
    \left( \mathbold{S}(k,\lambda) - \mathbold{I} \right)^{-1},$$ where $\mathcal{K}$ is a diagonal matrix with elements $(\mathcal{K})_{ii} = k^{2l_{i}+1}$ on the diagonal and where we extract the $S$-matrix using from the numerical solution of , and . The function used for continuation then is $$\label{eq:det_multi_channel_regularized}
    \det(\mathbold{F}(k,\lambda)) = \left( \prod_{l=0}^{l_{\text{max}}} k^{2l+1} \right) \bigg/ \det(\mathbold{S}(k,\lambda)-\mathbold{I}).$$

In appendix \[app:regularization\] we show that $\det(\mathbold{F}(k,\lambda))$ has a zero if and only if $\det(\mathbold{S}(k,\lambda)^{-1})$ has a zero. This means that we can look for the bound states and resonances via and investigate the application of continuation to equation $\det(\mathbold{F}(k,\lambda))=0$ . In addition, the coalescing of zeros and poles that occurs in the $S$-matrix is now absent. Unfortunately, we cannot derive a bound for $\det(\mathbold{F}(k,\lambda))$ itself as was possible for . See appendix \[app:regularization\] for details.

Numerical Continuation methods {#sec:NumCont}
==============================

In the study of dynamical systems one is often interested in the dependence and evolution of solutions of nonlinear problems in terms of some system parameter. Generally these solution sets can only be approximated numerically and finding those requires efficient and robust techniques. More specifically, numerical continuation techniques are concerned with approximating the solution set $\left\{ \mathbold{x} | F(\mathbold{x}) = \mathbold{0} \right\}$ of an underdetermined system of nonlinear equations: $$\begin{aligned}
    F : \mathbb{R}^{n+1} &\longrightarrow \mathbb{R}^{n} \\
    \mathbold{x} = (\mathbold{u},\lambda) &\longmapsto F(\mathbold{u},\lambda),\end{aligned}$$ that depends on some system parameter $\lambda$. Such equations arise for instance in the study of dynamical systems or parametrized systems of ODEs but virtually every problem that can be written in the above formulation and that satisfies certain smoothness conditions can be subjected to those methods.

Due to the computationally intensive nature of many of these problems, efficiency is one of the main concerns in the numerical study of such dynamical systems. The number of evaluations of the function $F(\mathbold{x})$ should therefore be kept to a minimum. Additionally, solution sets often exhibit complex geometries with intersecting and bifurcating solution branches. Generally, the study of bifurcating branches involves rigorous stability analysis of the underlying solutions and is a complicated subject on its own. In this work the only type of bifurcations that can occur, are so called “simple” bifurcation points. These manifest themselves as two intersecting solution branches and are characterized by the dimension of the null-space of the Jacobian $F_{\mathbold{x}}(\mathbold{x}_{t})$ being 2 in a point $\mathbold{x}_t$, which is called a “branch point”. In many problems, however, the bifurcations are much more involved and a thorough treatise on bifurcations and stability of solutions can be found in [@Seydel1994; @Doedel2007; @Allgower1990; @Keller1977].

In the last few decades various techniques for numerical continuation have been developed. Many of those fall in the category of so called “predictor-corrector” methods: starting from a known solution point $\mathbold{x}_{i}$ (such that $F(\mathbold{x}_{i})=\mathbold{0}$), the next one is found by first taking a guess based on the previous point and correcting that guess in a second step. As such, starting from one known initial point $\mathbold{x}_{0}$ in the solution set an approximation $\left\{ \mathbold{x}_{i} | i=0,\ldots,n \text{ with } F(\mathbold{x}_{i})=\mathbold{0} \right\}$ of that solution set is constructed by finding subsequent points on the curve.

Many choices can be made for the implementation of both steps. A good overview can be found in [@Seydel1994; @Doedel2007; @Allgower1990; @Keller1977]. The popular “pseudo arclength” method [@Keller1969] takes a predictor-corrector approach. The predictor step uses the tangent $\mathbold{t_{i}}$ to the solution curve at $\mathbold{x}_{i}$ to construct a prediction $\mathbold{x}_{i+1}^{p}$ for the next point $\mathbold{x_{i+1}}$. The corrector step consists of Newton iterations on the system of equations $F = \mathbold{0}$ augmented with an additional equation that constrains the corrector step to a plane orthogonal to the tangent $\mathbold{t_{i}}$. The Newton iterations then converge to the intersection of that plane and the solution curve, guaranteeing a unique next point $\mathbold{x_{i+1}}$ in the continuation [@Keller1969].

The pseudo arclength algorithm has been implemented in the AUTO [@Doedel1981; @AUTO2007] library which we have used in this work. Other well-known implementations of numerical continuation algorithms include LOCA [@Salinger2002] which is a part of the Trilinos framework [@Heroux2005], <span style="font-variant:small-caps;">Matcont</span> [@Dhooge2006] (a Matlab implementation) and Multifario [@Henderson2002] which allows multi-parameter continuation.

Numerical Results {#sec:results}
=================

The results of this section have been obtained with a renormalized Numerov solver [@Johnson1977; @Johnson1978] for the calculation of the coupled channel solutions of equation . However, these results are independent of the particular numerical method used to solve the equations as long as the method is able to extract the $S$-matrix e.g. through equation . Indeed, the regularized function $\det(\mathbold{F}(k,\lambda))$ is the only required input for the continuation routines.

To demonstrate our approach we have applied numerical continuation to a number of model problems involving well-behaved short-range potentials. The numerical implementation of the continuation algorithm itself is provided by the AUTO-07p Fortran library [@AUTO2007]. The coupled channel renormalized Numerov solver and the code to obtain the regularized function $\det(\mathbold{F}(k,\lambda))$ are used as driver routines for AUTO and we have implemented them in C++.

Comparison with analytically solvable models {#sec:model1}
--------------------------------------------

As a first testing example for the study of numerical code and regularization we have compared our results thoroughly with a model problem from [@Vanroose2002]. There, among other examples, a short-range two-channel coupling model for the low-energy electron-CO$_{2}$ scattering is studied. An analytically solvable expression with square-well potentials is presented together with resulting bound/resonant pole trajectories.

The application of numerical continuation on both the analytical expressions and the numerical solutions of this model matched the results previously published by the authors.

Gauss potential with $s$-wave and $p$-wave coupling {#sec:model2}
---------------------------------------------------

In this section we look at a system of two coupled Gaussian potential wells with angular momenta $l_{1}=0$ (referred to as $s$-wave in physics) and $l_{2}=1$ (referred to as $p$-wave). The mass is chosen to be $\mu=1$. The two channel potentials and the coupling potential are set to: $$V_{ij}(r,\lambda) = -\lambda_{ij}e^{-r^{2}},$$ where $i,j\in\{1,2\}$ and $V_{ij}(r,\lambda) = V_{ji}(r,\lambda)$.

Let us first consider each of the two channels separately, i.e. uncoupled ($\lambda_{12}=0$). Setting $\lambda_{11}$, the strength of the first-channel potential, to 7 gives a system with one bound state shown in table \[tab:gauss\_bs\]. Similarly, bound states in the second channel with $\lambda_{22}=20$ are shown in table \[tab:gauss\_bs\].

These results were calculated through Newton iterations on the function $\det(\mathbold{F}(k,\lambda))$. They compare, see table \[tab:gauss\_bs\], favorably with results obtained with <span style="font-variant:small-caps;">Matslise</span> [@Ledoux2005].

    $i$   $k_{i}$ (<span style="font-variant:small-caps;">Newton</span>)   $k_{i}$ (<span style="font-variant:small-caps;">Matscs</span>)   $k_{i}$ (<span style="font-variant:small-caps;">fin. diff.</span>)
  ----- ---------------------------------------------------------------- ---------------------------------------------------------------- --------------------------------------------------------------------
      0                                                  `3.623677e+00i`                                                  `3.623551e+00i`                                                      `3.623677e+00i`
      1                                                  `2.178012e+00i`                                                  `2.178012e+00i`                                                      `2.195645e+00i`
      2                                                  `9.035406e-01i`                                                  `9.035275e-01i`                                                      `9.039594e-01i`

  : Coupled channel bound states in the $s$- and $p$-wave example, with Newton iteration (<span style="font-variant:small-caps;">Newton</span>), coupled channel extension of CPM{$P$,$N$} (<span style="font-variant:small-caps;">Matscs</span>) and finite differences (<span style="font-variant:small-caps;">fin. diff.</span>)[]{data-label="tab:gauss_l0_l1_c05_bs"}

The coupled system on the other hand has bound states at positions deviating from those in table \[tab:gauss\_bs\] as the coupling strength increases. In this case, with $\lambda_{12} = 0.5$, our results, see table \[tab:gauss\_l0\_l1\_c05\_bs\], are in close agreement with those obtained through the recently extended CPM{$P$,$N$} method for coupled channel systems [@Ledoux2006; @Ledoux2007] implemented in <span style="font-variant:small-caps;">Matscs</span> [@Matscs]. Reference values obtained with eigenvalue calculations based on finite difference discretization, although not very accurate, are also provided.

It is now interesting to see how these bound states evolve in terms of variations in $\lambda_{22}$. Applying numerical continuation with starting values $k_{i}$ (<span style="font-variant:small-caps;">Newton</span>) taken from table \[tab:gauss\_l0\_l1\_c05\_bs\] gives results shown in figure \[fig:gauss\_l0\_l1\_c05\_continuation\]. The continuation was performed both in the direction of the origin and away to show as much of the resonant trajectories as possible.

Details of the underlying numerical procedures are as follows. The integration interval of the coupled channel Numerov method spans from 0 to 4.6 and contains 4096 points. The tolerance for Newton iterations used by AUTO is of the order $10^{-6}$ whereas the length of the continuation steps might vary adaptively between $10^{-2}$ and $10^{-4}$.

Each of the three continuation branches originates from a bound state and follows a path towards the resonant region (i.e. the lower half of the complex $k$-plane) as the potential strength $\lambda_{22}$ is decreased. Near the origin of the complex plane the second channel potential becomes too weak to hold the bound state and reaches a threshold value where the state transforms into a resonance. Looking at the continuation curves in figure \[fig:gauss\_l0\_l1\_c05\_continuation\] we can therefore distinguish three regions: the bound state regime on the positive imaginary axis, the transition into a resonance around the origin and the resonant regime below the real axis. Their properties are summarized as follows. (1) The bound states are situated on the positive imaginary axis and correspond to real, negative (rel. to potential asymptote) energies. In the neighborhood of the bound states of the first channel, in this case $k\approx 2.2i$, “avoided crossings” can be seen where a relatively big decrease in $\lambda_{22}$ has a minor influence on $k_{i}$. However, as soon as the lower state $k_{i-1}$ approaches, the higher one gets pushed away from this position and makes place for $k_{i-1}$. (2) As $\lambda_{22}$ passes certain threshold values, the potential becomes too weak to hold a bound state $k_{i}$ which transitions into a resonance. For $s$-waves these transitions happen on the negative imaginary axis not far from the origin, whereas for higher angular momenta this branching occurs exactly at 0. In the case of $s$/$p$-wave coupling, both effects influence each other depending on the coupling strength and bound states of the separate channels. The threshold values for $\lambda_{22}$ with the corresponding branch points for the three states in our example are summarized in table \[tab:branchpoints\]. (3) In the resonant regime ($\text{Im}(k)<0$) the states gain a real part which corresponds to a positive scattering energy. The negative imaginary part of these states is related to the time delay in the time-dependent picture. These branches are symmetric but only those in the fourth quadrant have physical significance. On the negative imaginary axis one also finds so called virtual states.

    $i$   $\text{Re}(k_{i})$   $\text{Im}(k_{i})$   $\lambda_{22}^{t}$
  ----- -------------------- -------------------- --------------------
    $1$      `-1.227997e-12`      `-7.999927e-04`       `6.091215e+00`
    $2$       `6.280509e-13`      `-1.589331e-02`       `1.742094e+01`

  : Branch points at which bound states transition into a resonance for the Gauss $s$/$p$-wave example. The real part of $k$ is negligible whereas the imaginary part indicates slight off-zero branching discussed in \[sec:model2\].[]{data-label="tab:branchpoints"}

The proposed regularization plays a significant role near and at threshold values. Indeed, we have proven in appendix \[app:regularization\] the absence of additional, spurious zeros of $\det(\mathbold{F}(k,\lambda))$. All $\lambda_{22}$-dependent poles (which would otherwise coalesce with zeros in the origin at threshold values) are eliminated as well which makes continuation through the origin possible. We can not, however, guarantee smoothness of $\det(\mathbold{F}(k,\lambda))$ over the whole domain as demonstrated by the two stationary poles in figure \[fig:gauss\_l0\_l1\_c05\_smatrix\], where we show a comparison of $\det(\mathbold{S}(k,\lambda))$ and $\det(\mathbold{F}(k,\lambda))$.

\

\

\

Gauss potential with $p$-wave and $d$-wave coupling {#sec:model3}
---------------------------------------------------

The second model problem also uses Gauss potentials in both channels as well as the coupling and $\mu=1$. However, angular momenta for channel 1 and 2 were chosen to be $l_{1}=1$ and $l_{2}=2$. This choice is motivated by the previously mentioned distinction between $s$-waves and waves with higher angular momenta in \[sec:model2\]: In the $l=0$ case the bifurcation does not take place in the origin of the $k$-plane but on the imaginary axis below the origin, whereas for $l>0$ it does branch off in the origin [@Taylor2006]. By taking both angular momenta higher than zero we enforce branching of the coupled channel system in the origin and therefore an essentially different geometry for testing the regularized formulation than in the previous example.

Setting $\lambda_{11}=10$, $\lambda_{22}=30$ and $\lambda_{12}=0.3$ gives a system with 3 bound states. The approximations, analogous to the first example, are given in table \[tab:gauss\_l1\_l2\_c03\_bs\].

    $i$   $k_{i}$ (<span style="font-variant:small-caps;">Newton</span>)   $k_{i}$ (<span style="font-variant:small-caps;">Matscs</span>)   $k_{i}$ (<span style="font-variant:small-caps;">fin. diff.</span>)
  ----- ---------------------------------------------------------------- ---------------------------------------------------------------- --------------------------------------------------------------------
      0                                                  `3.796532e+00i`                                                  `3.796472e+00i`                                                      `3.796532e+00i`
      1                                                  `1.600083e+00i`                                                  `1.600083e+00i`                                                      `1.600127e+00i`
      2                                                  `6.599123e-01i`                                                  `6.599112e-01i`                                                      `6.603354e-01i`

  : Bound states of the coupled channel $p$- and $d$-wave example. We compare values obtained by Newton iterations (<span style="font-variant:small-caps;">Newton</span>) on $\det(\mathbold{F}(k,\lambda))$, (<span style="font-variant:small-caps;">Matscs</span>) and finite differences (<span style="font-variant:small-caps;">fin. diff.</span>)[]{data-label="tab:gauss_l1_l2_c03_bs"}

Again, these values were used as starting points for numerical continuation. The resulting curves are shown in figure \[fig:gauss\_l1\_l2\_c03\_continuation\]. The details of the coupled channel renormalized Numerov solver and AUTO continuation parameters were set to the same values as in the first example.

Also in this case the proposed regularization provides a viable way of removing poles in the origin at thresholds, hereby allowing continuation of the relevant zeros through the origin. This is shown in figure \[fig:gauss\_l1\_l2\_c03\_smatrix\]. This figure also shows the sharpness of nearby poles and zeros of $\det(\mathbold{S}(k,\lambda))$ compared to much smoother behavior of $\det(\mathbold{F}(k,\lambda))$ which positively influences the convergence of Newton iterations in the continuation process.

\

Discussion and Conclusions
==========================

In this article we have applied numerical continuation to track resonant states in a coupled channel Schödinger equation. In our approach, we track the zeros of a function that can be extracted from the numerical solutions of the coupled channel equation. This function, related to the $S$-matrix by a transformation, has a zero wherever the $S$-matrix has a pole and is amenable to numerical continuation. The analysis has also revealed that the function can have accidental poles which may slow down the convergence of the Newton iterations. However, for our model problems, we have not seen any difficulties caused by this.

We have tested this approach and solved a range of two channel model problems with robust results. The method automatically tracked all the resonant states and detected the bifurcation points. It also successfully traversed the threshold which indicates a successful regularization.

Our motivation to develop this method originates from chemical reactions that are modeled by the Born-Oppenheimer approximation. It is based on slow-fast separation between the motion of nuclei and electrons that is justified because of their mass difference. In this approximation the modeling is a two step process. In the first step the scattering of the electrons in the field of the nuclei is solved with the positions of the nuclei as parameter in this scattering problem. It is especially important to identify the resonances in the scattering problem. In the second step the molecular dynamics of the nuclei, that move on the potential surface given by the resonances determined in the first step, is solved.

This paper is focused on the first step where the resonances need to be identified in a family of scattering problems each with a different parameter, which is the position of the nuclei in these problems.

In our future work we aim to create a robust and reliable tool that scales to realistic systems.

Acknowledgments {#acknowledgments .unnumbered}
===============

We gratefully acknowledge the support from FWO-Vlaanderen through the project number G.0120.08.

Spherical functions {#app:spherical}
===================

Various solutions of the free radial wave equation are intensely used throughout this paper as they can describe the asymptotic behavior of the wave functions that solve the (non-free) Schrödinger equation. In an attempt to avoid ambiguities we summarize the used conventions in this appendix. For a far more thorough treatment we suggest [@Abramowitz1972] and [@Taylor2006], where similar conventions are used.

Spherical Riccati-Bessel functions are defined as solutions of the free radial Schrödinger equation $$\label{eq:radial_schrodinger_free}
    \left(\frac{d^{2}}{dr^{2}} - \frac{l(l+1)}{r^{2}} + k^{2} \right) u(r) = 0, \qquad k=\sqrt{2\mu E},\quad l\in\mathbb{N}.$$ The different, commonly used particular solutions of this equation are: ($x=kr$)

-   $\hat{j}_{l}(x) \ = \ x\,j_{l}(x) \ = \ \sqrt{x\frac{\pi}{2}}J_{l+1/2}(x)$, regular solutions, spherical Riccati-Bessel functions. E.g.: $\hat{j}_{0}(x) = \sin(x),\ \hat{j}_{1}(x) = \frac{1}{x}\sin(x)-\cos(x),\ \ldots$

-   $\hat{n}_{l}(x) \ = \ x\,n_{l}(x) \ = \ (-1)^{l+1}\sqrt{x\frac{\pi}{2}}J_{-l-1/2}(x)$, irregular, spherical Riccati-Neumann functions. E.g.: $\hat{n}_{0}(x) = -\cos(x),\ \hat{n}_{1}(x) = -\frac{1}{x}\cos(x)-\sin(x),\ \ldots$

-   $\hat{h}_{l}^{\pm}(x) \ = \ -\hat{n}_{l}(x) \pm i\hat{j}_{l}(x)$ Incoming ($-$) and outgoing ($+$) spherical Riccati-Hankel functions. E.g.: $\hat{h}_{0}^{\pm}(x) = e^{\pm ix},\ \hat{h}_{1}^{\pm}(x) = (\frac{1}{x} \mp i)e^{\pm ix},\ \ldots$

Regularization lemma {#app:regularization}
====================

This appendix technically describes the proof that the application of the regularizing procedure in section \[sec:Regularization\] does not introduce any zeros of the continuation function $\det(\mathbold{F}(k,\lambda))$ that can not be attributed to bound or resonant states. The proof involves the use of some well-known definitions and results from quantum mechanical scattering and gives a thorough explanation of the procedure yet is not necessary for the understanding and reproduction of the main results of our contribution. The conventions and definitions used are mostly those found in [@Taylor2006] and, on occasion in [@Newton1982; @Newton1960]. Section \[appsec:jost\] summarizes briefly a necessary minimal subset of definitions for the proof whereas section \[appsec:boundsproof\] proves why the regularizing procedure does not introduce false resonances or bound states.

For clarity we omit the potential parameter $\lambda$ in our notations as it does not play a role in the derivations; unless it is really necessary or clarifying. Following [@Taylor2006] we do however, indicate the explicit dependence of the wave functions on $k$, e.g. $\mathbold{\Psi}(k,r)$.

The Jost matrix and the $S$-matrix {#appsec:jost}
----------------------------------

The “regular solution” $\mathbold{\Phi}(k,r)$ of equation  is defined to behave exactly as $\hat{\mathbold{j}}(kr)$ at $r=0$. This fixes both $\mathbold{\Phi}$ and its normalization uniquely. At infinity the regular solution must tend to a linear combination of two other linearly independent solutions such as the free solutions $\hat{\mathbold{h}}^{\pm}$, hence allowing us to write $$\mathbold{\Phi}(k,r) \xrightarrow{r\to\infty} \frac{i}{2}\Big( \hat{\mathbold{h}}^{-}(kr)\mathbold{f}(k) \, - \, \hat{\mathbold{h}}^{+}(kr){\mathbold{f}(k)}^{*} \Big),$$ which defines the “Jost matrix” $\mathbold{f}(k)$ as a diagonal matrix that bundles so called “Jost functions” $f_{l}(k)$ $$\mathbold{f}(k) = \begin{pmatrix}
        f_{0}(k) \\
        & \ddots \\
        && f_{l_{\text{max}}}(k)
    \end{pmatrix}.$$ (it is known that ${f_{l}(k)}^{*} = f_{l}(-k)$ and hence ${\mathbold{f}(k)}^{*} = \mathbold{f}(-k)$) Since we have already written the asymptotic form of the “physical solution” $\mathbold{\Psi}(k,r)$ (which is calculated by the numerical solver) in terms of the $S$-matrix as in $$\mathbold{\Psi}(k,r) \xrightarrow{r\to\infty} \frac{i}{2}\left( \hat{\mathbold{h}}^{-}(kr) - \hat{\mathbold{h}}^{+}(kr)\mathbold{S}(k) \right),$$ it follows that the $S$-matrix can be expressed as the ratio of the Jost matrices $$\mathbold{S}(k) = \mathbold{f}(-k){\mathbold{f}(k)}^{-1}$$ and that $\mathbold{\Phi}$ and $\mathbold{\Psi}$ are related through the well-known expression[^1] $$\label{eq:physical_vs_regular}
    \mathbold{\Psi}(k,r) = \mathbold{\Phi}(k,r){\mathbold{f}(k)}^{-1}.$$ The previous equations illustrate a widely accepted result of bound and resonant states being represented, equivalently, by poles of the $S$-matrix and zeros of the Jost function. As we can not extract the Jost matrix directly from the numerical solution of but only a function that is proportional the Jost matrix, we have to use the poles of the $S$-matrix as the characterization of bound and resonant states.

The physical solution $\mathbold{\Psi}$ can also be written as the solution of the coupled channel partial wave version of the Lippmann-Schwinger equation $$    \label{eq:lippmann-schwinger}
    \mathbold{\Psi}(r) = \hat{\mathbold{j}}(kr) + \int_{0}^{\infty}dr'\, \mathbold{G}(r,r';k)\mathbold{U}(r')\mathbold{\Psi(r')},$$ where $\mathbold{G}(r,r';k)$ denotes the diagonal matrix of one-channel Green’s functions $G_{l_{i}}^{0}(r,r';k) = \frac{-1}{k}\hat{j}_{l}(kr_{<})\hat{h}_{l}^{+}(kr_{>})$, with $r_{<} = \min(r,r')$ and $r_{>} = \max(r,r')$ and where $\mathbold{U}=2\mu\mathbold{V}$. If we now combine with the asymptotic behavior of we find $$\mathbold{S}(k) - \mathbold{I} = \frac{-2i}{k} \int_{0}^{\infty} dr'\, \hat{\mathbold{j}}(kr') \mathbold{U}(r') \mathbold{\Psi}(r').$$ Inside the integral we replace $\mathbold{\Psi}$ with the regular solution $\mathbold{\Phi}$ through the identity  which leads to the equation $$\label{eq:145907062010}
    \mathbold{S}(k) - \mathbold{I} = \frac{-2i}{k}\left( \int_{0}^{\infty} dr'\, \hat{\mathbold{j}}(kr') \mathbold{U}(r') \mathbold{\Phi}(r') \right){\mathbold{f}(k)}^{-1}.$$

Bounds on the regularization {#appsec:boundsproof}
----------------------------

Equation  states that the $(l_{\text{max}}+1) \times (l_{\text{max}}+1)$ matrix $\mathbold{S}(k)-\mathbold{I}$ is proportional to a matrix of integrals multiplied by the inverse of the Jost matrix. To guarantee that the regularizing procedure from section \[sec:Regularization\] defined by eq.  as $$\label{appeq:det_multi_channel_regularized}
    \det(\mathbold{F}(k,\lambda)) = \left( \prod_{l=0}^{l_{\text{max}}} k^{2l+1} \right) \bigg/ \det(\mathbold{S}(k,\lambda)-\mathbold{I}),$$ does not introduce false bound states or resonances we must prove that zeros of the Jost matrix $\mathbold{f}(k)$ are the only contributors to the poles of expression  and as such, to the zeros of $\det(\mathbold{F}(k,\lambda))$. This requires the matrix of integrals in equation  to be bounded which is guaranteed by the following lemma.

An important remark has to be made ragarding the regularity constraints of the potential matrix $\mathbold{V}$. The proof of lemma \[lem:integralbound\] is a generalization of the one channel case where coupling terms in the potential do not occur. There, the boundary conditions that define the regular solution $\varphi_{l}$ and the iterative procedure of solving the one channel analogue of integral equation  for $\varphi_{l}$ indeed would give rise to proper bounds. Contrary to the one channel case however, the off-diagonal terms in $\mathbold{V}$ that couple channels with different angular momenta result in serious singularities in the origin $r=0$, hereby destroing the possibility for a bound without very restrictive assumptions on the potentials that would cancel out the singular behavior in the origin [@Newton1982; @Newton1960]. These restrictions can be lifted to the much more liberal requirement for the off-diagonal terms $\int_{0}^{\infty}dr\,r|V_{i\neq j}(r)| < \infty$ by rewriting the integral equation  with additional terms. The iterative procedure to solve the integral equation then gives the same regularity properties of $\mathbold{\Phi}$ as in the one channel case.

To retain clarity we demonstrate the bound using the more restrictive assumptions for the off-diagonal potential terms, keeping in mind that a similar bound can be found in the more liberal case of, for instance, the potentials studied in the numerical examples in section \[sec:results\].

\[lem:integralbound\] Let $\mathbold{\Phi}$ be the regular solution of eq.  where $\mathbold{U} = 2\mu\mathbold{V}$ the potential matrix satisfies the (restrictive) short range conditions $$\label{eq:192007062010}
        \begin{cases}
            \int_0^\infty dr \, r^{l_{i}+l_{j}+2}|V_{ij}(r)| < \infty & \forall i,j\\
            \int_{0}^{\infty} dr\, r^{-l_{i}+l_{j}+1}|V_{ij}(r)| < \infty & \forall i,j,
        \end{cases}$$ and let $\hat{\mathbold{j}}(kr)$ be the diagonal matrix of spherical Riccati-Bessel functions with corresponding $l_{i}$, then the integral $$\mathbold{\mathcal{I}} = \int_0^\infty dr \, \hat{\mathbold{j}}(kr)\mathbold{U}(r) \mathbold{\Phi}(r) \label{eq:integral_riccati_potential_regular},$$ is bounded.

Because of the particular boundary conditions for $r \rightarrow 0$ and the restrictions on the potential matrix $\mathbold{V}$ the regular solution $\mathbold{\Phi}$ fits the matrix integral equation (20.19) in [@Taylor2006] $$    \mathbold{\Phi}(r) = \hat{\mathbold{j}}(kr) + \int_0^r dr^\prime \, \mathbold{g}(r,r^\prime) \mathbold{U}(r^\prime) \mathbold{\Phi}(r^\prime),$$ where $\mathbold{g}(r,r^\prime)$ is a diagonal matrix composed of Green’s functions $$    g_{l_i}(r,r^\prime) = (1/k) \left( \hat j_{l_{i}}(kr)\hat n_{l_{i}}(kr^\prime) -\hat n_{l_{i}}(kr) \hat j_{l_{i}}(kr^\prime)\right),$$ on the diagonal.

This regular solution matrix can be written as a series $\mathbold{\Phi}(r) = \sum_n \mathbold{\Phi}^{(n)}(r)$ with the first term $\mathbold{\Phi}^{(0)} = \hat{\mathbold{j}}(kr)$ and the following recurrence relation between $\mathbold{\Phi}^{(n)}$ and $\mathbold{\Phi}^{(n-1)}$ $$\mathbold{\Phi}^{(n)}(r) = \int_0^r dr^\prime \, \mathbold{g}(r,r^\prime)\mathbold{U}(r^\prime) \mathbold{\Phi}^{(n-1)}(r^\prime).$$

The integral can then be written as $$\begin{aligned}
    \mathcal{I} &= \sum_{n=0}^{\infty} \int_0^\infty dr \, \hat{\mathbold{j}}(kr) \mathbold{U}(r) \mathbold{\Phi}^{(n)}(r) \\
    &= \int_0^\infty dr \, \hat{\mathbold{j}}(kr) \mathbold{U}(r) \hat{\mathbold{j}}(kr) + \int_0^\infty dr \, \hat{\mathbold{j}}(kr) \mathbold{U}(r) \int_0^r dr_{1} \, \mathbold{g}(r,r_{1}) \mathbold{U}(r_{1})\hat{\mathbold{j}}(kr_{1}) \\
    &+ \cdots + \int_0^\infty dr \, \hat{\mathbold{j}}(kr) \mathbold{U}(r) \int_0^r dr_1 \, \mathbold{g}(r,r_1) \mathbold{U}(r_1) \int_0^{r_1} dr_2 \, \cdots \int_0^{r_{n-1}} dr_n \, \mathbold{g}(r_{n-1},r_n)\mathbold{U}(r_n) \hat{\mathbold{j}}(kr_n) + \cdots\end{aligned}$$

Without loss of generality, let us look for example at a $2\times2$ coupled channel problem with $l_1$ in the first channel and $l_2$ in the second channel. The integral above is then $$\begin{aligned}
    \mathcal{I} &= \int_{0}^{\infty} dr \,
    \begin{pmatrix}
        \hat{j}_{l_1}(kr) & 0 \\
        0 & \hat{j}_{l_2}(kr)
    \end{pmatrix}
    \begin{pmatrix}
        U_{11}(r) &U_{12}(r) \\
        U_{21}(r) &U_{22}(r)
    \end{pmatrix}
    \begin{pmatrix}
        \hat{j}_{l_1}(kr) & 0 \\
        0 & \hat{j}_{l_2}(kr)
    \end{pmatrix} \\
    &+\int_{0}^{\infty} dr \,
    \begin{pmatrix}
        \hat{j}_{l_1}(kr) & 0 \\
        0 & \hat{j}_{l_2}(kr)
    \end{pmatrix}
    \begin{pmatrix}
        U_{11}(r) & U_{12}(r) \\
        U_{21}(r) & U_{22}(r)
    \end{pmatrix} \\
    &\times \int_0^r dr_1 \,
    \begin{pmatrix}
        g_{l_1}(r,r_1) & 0 \\
        0 & g_{l_2}(r,r_1)
    \end{pmatrix}
    \begin{pmatrix}
        U_{11}(r_1) &U_{12}(r_1) \\
        U_{21}(r_1) &U_{22}(r_1)
    \end{pmatrix}
    \begin{pmatrix}
        \hat{j}_{l_1}(kr_1) & 0 \\
        0 & \hat{j}_{l_2}(kr_1)
    \end{pmatrix} + \ldots .\end{aligned}$$

When we look at a particular element of this two by two example we see the contributions $$\begin{gathered}
    \mathcal{I}_{11} = \int_0^\infty dr \, \hat{j}_{l_1}(kr) U_{11}(r) \hat{j}_{l_1}(kr)
    + \int_0^\infty dr \, \hat{j}_{l_1}(kr) U_{11}(r) \int_0^r dr_1 \, g_{l_1}(r,r_1) U_{11}(r_1) \hat{j}_{l_1}(kr_1) \\
    + \int_0^\infty dr \, \hat{j}_{l_1}(kr) U_{12}(r) \int_0^r dr_1 \, g_{l_2}(r,r_1) U_{21}(r_1) \hat{j}_{l_1}(kr_1) + \ldots .\end{gathered}$$ We note that Riccati-Bessel functions that appear in the integral both have $l=l_1$. In a similar way, the matrix element $\mathcal{I}_{12}$ will have the $\hat{j}_{l_1}$ in the beginning of the integral and $\hat j_{l_2}$ at the end.

To show that this integral $\mathcal{I}_{11}$ is finite we make use of the bounds fond in [@Newton1982; @Newton1960; @Taylor2006]. There exist constants $C_{1}$ and $C_{2}$ such that $$\begin{aligned}
    |\hat{j}_l(kr)| &\le C_1 e^{|\text{Im}(k)|r}\left(\frac{|k|r}{1+|k|r}\right)^{l+1} \\
    |g_l(r,r^\prime)| &\le C_2 e^{|\text{Im}(k)|(r-r^\prime)}\frac{1}{|k|} \left(\frac{|k|r}{1+|k|r} \right)^{l+1}\left(\frac{|k|r^\prime}{1+|k|r^\prime} \right)^{-l}.
    \label{eq:bounds}\end{aligned}$$ After substitution we find $$\begin{aligned}
        |\mathcal{I}_{11}| \le k^{2l_1+2} &\bigg( C_1^2 \int_0^\infty dr \, e^{2|\text{Im}(k)|r}\left(\frac{r}{1+|k|r}\right)^{2l_1+2}|U_{11}(r)| \\
        & + C_1^2 \int_0^\infty dr \, e^{|\text{Im}(k)|r}\left(\frac{r}{1+|k|r}\right)^{l_1+1} |U_{11}(r) |
        C_2 e^{|\text{Im}(k)|r}\left(\frac{r}{1+|k|r}\right)^{l_1+1} \\
        & \times \int_0^r dr_1 \, C_2 e^{-|\text{Im}(k)|r_1}\left(\frac{r_1}{1+|k|r_1} \right)^{-l_1} |U_{11}(r_1)| e^{|\text{Im}(k)|r_1}\left(\frac{r_1}{1+|k|r}\right)^{l_1+1} \\
        &+ C_1^2 \int_0^\infty dr \, e^{ |\text{Im}(k)|r}\left(\frac{r}{1+|k|r}\right)^{l_1+1} |U_{12}(r) |
        C_2 e^{|\text{Im}(k)|r}\left(\frac{r}{1+|k|r}\right)^{l_2+1} \\
        & \times \int_0^r dr_1 \, C_2 e^{-|\text{Im}|(k)r_1}\left(\frac{r_1}{1+|k|r_1} \right)^{-l_2} |U_{21}(r_1)| e^{|\text{Im}(k)|r_1}\left(\frac{r_1}{1+|k|r}\right)^{l_1+1} \\
        & + \ldots \bigg) \\
        = k^{2l_1+2}& C_{11}, \qquad \text{for some constant } C_{11}.
    \end{aligned}$$ This series can be summed up because the potentials $V_{11}$, $V_{12}$, $V_{21}$ and $V_{22}$ satisfy the regularity condition . In a similar way we can find bounds for the other elements of $\mathcal{I}$.

In a more general case we have the bound $$\begin{pmatrix}
         |\mathcal{I}_{11}| & |\mathcal{I}_{12}| & \ldots \\ 
        |\mathcal{I}_{21}| & |\mathcal{I}_{22}| &\\
        \vdots  &  & \ddots\\
     \end{pmatrix}
    \le 
    \begin{pmatrix}
             k^{2l_1+2} C_{11} &  k^{l_1 + l_2 +2} C_{12} & \ldots \\
        k^{l_2 + l_1 +2} C_{21} &  k^{2 l_2 +2} C_{22} &  \\
        \vdots & & \ddots
     \end{pmatrix}.$$ This indeed means that the only way expression  can have a zero is through a zero of the Jost function, which ensures that no additional false “resonances” or “bound states” are introduced by using zeros of this function to locate them. Naturally every zero of $\mathbold{f}(k,\lambda)$ will result in a zero of $\det(\mathbold{F}(k,\lambda))$ which means we are able to track all resonances and bound states through the use of . This proves the equivalence of bound and resonant states and the zeros of $\det(\mathbold{F}(k,\lambda))$.

[^1]: Please note the different approaches studied in various reference works on the subject. Some authors such as [@Taylor2006] prefer the simpler approach of expression , whereas others like [@Newton1982] explicitly write certain factors without incorporating them, for instance, in the Jost function. In the single channel variant of for instance this would result in $\psi_{l}(k,r) = \frac{k^{l+1}}{(2l+1)!!}\frac{\phi_{l}(k,r)}{f_{l}(k)}$. Unless explicitly stated otherwise, we choose to follow the former convention.
