Fitting ultracold resonances without a fit

We present a numerical procedure allowing one to extract Feshbach resonance parameters from numerical calculations without relying on approximate fitting procedures. Our approach is based on a simple decomposition of the reactance matrix in terms of poles and residual background contribution, and can be applied to the general situation of inelastic overlapping resonances. A simple lineshape for overlapping inelastic resonances, equivalent to known results in the particular cases of isolated and overlapping elastic features, is also rigorously derived.


I. INTRODUCTION
Feshbach resonances are an ubiquitous phenomenon taking place when the energy of colliding particles is nearly degenerate with a quasi-bound level of the system. Resonances are traditionally observed as a function of collision energy or also by varying applied electromagnetic fields. In ultracold gases, it has been long realized that resonances occurring at extremely low energies can be induced by tuning an external magnetic field [1], allowing one to vary in a controlled way the effective interatomic interaction. The possibility to control the interaction using resonances is at the hearth of a series of groundbreaking experiments in the field of ultracold gases; See e.g. [2].
Feshbach resonances have been studied both theoretically and experimentally in a number of ultracold atomic systems [3]. Theory usually conveys information extracted from numerical models in a small number of lineshape parameters, namely position, width and the non-resonant (or background) contribution, extracted by fitting the numerical result to an analytical form. A celebrated example is the case of an isolated elastic resonance [1] where the s-wave scattering length diverges both positively and negatively as a function of the magnetic field B across the center B 0 : Here ∆ represents the magnetic width and the background scattering length a bg is the off-resonance value of a(B). A more general lineshape describing a pair and an arbitrary number of elastic overlapping resonances has been derived using formal scattering theory and quantum defect analysis in the Supplemental Material of [4] and in [5], respectively. In this case, the scattering length is still a purely real quantity showing a pair or a series of nearby divergences.
Isolated resonances in the presence of inelastic processes have been discussed in [6] and studied in greater detail in [7,8]. The main influence of inelastic loss channels is to make the scattering length a complex quantity and to suppress the pole in its real part, that no longer diverges [7,9]. In the presence of inelastic processes both the lineshape and the extraction of related resonance parameters become more involved. A relatively elaborate fitting strategy based on a three-points iterative procedure has been developed in [10]. However, like any other approach relying on fitting the lineshape to an analytical expression, the latter procedure presents few intrinsic drawbacks.
First, the precise value of the fitting parameters depends on the chosen fitting interval and bears an error bar. Second, the analytical lineshape can be a poor approximation to the fitted data, even to the extent that the fitting procedures may not converge [10]. Even if the least square iterations do converge, the physical meaning of the fitting parameters will become unclear if the analytical lineshape poorly describes the data. For instance, if the background scattering length presents an unusually rapid variation, defining such quantity as the off-resonance value of a(B) is not a well defined prescription. Finally, the analytical lineshape for overlapping resonances in the presence of inelastic processes is not yet available. This situation which indeed occurs frequently for collisions in excited states [11] is still missing a systematic procedure allowing one to summarize the numerical data in few theoretical parameters. A simple approach consisting in presenting graphics of the numerical results without accompanying tables of resonance parameters is not fully satisfactory, since reading quantitative information from a plot without having access to the rough data can be difficult or impossible.
In this paper, we present a formalism designed to overcome all of the above difficulties. Our approach allows one to define and extract in an algorithmically well defined way resonance parameters even in cases where the lineshape cannot be modeled by a simple expression. For the sake of better physical insight and clarity, our discussion is based on Feshbach resonance theory [12]. However, at a fundamental level the present strategy is based on a general pole approximation that could be taken as starting point in a similar fashion as scattering resonances can be defined as scattering-matrix poles [13]. Special attention is granted to stable ways to implement such pole decomposition, a numerically non trivial task in the presence of naturally diverging quantities such as the scattering length on resonance. As a by-product of the formalism, we also demonstrate that the complex scattering length can always be parametrized in a simple and intuitive way in terms of complex magnetic field locations and magnetic widths.
The paper is organized as follows. In section II we introduce the formalism building on Feshbach resonance theory. Section III proves that the present approach reduces under suitable conditions to results previously appeared in the literature.
In section IV we present an application to two situations of physical interest for 39 K, the case of overlapping and inelastic s-wave resonances and of finite-energy p-wave resonances.
A short conclusion ends this work.

II. FORMALISM
A. Basic ideas for single channel scattering The main idea can be illustrated with reference to the simple single channel resonant expression of the s-wave scattering length in the presence of an isolated resonance. Here we consider the most common case that the scattering length is tuned via an external magnetic field B, but the formalism applies equally well to a general external fields or variation of any physical parameter in the Hamiltonian capable of inducing resonances. In the case of a magnetic field, recalling (1), one can write : where we have defined p = −a bg ∆ the field-independent pole strength. Any variation of a with B beyond the pole contribution is included in the background term, which will therefore be considered a function of B.
The pole strength can be extracted from a theoretical calculation of a as It will be shown below how this limit can be performed in a numerically stable way; See equation (16). Once the p parameter is determined, one can compute the on-resonance value of the background scattering length according to Again, this limit can be computed in a stable way; See (19). We stress that the usual interpretation of the background scattering length as the limit value of a far away from resonance may become ambiguous in cases where a bg is not strictly constant or, even worse, is rapidly varying. On the converse, the quantity (4) is well defined in any case and denotes simply the value of a(B) at the resonance location after subtraction of the divergent part. We now generalize these ideas to multichannel scattering.

B. Pole decomposition for multichannel scattering
Les us consider a general scattering problem with N channel, N c closed channels, and N op = (N − N c ) open channels. We will further partition the N op open channels into N o degenerate (also termed in scattering theory super-elastic) channels with energy equal to the incoming channel, and Nō inelastic channels with energy release (or absorption). We use Feshbach theory, that formally partitions channel space into open and closed channel subspaces P and Q [14]. Scattering in the P subspace is supposed to be known, and the influence of the Q subspace is taken into account through the introduction of an effective Hamiltonian. The nature of P and Q subspaces depends on the specific problem at hand and for general formal manipulations does not even need to be specified [15]. A resonance occurs when the energy of bound states in Q subspace is quasi-degenerate with the energy of the colliding particles. Under such conditions the influence of the Q subspace on the total scattering amplitude becomes dramatic.
Feshbach formalism is most often applied to the scattering matrix [14], but also provides a partition of the reactance matrix [12]. Working with the reactance matrix is convenient since the resulting scattering matrix is automatically unitary, algebra is real, and the number of fitting parameters is minimal [12]. According to Feshbach theory, the reactance matrix at total energy E can be written where 1 P is the identity in P space and the resonant part of the reactance matrix has the form Here y is the N op × N c matrix formed by half predissociation amplitudes and the operator (E1 Q − E b − V) acts in the closed channel subspace, the Q-space of Feshbach theory. Zero energy is taken at the energy of the separated atoms at the given value of the applied field. The diagonal matrix E b represents the energy levels of the bare molecular states defining the Q-space. We will make the usual assumption that the position of such levels depends linearly on the applied field E b = −µ(B1 Q − B r ) with µ a diagonal matrix whose entries µ α are the differential magnetic moment of state α with respect to the one of two separated atoms in the initial channel. The entries B r,α of the diagonal matrix B r represent the bare resonance magnetic field, i.e the field at which level α crosses the incoming channel dissociation threshold. The coupling V contains in general an energy-independent direct coupling U between molecular states as well as indirect coupling through the continuum via the level shift operator δ E at energy E [16,17]. At this stage it is not necessary to distinguish the two contributions, only the Hermitic nature of V being of importance.
To proceed further with the derivation, inspired by the single-channel case, we factor out the magnetic moment in the denominator [1]. To this aim, we defineỹ = µ −1/2 y and B = B r + µ −1/2 (V − E1 Q )µ −1/2 . The resonant reactance matrix is expressed in terms of these redefined quantities as : If the magnetic moments all share the same sign the matrix µ 1/2 will be either purely real or purely imaginary. In this case, by its definition B will be real and symmetric, thus with real eigenvalues. If the magnetic moments present both negative or positive values the matrix B will in general be complex, yet still symmetric. We now insert (8) into (6). After some algebra, see proof in Appendix A, one can recast the K matrix in the form where the shifted resonance position matrixB is defined asB = B +ỹ t tan(ξ bg )ỹ. If there are no inelastic channels and if threshold collisions are considered (E → 0), all elements of tan(ξ bg ) as well as of the predissociation amplitude y tend to zero, henceB → B. TheB operator is still complex symmetric, is real if B is real, but unlike B in general it varies with magnetic field because of the background phase ξ bg appearing in its definition. It can therefore be brought a diagonal formb through a B-dependent complex orthogonal transformation C such thatB = CbC t and C t C = CC t = 1 [18]. With this transformation, equation (9) can be rewritten : in terms of a redefined predissociation amplitudeŷ = cos(ξ bg ) −1ỹ C t , which also depends on B.
The reactance matrix clearly diverges at the solutions, say b α , of the nonlinear equation B −b α (B) = 0. The following manipulations will turn out to be simpler if one isolates pole contributions with a field-independent pole residue. To this aim, let us define the constant amplitude vectorŷ (0) α =ŷ α (b α ) and rewrite (10) where the modified background is still regular at the poles and reduces to tan(ξ bg ) if ξ bg is strictly constant. According to the remark below (8) only real values of b α are possible if all entries in matrix µ have the same sign. If this is not the case, it can be shown that solutions b α are either real or occur in pairs of complex conjugate quantities. On a physical ground, such complex solutions arise when a pair of molecular levels with opposite slopes avoid each other exactly at the collision threshold, such that neither level crosses the latter. For purely elastic scattering, this unlikely case will result in a scattering length that may have a rapid variation, but will not diverge as B varies on the real axis. We will not consider this situation as a resonance condition and will henceforth assume that b α is always real, with the implicit assumption that the contribution from complex solutions b α , if any, is included in the background.
Since we are mainly interested in scattering near a collision threshold, it is convenient to factor out the low-energy limit behavior of the reactance matrix for channels with orbital angular momentum ℓ as follows [19]: The diagonal normalization matrix q −(ℓ+1/2) has by definition elements q −(ℓi+1/2) i for all channels degenerate with the incoming one, and to 1 for inelastic ones. The resulting normalized K-matrix is thus finite at threshold and is the main quantity needed to compute zero-energy scattering lengths. Performing such channel normalization and introducing the background scattering length matrix A = q −(ℓ+1/2) tan(ξ bg )q −(ℓ+1/2) as well as the resonance amplitudeȳ = q −(ℓ+1/2)ŷ(0) , equation (11) transforms intō Note that according to the Wigner laws [? ], by construction both A andȳ are finite at threshold. With the assumption above, the matrixK diverges at the real magnetic field locations b α . This should be contrasted with the behavior of the complex scattering length, that even in the presence of a pole inK on the real axis can in some cases present barely no structure. In order to locate the poles numerically, we find it more convenient to consider the inverse reactance matrix M =K −1 and to search for the zeros in the determinant det M. We first bracket zeros by a coarse search, then quickly refine their position using Ridder's method. Such nonlinear root-finding algorithm has as main advantages a superlinear convergence and that of keeping the root bracketed [20].
We now show that the residue as well as the background contribution at the pole position (recall that A will, in general, have a dependence on B) can be accurately determined. Concerning the residue, by definition : In order to evaluate the indeterminate 0 · ∞ product, we writeK = M −1 = adj M/ det M and use l'Hôpital rule.
Keeping into account continuity of the adjugate, one promptly obtains Equation (16) gives a numerically stable way to evaluate the residue in terms of M. The needed derivatives of the determinant can be evaluated using Jacobi's formula : The local background at the pole is defined by subtracting the divergent contribution from the global reactance matrix : Note that A loc contains contributions from the global background A as well as from all resonances but the α-th one. The lhs presents itself in a ∞ − ∞ form. In order to evaluate the limit, we write againK in terms of the adjugate and determinant of M. After bringing to common denominator and Taylor expanding to second order one obtains One proceeds similarly for the first derivative, which will also be needed. One starts from the definition and evaluates the indeterminate limit using Taylor expansion to get The required derivatives of the adjugate can be evaluated keeping into account that elements of the adjugate are in turn determinants of minors, and using again (17). Higher order derivatives of a determinant (and thus of the adjugate) follow by successive differentiation of (17). For instance, the second derivative can be calculated as Once the local background matrix (which includes contributions from nearby resonances) has been determined, the usually smooth global background and its derivative are respectively given by and In practice, it is clear that in these expansions one will not need to include all poles of K. One rather focuses on resonances occurring in a predefined magnetic field region and includes all poles such that the addition of additional ones has negligible influence on the parameters of the resonances of interest. In simpler terms, only resonances significantly overlapping with the ones under investigation must be taken into account. Is is also clear that without a stable and accurate algorithm to compute the M matrix and its derivatives the present formalism would be of little practical use. Our numerical approach for computing such derivatives is the subject of the following subsection.

C. Computing the derivatives
On a general ground, a basis representation of the wavefunction is particularly suited to compute derivatives of any order, since the latter can then be obtained in an essentially analytical fashion. In fact, while propagation algorithms also can be adapted for obtaining derivatives (see e.g. [22]), algorithmic complexity and computational cost both increase with the derivative order. In this work we exploit the pseudospectral representation termed spectral element method, a spatial grid representation presented for instance in [21] and references therein. In the present work we impose asymptotic R-matrix boundary conditions, but one could equally well work with log-derivative or scattering boundary conditions. In the discrete spatial grid representation, the time-independent multichannel Schrödinger equation reads : where the constant vector c on the rhs equal to zero everywhere but at last grid point r max , where it equals the unit matrix of dimension the total number of channels [21]. Elements of the Hamiltonian matrix in the spectral element representation can be found in [21]. By definition, the multichannel wavefunction Ψ at r max equals the R-matrix R. Once Ψ has been determined from equation (25), first derivative of the matrix wavefunction is obtained from differentiating (25) with respect to parameter B, and solving the resulting linear system: Again, the derivative of the R-matrix with respect to B corresponds to the computed ∂Ψ/∂B at r max . The K-matrix is ordinarily computed by imposing the asymptotic form of the wavefunction Ψ ∼ f − gK in terms of regular and irregular solutions f and g for the free particle problem. Since in our algorithm K is not needed, we rather determine directly its inverse M by matching to Ψ ∼ f M − g. The linear matching equations for M simply follow from equation (25) of [21] and read : the prime denotes the spatial derivative and all spatial arguments are tacitly taken at last point of the discretization grid. The first derivative of M with respect to the field is solution to the linear system obtained taking the derivative of (27) in which M, R, and ∂R ∂B are known. Derivatives of the reference functions with respect to the field are transformed into derivatives with respect to energy using chain rule and the known behavior of asymptotic channel energies with respect to B. The latter are calculated using a stable formula for both open and closed channels from [22], equations (11) and (12).
Formal generalization to the n-th derivative of Ψ is straightforward. One differentiates n times the Schrödinger equation, and solves for the highest order derivatives of the wavefunction k n k knowing the derivatives of Ψ up to order n − 1. Similarly, to determine ∂ n M/∂B n one differentiates n times the matching condition (27) k n k and solves for the highest order derivative knowing the derivatives of M up to order n − 1 as well as the derivatives of R up to order n. Higher order derivatives of reference functions with respect to B are obtained by successive differentiation of equations (11) and (12) in [22].

D. Observables : the energy-dependent scattering length
We now turn to the calculation of observable quantities. The scattering cross-sections at arbitrary energy can be expressed exactly as a function of the so-called energy-dependent scattering length [23]. Its definition for a single channel in terms of the diagonal scattering matrix element S jj : where q o is the channel wavevector. We consider here the general situation of a number N o of degenerate channels occurring for example when several partial waves contribute to the scattering process or in the presence of magnetic Zeeman sublevels in a vanishing magnetic fields. Hence we define a square complex scattering length matrix The off-diagonal elements of a are particularly relevant to parametrize anisotropic pseudopotentials [24] or in the context of spinor Bose-Einstein condensates [25]. We now write the relevant square matrices X = K, A, 1 as well as the generally non-square predissociation amplitudes y in blocked form for open degenerate and non-degenerate channels : With this decomposition and applying the projection technique of [19] one can express the middle matrix ratio on the rhs of (32) as a function of the blocks of the reactance matrix Transforming to the normalized quantities according to (13) and with the help of (14) we obtain : Making use of the Woodbury matrix identity [26] one can further write : Substituting (37) into (36), defining a new predissociation matrix uō = (1ōō − iAōō) −1ȳō and the complex resonance (36) can be recast in the more compact form With some additional algebra and the definitionū o = iA oō uō +ȳ o it is possible to put (38) in a form that clearly exhibits the divergence of the complex scattering length at the complex eigenvalues of the matrix B c : It should be noted that the matrix function B c and therefore its eigenvalues B c,i depend on the magnetic field B through the dependence of B c on Aōō, which in general vary with B. In order to carry out the pole expansion one can proceed in a similar way as it has been done for the K-matrix. First one writes where the determinant can be expressed as the product of the eigenvalues of its matrix argument Using regularity of the adjugate and omitting for notational simplicity the common energy dependence, the residue at the pole can be extracted as the limit where we have defined D α (B) = β =α (B − B c,β (B)). The local complex background at the pole b c,α is defined by subtracting the contribution of pole α from the full complex scattering length : Substituting (42) for the pole strength, bringing to common denominator and applying l'Hôpital's rule twice, one obtains This quantity contains contributions from all poles but α-th one. The global complex background at the pole b c,α is defined by subtracting all other pole contributions : resulting in where the superscript (k) denotes perturbation terms of order k on eigenvalues and eigenfunctions. Equating η powers order by order, one can compute terms up to perturbative order k = 2 in the eigenvalues, and identify B ′ c,α . The only difference with the standard procedure is that matrix B c is complex symmetric rather than self-adjoint such that a generalization of perturbation theory to non-Hermitian operators must be employed in order to calculate b (1) c,α and b (2) c,α . It should be stressed that the only approximation involved in this procedure is the neglect of the second derivative in the Taylor expansion of B c (B) near b c,α . This is justified since this term is usually small and it could in principle be included if needed, even if at the expenses of significant added computational complexity to obtain A ′′ oō . Further, derivatives of the adjugate C are evaluated as explained below equation (21). The termsū ′ o are obtained based on the definition ofū o and from the knowledge of A and A ′ . Also note the useful formula allowing the ratio in last but one term on the rhs of (46) to be easily determined.
Once the pole strength p α known for each resonance, one can cast the energy-dependent complex scattering length (39) in the intuitive form a quantity by construction regular at the poles. If A is strictly constant, it is clear that last two terms exactly cancel and the background reduces to Generalizing the usual definition from the single-channel case [1], we define the partial magnetic widths associated to pole α via the equation In the absence of inelastic channels, the matrix ∆ α is real, its imaginary part arises from inelastic processes. As a final note, we have undertaken the approach to collect all degenerate channels in the same open-channel manifold, which led us naturally to introduce matrices of scattering lengths and of magnetic widths. This description is most useful in the presence of degenerate partial waves or in equilibrium situations where all degenerate internal states are equally populated. In other experiments one is more interested in polarizing the atoms in a single quantum state and in either monitoring the evolution of the population or studying the collective behavior of the gaz in such specific state. In this case, it will be more useful to identify the open channel manifold with the unique state of interest, and to treat super-elastic channels on the same footing as the inelastic ones. All derivation above holds, the main difference being that the complex scattering length will not be real even in the absence of inelastic channels, since the initial quantum state will in general decay to degenerate channels.

III. COMPARISON WITH OTHER APPROACHES
The simplest case of an isolated elastic resonance in an elastic channel is trivially recovered from equation (49) by dropping the α index, putting p = −a bg ∆ and realizing that the resonance center b c is a real quantity.
For the case of several overlapping resonances in an elastic channel, simplicity of equation (49) may seem at first sight surprising as compared to equation (26) of Jachymski and Julienne [5]. In order to show that the present approach reduces to the quantum defect analysis presented in the latter work, we give a brief alternative derivation of their result starting from the present standard scattering theory. To this aim, let us specialize (5) to the case of a single open channel and consider by the sake of simplicity the low-energy limit : In the approach of Jachymski and Julienne, one first factors the energy dependence out of the predissociation amplitude y using the energy-dependent quantum defect amplitude C as y = C −1ȳ . Then one recognizes that the shift matrix in the denominator of K res contains an energy-dependent part proportional toȳȳ t via the quantum defect function tan(λ) [27]. This leads to : where D contains the Zeeman shift −µ (B1 Q − B r ) and all energy-independent couplings between channels in Q space. Making use of Woodbory's matrix identity to express the inverse, equation (54) is seen to be equivalent to Assuming absence of coupling between molecular states near threshold, matrix D takes the simple diagonal form D = −µ (B1 Q − B r ). Replacing in (55) leads after simple manipulations to equation (26) of [5]. Note that in the present work one first diagonalizes the denominator of K res at each value of the magnetic field. Therefore, our equation (49) can be considered in a sense the adiabatic version of equation (26)   We now consider an inelastic isolated resonance. With reference to the notation of equation (1) in [6] their lineshsape is recovered by first writing our complex pole p in the form p = −∆a bg , then by expressing the complex width ∆ in polar form as e 2iΦR ∆ el . Finally, one identifies a bg with a ∞ and the complex magnetic field location b c in the denominator of (49) with the term B 0 − i ∆ inel 2 in [6].
Barring from the choice of the notation, it is also easy to verify the equivalence of the lineshapes in [7] and [6] and thus in the present work.

IV. APPLICATIONS
A. Overlapping inelastic s-wave resonances As a first application, let us consider collisions of ultracold 39 K atoms, a system whose scattering properties have been determined very precisely first based on photoassociation then on magnetic Feshbach resonance spectroscopy [28]. We consider here the case of s-wave collisions for atoms in Zeeman sublevels correlating in zero magnetic field with hyperfine quantum numbers (F = 1, M = 1) and (F = 1, M = −1), respectively. This combination is not stable under spin-exchange collisions, since it can decay to a pair of (F = 1, M = 0) atoms. The complex scattering length shows two features, which have been reported before [29]. Based on the much stronger variation of the real part of a near 475 G than near 500 G, the authors of [29] concluded that only one feature is of resonant origin. A closer analysis of figure 1 shows however that the real part of a does not diverge near either feature, as expected on a general ground in the presence of inelastic processes. While the feature at higher field is indeed significantly more pronounced, one may still suspect the presence of two quasi-bound states underlying the observed structure.
A complimentary picture is provided by the trajectory of a(B) in the complex scattering-length plane as a function of the real parameter B; See figure 2. It it easy to verify that in the ideal case of an isolated resonance with strictly constant background such trajectory is a circle [7]. One can actually observe in the figure one large circle associated with the strong resonance, whereas closer inspection near the origin (see inset) highlights the presence of a second deformed loop described as the magnetic field crosses the weak feature, suggesting again the presence of a resonance pair. One can also observe that the global trajectory does not exactly close, since the values of the background scattering length away from resonance are not identical on the left and right resonance sides. Conclusive evidence about the nature of this feature is obtained by inspecting the elements of the full K-matrix depicted in figure 3. In agreement with equation (11), all elements of the reactance matrix diverge on resonance at the same location. Two poles are clearly observed, confirming that the relevant Q subspace contains two bare bound states coupled to the incoming channel. Following the numerical procedure outlined in section II, it is now easy to extract well defined lineshape parameters without resorting to a fit. Resulting numerical values for the complex background scattering length, field location and magnetic width can be found in table I.
It is also interesting to compare graphically the exact numerical result with the pole approximation (49) represented by the dots in figure 1. Note from the numerical values in the table that the background is nearly constant and can be represented locally by an interpolating linear function. As expected, the agreement with the non trivial behavior of the full numerical calculation is excellent, confirming that all elements of the theory have been computed and integrated properly. B. Energy-dependent p-wave scattering volume As a further example, let us consider collisions in higher-order partial waves. Specifically, we consider p-wave collisions between a pair of 39 K atoms prepared respectively in states |1 1 and |1 0 at a finite collision energy of E/k B = 1µK. We consider the projection M = 2 of the total axial angular momentum, corresponding to a pair of atoms with initial axial mechanical angular momentum m = 1. In this symmetry block, if p waves only are included, scattering is purely elastic. The energy-dependent scattering volume is defined according to the standard definition (32).
One can observe in figure 3 two series of overlapping resonances in the range of magnetic fields below 1000G, a doublet near 400 G and a triplet near 750 G. The largest resonance at 379 G is sufficiently broad to modify the local background near 750 G. Visually, the background appears to vary in a non trivial way with the magnetic field, in particular near the origin. The present approach has the advantage that no modeling is required in order to extract the resonance parameters. It is sufficient to locate precisely the resonances, extract the corresponding pole residues,  table I. If needed, the field-dependent background scattering length can be extracted by subtracting the pole contribution from the full scattering volume, i.e. last term on the rhs of (49) from the lhs. In spite of the fact that such operation involves taking the difference of large numbers and as such is in principle numerically flawed very close to the poles, we obtain a smooth background that, as anticipated, varies relatively fast in a non trivial way from negative to positive values.
Without performing a detailed analysis that is not the aim of the present methodological work, it is interesting to stress the effect of the spin-spin dipolar interaction. To this extent we set artificially to zero the dipolar coupling constant and perform the same steps as above. In the absence of anisotropic interactions, not only M but also M f = M a + M b is a good quantum number. This implies that resonances associated with molecular states with M f different from the one of the incoming atoms are completely uncoupled, i.e. have zero width. Indeed, one can remark  (49), with powers of tens indicated in brackets. Sample resonance parameters for ultracold collisions of K atoms in hyperfine states and partial wave specified in the first column are reported. For comparison, calculations with and without the spin-spin interaction are reported for p-wave scattering.
Channel   table) differ by no more less than 1%. Similarly, the spin-spin interaction only weakly (below 300 mG) shifts the resonance position. On the other side, the background scattering volume (and a-fortiori the magnetic width), has a significantly stronger variation with magnetic field in the presence of the dipolar interaction, highlighting the strong influence of the spin-spin interaction on the incoming wavefunction.

V. CONCLUSIONS
We have presented a formalism allowing resonance parameters to be extracted from a numerical calculation in a well determined algorithmic way. In our opinion the present approach fully solves a longstanding difficulty in the analysis of resonance spectra involving multiple overlapping features or a non trivial behavior of background scattering. Our results also provide a simple and intuitive lineshape for inelastic overlapping resonances that has been demonstrated starting from first principles. The machinery presented in this work can be incorporated in existing computer packages, provided that the underlying algorithm to solve the Schrödinger equation allows the wavefunction derivatives to be easily evaluated. As a complimentary step to the present paper that focuses on the lineshape parameters, it could be interesting to analyze the wavefunction at the poles of the reactance matrix in order to gain information on the nature of the resonant states.
To prove (39), we start by expanding (38) as Let us now transform second line. Using definitions of B c and uō, one can replace the termsȳ t o uō, u t oȳō , and y t o (1ōō − iAōō) −1ȳō by (B c − b)/i. Regrouping and simplifying, second line transforms into : Adding last term in first line of (A4) to result (A5) yields after simple matrix algebra similar to the manipulations in (A3) Along the same lines, one can simplify third and last line in (A4). For instance, for the third line one has where first equality follows upon using the definition of uō and using againȳ t o uō = −i(B c − b), and collecting common terms. Second equality results from simple matrix algebra. Finally, last line in (A4) can be expressed as Substituting equations (A5), (A7) and (A8) into (A4), collecting terms with common factor (B1 Q − B c ) −1 and defining, according to the main text,ū o = iA oō uō +ȳ o , (39) promptly follows.