Integral representation for scattering phase shifts via the phase-amplitude approach
D. Shu I. Simbotin R. Côté
physics.atom-ph
















Department of Physics, University of Connecticut,
   2152 Hillside Rd., Storrs, CT 06269-3046, USA






A novel integral representation for scattering phase shifts is
obtained based on a modified version of Milne's phase–amplitude
approach [W. E. Milne, Phys. Rev. 35, 863 (1930)].  We
replace Milne's nonlinear differential equation for the amplitude
function y with an equivalent linear equation for the envelope
ρ=y^2, which renders the integral representations highly amenable
to numerical implementations.  The phase shift is obtained directly
from Milne's phase function, which in turn is expressed in terms of
the envelope function.  The integral representation presented in this
work is fully general and it can be used for any type of scattering
potential, including the Coulomb potential.





xx.yy.zz, aa.bb.cc, etc.etc.etc.....






    
        Integral representation for scattering phase shifts 
        via the phase–amplitude approach
       
    R. Côté
    December 30, 2023
==============================================================================================================









§ INTRODUCTION




The phase–amplitude approach for Schrödinger's radial (or
one-dimensional) equation was pioneered by Milne <cit.>, and 
has since been used extensively in atomic and molecular physics
<cit.>,
in chemical physics
<cit.>,
and in other areas of physics
<cit.>.
Although it was originally intended for tackling bound states, the
phase–amplitude method is also applicable for scattering problems;
indeed, Milne's approach is especially suitable in the framework of
many-channel quantum defect theory
<cit.>
because it makes it possible to construct optimal reference functions
in each scattering channel.



The virtues of Milne's approach stem from the fact that phases and
amplitudes are quantities which are well behaved in the energy domain,
even across channel thresholds.  Moreover, the phase–amplitude method
allows for highly efficient numerical implementations, because the
direct computation of highly oscillatory wave functions can be avoided
entirely; instead, any solution of the radial equation is 
evaluated accurately in terms of the amplitude and phase functions,
which have a simple radial dependence.



In this work we use a modified version of Milne's method to derive
a novel integral representation which yields the true value of the
scattering phase shift, despite the modulo π ambiguity inherent in
its usual definition.  We begin with a simple integral representation
in Sec. <ref>, which we generalize for the Coulomb case in
Sec. <ref>.  We give a brief discussion of envelope and
phase functions in Sec. <ref>, with additional details
of the phase–amplitude approach presented in the Appendix.  In
Sec. <ref> we derive a two-envelope formula for phase
shifts, which is essential for high partial waves.  The computational
approach is discussed in Sec. <ref>, and we illustrate the
practical usefulness of our integral representations with numerical
applications in Sec. <ref>.  Finally, we give concluding
remarks in Sec. <ref>.









§ THEORY




We consider the radial Schrödinger equation for a spherically
symmetric potential V(R),

    ψ” = Uψ,            U=2μ(V_ eff-E),

where V_ eff(R)=V(R)+ℓ(ℓ+1)/2μ R^2 is the effective
potential,  μ is the reduced mass of the two particles undergoing
scattering, and  E>0 is the energy in the center-of-mass frame.  Atomic
units are used throughout.



Assuming that the interaction potential V(R) vanishes faster than
R^-1 asymptotically, the physical solution ψ(R) has the well
known asymptotic behavior

    ψ(R) sin(kR-ℓπ/2+δ_ℓ),

which yields the scattering phase shift δ_ℓ for partial wave
ℓ.  In the equation above, k=√(2μ E) is the momentum for
the relative motion.  The phase shift is usually obtained from
the matching condition,

    ψ(R) = A f(R) + B g(R),

where f and g are exact solutions of the radial equation,
specified by their asymptotic behavior,

    f(R)        sin(kR-ℓπ/2)
    
    
    g(R)        cos(kR-ℓπ/2).

From the equations above we have δ_ℓ=arctan(B/A), which
yields the phase shift modulo π, i.e.,
δ_ℓ∈[-π/2, π/2].  To simplify our notation,
we shall omit the label ℓ for ψ, f, g and other related
quantities, except for the phase shift δ_ℓ.







    

 §.§ Integral representation for the full phase shift

                  


We now derive an expression for δ_ℓ which does not rely on
the explicit evaluation of wave functions; instead, the phase shift
will be extracted from an R-dependent phase function.  Our approach
follows Milne's phase–amplitude method <cit.>, which we
review in appendices A, B, C.  We emphasize that the true
value of δ_ℓ will be obtained, despite the modulo π
ambiguity inherent in its customary definition.



We first introduce the envelope function

    ρ = f^2 + g^2,

with f and g  exact solutions of the radial
equation (<ref>) obeying the asymptotic
behavior (4, 5).  The phase function θ is constructed by
integrating

    θ' ≡dθ/dR = k/ρ.

We remark that θ(R) is defined up to an integration constant,
which can be chosen freely;  a judicious choice guided by the
computational strategy shall be made in Sec. <ref>.



As shown by Milne <cit.>, the general solution of the radial
equation (<ref>) can be represented exactly in terms
of ρ and θ.  In particular,
the physical solution reads

    ψ(R) = √(ρ(R))sin[θ(R)-θ(0)].

Note that ψ(R) vanishes explicitly at R=0, while
Eqs. (<ref>, <ref>) ensure a very simple
asymptotic behavior for ρ and θ,

    ρ(R)            1  
    
    θ(R)          kR + const.

We now define the reduced phase

    θ̃(R)  ≡θ(R) -kR,

and use Eqs. (<ref>, <ref>) to
find the asymptotic behavior of Milne's
parametrization (<ref>),

    ψ(R)  sin[kR+θ̃(∞)-θ(0)],

which is identical to the asymptotic behavior of ψ in
Eq. (<ref>).  Consequently, we obtain

    δ_ℓ -ℓπ/2 = θ̃(∞)-θ(0).

Making use of Eq. (<ref>) we have
θ(0)=θ̃(0), and
Eqs. (<ref>, <ref>) yield
θ̃'=k/ρ-k.  Hence,
Eq. (<ref>) can be recast as an integral
representation,

    δ_ℓ -ℓπ/2   =   θ̃(∞)-θ̃(0)
    
       =     k∫_0^∞ dr [1/ρ(r) -1].

We remark that the reduced phase θ̃(R) defined in
Eq. (<ref>) cannot be regarded as the nontrivial
phase contribution, as it also includes the Bessel contribution
(due to the centrifugal term).  Therefore, strictly speaking, the
equation above yields the full phase shift
(δ_ℓ-ℓπ/2).







 §.§ Generalization to potentials with a  Coulomb term

           

     


So far, we assumed that the potential V(R) vanishes faster than
R^-1 asymptotically.  In this section we consider the general case
when V(R) contains a Coulomb term,

    V_C(R) = Z_1Z_2/R.

with Z_1,2 the electric charges of the two colliding particles.
The remainder (V-V_C) of the interaction potential is responsible
for the phase shift δ_ℓ, which is obtained from the well
known asymptotic behavior

    ψ(R)  sin[kR-C/kln(2kR)
              +η_ℓ -ℓπ/2 +δ_ℓ],

where η_ℓ=Γ(ℓ+1 + iC/k) is the Coulomb phase
shift <cit.>, and C≡μ Z_1Z_2.
Following the same steps as in the previous section, we use again
Milne's parametrization (<ref>) for the physical wave
function, namely ψ(R)=√(ρ(R))sin[θ(R)-θ(0)],
with the phase θ behaving asymptotically as

    θ(R)    
                kR-C/kln(2kR) + const,

which is the generalized version of Eq. (<ref>).
Accordingly, the reduced phase is again defined such that it is finite
asymptotically,

    θ̃(R) ≡θ(R) -kR+C/kln(2kR),

and thus the full phase shift reads

    δ_ℓ - ℓπ/2 + η_ℓ
     = θ̃(∞)-θ(0),

which is the general form of Eq. (<ref>).



In the Coulomb case, an integral representation can only be written if
we divide the radial domain in two intervals.  Indeed, unlike the
previous section, θ̃ now diverges logarithmically when
R→0.  Thus, we shall employ θ(R) for R≤ R_0 and
θ̃(R) for R≥ R_0, with R_0 fixed arbitrarily.
Specifically, we have

    θ(R_0)-θ(0)    =    k∫_0^R_0dr/ρ(r)
    θ̃(∞)-θ̃(R_0)    =   
    k∫_R_0^∞dr[1/ρ(r)-1+C/k^2r].

Adding the two equations above, and  making use of
Eq. (<ref>) at R=R_0, we find

    δ_ℓ - ℓπ/2 + η_ℓ   =   θ̃(∞)-θ(0)
    
       =    k∫_0^R_0dr[1/ρ(r)-1] + C/kln(2kR_0)
    
       +    k∫_R_0^∞dr[1/ρ(r)-1+C/k^2r],

which represents the generalization of Eq. (<ref>).
Indeed, in the absence of the Coulomb term, i.e., setting C=0 in the
equations above, we have η_ℓ=0, and we recover the results of
Sec. <ref>.  We emphasize that the expression in
Eq. (<ref>) is independent of R_0, which we
illustrate with numerical results in Sec. <ref>.
Finally, we remark that
Eqs. (<ref>, <ref>) yield the true value
of δ_ℓ unambiguously (not modulo π), and in the case of
a purely Coulombic potential, Eq. (<ref>) yields the
true value of η_ℓ.








 §.§ Envelope and phase functions

              



In order to use the approach outlined above in numerical applications,
it is necessary to devise a reliable method for computing the envelope
directly, rather than using Eq. (<ref>).  As shown in
App. <ref>, ρ obeys a linear differential
equation,

    ρ”' - 4Uρ' - 2 U'ρ = 0.

Therefore, we now regard the envelope ρ=f^2+g^2 in
Eq. (<ref>) as a particular solution of
Eq. (<ref>).  Namely, we impose the asymptotic boundary
condition ρ(R)→1, which makes the solution unique.  As we shall
see in Sec. <ref>, we initialize the envelope at R=∞
and we propagate it inward; accordingly, we also propagate the reduced
phase θ̃ inward from R=∞.



Recall that Eq. (<ref>) allows for an 
integration constant to be chosen freely when constructing the phase
θ(R) or θ̃(R).  The integration constant can be
fixed, e.g., by setting the value of θ(0), or the value of
θ̃(∞).  We prefer the latter, which is suitable when
employing the inward propagation mentioned above; specifically, we
choose

    θ̃(∞) =0,

and thus Eq. (<ref>) reads

    δ_ℓ = ℓπ/2 -θ(0).




The reduced phase is constructed by direct integration; in the absence
of a Coulomb interaction term, we have

    θ̃(R)
     = k∫_R^∞dr ρ̃(r)/ρ(r),

where ρ̃ denotes the reduced envelope

    ρ̃≡ρ -1.

In the Coulomb case we make use of
θ̃'(r)=k/ρ(r)-k+C/kr, see
Eqs. (<ref>, <ref>), and thus the
reduced phase reads

    θ̃(R)
     = k∫_R^∞dr [1-1/ρ(r)-C/k^2r].

In order to show that the integral above is well defined and yields a
reduced phase obeying θ̃(∞)=0, and also to justify
that the full phase θ(R) has the asymptotic
behavior (<ref>), we write the envelope as an
asymptotic series,

    ρ(R) = 1 +∑_n≥1b_n/R^n.

Assuming the potential has the long-range behavior

    V(R) = ∑_n≥1C_n/R^n,      with  C_1=C=μ Z_1Z_2,

we substitute the ansatz (<ref>) in
Eq. (<ref>), and we obtain the coefficients b_n.  In
particular, for n=1 we have b_1=C/k^2, and thus the asymptotic
behavior of the envelope reads

    ρ(R) ≈ 1+C/k^2R +b_2/R^2 + ⋯.

Substituting the result above in Eq. (<ref>) yields
θ'≈ k-C/kR+𝒪(1/R^2), which upon
integration confirms Eq. (<ref>), while the
reduced phase in Eq. (<ref>) has the
asymptotic behavior

    θ̃(R) ≈b̃_1/R + b̃_2/R^2
     + ⋯0,

where the coefficients b̃_n are expressed in terms of b_n,
e.g., b̃_1=b_2-b_1^2.



For computational purposes, it is advantageous to employ the reduced
envelope

    ρ̃≡ρ -1 -C/k^2r.

Thus, we recast Eq. (<ref>) in terms of
ρ̃,

    θ̃(R) = k∫_R^∞dr 1/ρ(r)[(1-C/k^2r)ρ̃(r)
    -(C/k^2r)^2].

The full phase shift in Eq. (<ref>) can now be expressed
in a form suitable for computation,

    δ_ℓ   -   ℓπ/2 + η_ℓ
     = -θ(0)
       =    k∫_0^R_0dr/ρ(r)
     -kR_0 + C/kln(2kR_0) -θ̃(R_0).



We emphasize that θ̃ should be computed using
Eq. (<ref>), because the integrand in
Eq. (<ref>) suffers from cancellation at
large r.  Thus, in the asymptotic region, ρ̃ should be
obtained directly, rather than ρ itself; namely, we substitute
ρ=ρ̃+1+C/k^2R in Eq. (<ref>) which becomes
an equation for ρ̃.  The numerical approach for solving the
envelope equation is described in Sec. <ref>, where we show
that the entire asymptotic region can be treated in a numerically
exact fashion by mapping it onto a finite interval and using a
spectral Chebyshev method.







    

 §.§ Two-envelope formula for phase shifts

                 


We now derive an integral representation which makes it possible to
compute the phase shift δ_ℓ directly.  Recall that the
simple integral representations (<ref>) and
(<ref>) yield the full phase shift (including
Bessel and Coulomb contributions).  As shown in
Sec. <ref>, accurate values of δ_ℓ at high
ℓ cannot be obtained using Eq. (<ref>).  Thus, in
this section we formulate a two-envelope approach which can be used to
compute accurate phase shifts for all partial waves.  For the sake of
generality, we consider two different potentials V_1 and V_2, each
containing the same Coulomb interaction term (if present).  For a
given scattering energy, E>0, and a partial wave ℓ, we make use
of Eq. (<ref>) for each potential, and we
employ the convenient choice
θ̃_1(∞)=θ̃_2(∞) to find

    δ_ℓ^(2) - δ_ℓ^(1) = θ_1(0) - θ_2(0),

Using θ'_1, 2=k/ρ_1, 2, see
Eq. (<ref>), the phase difference above can be recast
as an integral,

    δ_ℓ^(2) -δ_ℓ^(1)
     = k∫_0^∞dr[1/ρ_2(r)-1/ρ_1(r)].

Although both θ_1, 2(R) diverge when R→∞, the
integral above is finite because the phase difference
θ_1(R)-θ_2(R)=θ̃_1(R)-θ̃_2(R) vanishes
asymptotically.  Note that in the asymptotic region we have
ρ_1≈ρ_2, which can also hold in the inner region if
V_1≈ V_2.  Thus, Eq. (<ref>) will suffer from
catastrophic cancellation,  rendering it unsuitable for
numerical applications.  Nevertheless, we show next that  a
computationally robust integral representation based on the
two-envelope approach can be formulated.



We choose V_1≡ V_ ref as a reference potential (with the
corresponding effective potential including both the centrifugal and
Coulomb terms, see below), while V_2≡ V=V_ ref+V̂ is the
full interaction potential.  The reduced envelope and phase are now
defined relative to the corresponding reference quantities,

    ρ̂   =   ρ-ρ_ ref
    θ̂   =   θ-θ_ ref.

We  employ a nontrivial reference problem by setting

    U_ ref(R) = -k^2 + ℓ(ℓ+1)/R^2 +2C/R,
         C=μ Z_1Z_2,

and we use of Eqs. (<ref>, <ref>) with
U_1=U_ ref given above and U_2=U=U_ ref+2μV̂.
Thus, the Bessel phase shift (-ℓπ/2) and the Coulomb phase
shift η_ℓ are both eliminated, and the phase shift
δ_ℓ reads

    δ_ℓ   =    -θ̂(0)
       =    -k∫_0^∞drρ̂(r)/ρ(r)ρ_ ref(r).

The reference envelope is the solution of Eq. (<ref>),

    ρ_ ref”' - 4U_ ref ρ'_ ref - 2 U_ ref' ρ_ ref
     = 0,

with U_ ref(R) given in Eq. (<ref>), while the reduced
envelope obeys a non-homogeneous differential equation,

    ρ̂”' - 4Uρ̂' - 2U'ρ̂= 4Ûρ'_ ref+2Û'ρ_ ref,

which was obtained by combining Eqs. (<ref>)
and (<ref>).  In the equation above we used the notation
Û≡ U-U_ ref=2μV̂.
Note that in the absence of a Coulomb term, we use U_ ref=
-k^2+ℓ(ℓ+1)/R^2, with V_ ref=0 and V̂=V,
which is illustrated with numerical results in
Sec. <ref>.  Finally, we remark that in the absence
of a Coulomb term one can also use the trivial choice U_
  ref=-k^2 with θ_ ref(R)=kR, which yields θ̂
identical to θ̃ in Eq. (<ref>), thus
recovering the integral representation of the full phase shift given
in Sec. <ref>.



In practical applications, one first solves Eq. (<ref>) for
the reference envelope, which is subsequently used in
Eq. (<ref>).  The latter is solved to obtain
ρ̂, and thus the full envelope is obtained: ρ=ρ_
  ref+ρ̂.  We remark that the numerical approach used for the
homogeneous envelope equation, see Sec. <ref>, can also be
employed for the non-homogeneous differential
equation (<ref>).  The two-envelope approach is
fully general, but is especially useful when Û=2μV̂ is
small, such that |ρ̂|≪ρ≈ρ_ ref, and thus
θ̂ and δ_ℓ will also be small.







§ COMPUTATIONAL APPROACH


          



The integral representations derived in this article are valid in
general; however, they are useful in numerical applications only if
the integrands are well behaved.  Specifically, the envelope should
behave in a non-oscillatory fashion, which in turn ensures the
smoothness of the phase function.  In practice, there is considerable
difficulty in finding the smooth envelope
<cit.>, because the general solution of
the envelope equation has an oscillatory behavior, as mentioned in
App. <ref>.  In this section, a novel computational
strategy for obtaining the smooth envelope in the asymptotic
region is presented.  We emphasize that finding the unique, smooth
solution is critically important for the efficiency and accuracy of
numerical schemes using our integral representations, and for using
the phase–amplitude method in general.



For clarity, we assume that V(R) vanishes faster than R^-1
asymptotically.  As discussed in Sec. <ref> (see also
App. <ref>), we employ the initial condition ρ=1
at R=∞ when solving Eq. (<ref>).  However, rather
than using the radial variable, it is highly advantageous to
reformulate the envelope equation by mapping the asymptotic
radial domain onto a finite interval.  We developed a convenient and
efficient numerical implementation based on a simple change of
variable,

    x=1/R,

which allows to take fully into account the long-range tail of any
potential.  Thus, the infinite radial interval R_1<R<∞ is now
mapped on a compact interval, x_1>x>0, with x_1=1/R_1.  The
boundary R_1, i.e., the size of the interval [0, x_1], will be
chosen to ensure the desired level of accuracy.




We now regard the envelope as an x-dependent function, and we
present the computational approach for finding the smooth solution of
the envelope equation inside the interval [0, x_1].  First, we
rewrite Eq. (<ref>) using x=1/R as the
independent variable,

    x^4⃛ρ +6x^3ρ̈+6x^2ρ̇-4Uρ̇-2U̇ρ = 0,

where dots above symbols denote derivatives with respect to x,
e.g., ρ̇=dρ/dx.  Recall that U=2μ V_ eff-k^2.
Next, we define

    u ≡ρ̇,

and we regard it as the unknown function.  Making use of the
initial condition ρ=1 at x=0, we write

    ρ(x)=1+∫_0^xdt u(t),

which we substitute in the last term of Eq. (<ref>) to obtain

    ( 4k^2 -8μ V_ eff
        + x^4D_x^2 +6x^3D_x +6x^2 -4μV̇_ effS_x)u
      = 4μV̇_ eff.

The operators D_x and S_x read

    D_xu=u̇=du/dx,        S_xu=∫_0^xdt u(t).



We solve Eq. (<ref>) using a spectral Chebyshev
method <cit.>,
i.e., we employ a small number of Chebyshev polynomials T_n(x) with
n=0,1,2,…,N-1, which are mapped onto the interval [0,x_1].
We expand the unknown function u(x) in the Chebyshev basis, and the
operators D_x and S_x are represented as finite (N× N)
matrices <cit.>.  Thus, Eq. (<ref>) becomes a
simple linear system,

    𝐌 𝐀 = 𝐁,

where the column 𝐀  contains the  Chebyshev coefficients of our unknown function,

    u(x) = ∑_n=0^N-1 A_n T_n(x),

𝐁 contains the Chebyshev coefficients for the expansion

    4μV̇_ eff(x) =  ∑_n=0^N-1 B_n T_n(x),

and 𝐌 is the matrix of the operator in Eq. (<ref>),

    M=4k^2 -8μ V_ eff + x^4D_x^2 +6x^3D_x +6x^2 -4μV̇_ effS_x.

M is singular, but its matrix (𝐌) in the small Chebyshev
basis is well conditioned and thus yields highly accurate solutions:
𝐀=𝐌^-1𝐁.  The smooth envelope is obtained as the unique
solution, because all the other solutions oscillate infinitely fast
near x=0 (R→∞) and they are thus eliminated; indeed, highly
oscillatory behavior cannot be accommodated by the small number of
Chebyshev polynomials.  We emphasize that the linearity of the
envelope equation is crucially important for the feasibility of the
approach presented here.  Finally, the solution obtained inside the
interval [0, x_1] can now be used to initialize the propagation for
x>x_1, i.e., R<R_1.  A detailed description of our computational
approach will be published elsewhere.









The advantage of the phase–amplitude method combined with the change
of variable x=1/R is readily apparent in
Fig. <ref>, where we show the x-dependence of the
reduced phase θ̃(x) and the envelope ρ(x) along with
|ψ|^2 for ℓ=475 and E=0.002 a.u., for the potential energy
used in Sec. <ref>, see Eq. (<ref>).  Note that the
wave function in Fig. <ref> was evaluated numerically
using Eq. (<ref>), i.e.,
ψ=√(ρ)sin(θ-θ_R=0).



When the potential has a repulsive wall at short range, the
formulation of the phase–amplitude approach based on the change of
variable x=1/R can be employed for the entire radial domain, as
shown in Fig. <ref>.  Indeed, the repulsive wall makes
it possible to stop the inward propagation at R_ min>0, which
corresponds to a finite value x_ max=1/R_ min.  We
remark that, when the inward propagation of the reduced phase
θ̃ approaches the repulsive wall, it is convenient to
convert it to the full phase θ using
Eq. (<ref>) at a point R_0 just outside the inner
wall; for the remainder of the radial domain (R<R_0), one should
propagate θ(R) instead of θ̃(R), because the former
converges much faster than the latter.  Indeed, if V_ eff→+∞
when R→0, we have ρ→∞ and
θ'→0, while θ̃'→-k.  In other words,
θ(0) should be computed as

    θ(0) = kR_0
     +k∫_R_0^∞dr ρ̃(r)/ρ(r)
     -k∫_R_ min^R_0dr/ρ(r),

which is independent of R_0.  We remark that the integration need
not extend fully to R=0, because R_ min is chosen inside the
repulsive wall to ensure the contribution of the interval 0<R<R_
  min is entirely negligible.  Consequently, the radial domain can
be safely restricted to R_ min<R<∞, which is mapped onto
the compact interval x_ max>x>0.






§ EXAMPLES

  

We now apply the integral representations and show that they yield
highly accurate results.
Our first example is the Coulomb potential, which we use as a test
case for the integral representation (<ref>).  The
two-envelope formula, see Eqs. (<ref>)
and (<ref>), will be employed and tested in
Sec. <ref>.







 §.§ The Coulomb potential




In the case of a purely Coulombic potential, Eq. (<ref>)
yields the Coulomb phase shift,

    η_ℓ =ℓπ/2 -kR_0 + C/kln(2kR_0) 
     +k∫_0^R_0dr/ρ(r) -θ̃(R_0),

with θ̃ given Eq. (<ref>).
The result above is independent of R_0, as depicted in
Fig. <ref>.  Indeed, we show that our approach is robust
and accurate by comparing the value of η_ℓ obtained using
Eq. (<ref>) with the exact value
η^Γ_ℓ≡Γ(ℓ+1+iC/k) for C=-1,
k=0.1 and ℓ=5.  Our integral representation yields the value
η_ℓ=-20.22421961527, while the analytical expression gives its
value modulo 2π, namely η^Γ_ℓ=-1.3746636937335435.
Their difference equals an integer multiple of 2π to a high degree
of precision: (η_ℓ-η^Γ_ℓ)/2π=-3(1±10^-13).
Figure <ref> depicts the relative error
|(η_ℓ-η^Γ_ℓ+6π)/η_ℓ| as a function of
R_0.













 §.§ Results for an inter-atomic potential with long-range
  behavior of the type V(R)∼-C_3/R^3


              




For our second example we use both integral representations, i.e.,
Eq. (<ref>) for the full phase shift and the two-envelope
formula (<ref>) which yields the phase shift
directly.  We employ the potential energy

    V(R) = C_ wallexp(-R/R_ wall) 
     - C_3/R^3+R_ core^3,

with C_ wall=10, R_ wall=1, R_ core=5 and
C_3=18 (all in atomic units), and the reduced mass μ=m/2,
where m is the mass of ^88Sr.  We computed phase shifts for
E=0.01  a.u.≈0.272 eV for a wide range of partial waves.
The upper panel in Fig. <ref> depicts the
ℓ dependence of the phase shift, while the lower panel shows the
partial-wave terms of the elastic cross section,
σ_ℓ=4π/k^2(2ℓ+1)sin^2δ_ℓ.  An exceedingly
large number of partial waves contribute to the cross section; note
that the dominant contribution stems from very high partial waves
(ℓ>5000).  For a fully converged value of the elastic cross
section, we have computed phase shifts up to ℓ=10^5.










Recall that our integral representations yield the true value of the
phase shift (not modulo π) which has a rather simple
ℓ-dependence; this suggests that interpolation schemes could be
used to drastically reduce the number of partial waves for which phase
shifts need to be computed.  This added advantage is illustrated in
Fig. <ref>, where  we compare the true phase shift with
its modulo π version.  Moreover, our approach is not restricted to
integer values of ℓ, and makes it possible to use non-integer
values of ℓ as interpolation points; thus, highly accurate
interpolation methods with non-uniform grids, e.g., Chebyshev
interpolation, can be employed.













Regarding the practical aspects of the computation, some remarks are
in order.  We first emphasize that, using the computational approach
presented in Sec. <ref>, the simple integral
representation (<ref>) and the two-envelope
formula (<ref>) can be implemented numerically
such that they both yield accurate results.  However, when using the
integral representation (<ref>) for the full phase
shift, the value of δ_ℓ obtained from
Eq. (<ref>) will gradually lose precision at
very high ℓ.  Indeed, for ℓ→∞ we have
δ_ℓ→0 while θ(0)≈ℓπ/2.  Therefore,
the simple integral representation (<ref>) must be
avoided at high ℓ because of the catastrophic cancellation in
Eq. (<ref>), even though θ(0) can still
be computed accurately.  Consequently, when ℓ becomes extremely
large, δ_ℓ should instead be computed using the two-envelope formula
derived in Sec. <ref>.



Although the two-envelope formula is highly accurate for all partial
waves, Eq. (<ref>) has the advantage of much greater
simplicity and could be used at low ℓ, provided that it is
sufficiently accurate.  A simple rule of thumb exists for finding the
highest partial wave for which
Eqs. (<ref>, <ref>) yield accurate
results.  Namely, the two-envelope formula needs to replace the simple
formula only if ℓ is high enough for the centrifugal term to
become dominant over V(R).  However, at low ℓ the simple
formula is highly accurate, as we show next.











First, we performed a test for the Bessel case, V_ eff
=ℓ(ℓ+1)/2μ R^2, with V=0, in order to show that the
two-envelope formula is robust and accurate even for very high ℓ.
Namely, we computed directly the vanishingly small difference between
the Bessel phase shifts for ℓ_1 and ℓ_2 which are nearly
equal.  We used ℓ_1=ℓ and ℓ_2=ℓ-1/ℓ^2, for
which the exact value of
δ_ℓ=δ_B(ℓ_2)-δ_B(ℓ_1) is
δ_ℓ=π/2ℓ^2, where δ_B(ℓ)=-ℓπ/2
denotes the Bessel phase shift.  The relative error of the numerical
value δ_ℓ obtained with the two-envelope formula is shown in
Fig. <ref>, confirming that high accuracy is
preserved at high ℓ, despite the smallness of δ_ℓ.  In
contrast, if the exact values of δ_B(ℓ_1,2) are
subtracted numerically, the loss of accuracy is significant and
becomes catastrophic at very high ℓ, see
Fig. <ref>.  This also illustrates that the failure
of Eq. (<ref>) at high ℓ cannot be
avoided, as it is due to the cancellation of nearly equal quantities.
Nevertheless, Eq. (<ref>) is sufficiently
accurate at low ℓ.  Indeed, we performed a test for the simple
integral representation (<ref>).  Assuming the
two-envelope formula is numerically exact,  the relative error for
δ_ℓ computed using
Eqs. (<ref>, <ref>) for the
nontrivial potential energy (<ref>) is shown in
Fig. <ref>, which makes it readily apparent that the
simple formula is highly accurate for partial waves ℓ⪅
600, while significant loss of precision only occurs for much higher
partial waves.





§ CONCLUSIONS AND FUTURE WORK




In this work we presented a novel integral representation for
scattering phase shifts.  A simple version, Eq. (<ref>),
yields the so called full phase shift, while the two-envelope formula
can be used to obtain δ_ℓ directly.  We emphasize that,
unlike standard approaches which only yield δ_ℓ
mod π, our integral representations give the true value of the
phase shift.  The integral representations are very useful in
numerical applications, as shown in Sec. <ref>; indeed,
they can be easily implemented to obtain highly accurate results.
Moreover, our approach is based on the phase–amplitude method, which
avoids the explicit computation of highly oscillatory wave functions;
thus, our numerical implementations are very efficient, and are
consistently accurate even for very high partial waves.



We remark that our integral representations are always valid, but the
computational advantages mentioned above rely on the assumption that
the envelope is globally smooth.  However, a globally smooth solution
of the envelope equation (<ref>) does not always exist;
indeed, as is well known <cit.>, when a barrier
separates two classically allowed regions, the solution which is
smooth on one side of the barrier will oscillate on the other side.
To overcome this computational difficulty, it is necessary to develop
an optimization procedure for finding the locally smooth envelopes.
In addition, one needs to find a way to combine the different
solutions which are locally optimized in each classically allowed
region.  This is a nontrivial task, as it requires certain functional
relationships be established between the envelopes and phases
corresponding to disjoint regions which are separated by a barrier.
If such an approach could be rendered computationally feasible, it
would yield the true phase-shift even in the absence of a globally
smooth envelope.







This work was partially supported by the National Science
Foundation Grant No. PHY-1415560 (DS) and the MURI US Army Research
Office Grant No. W911NF-14-1-0378 (IS and RC).









     


     

§ THE PHASE–AMPLITUDE APPROACH


      

A brief overview of Milne's phase–amplitude approach is presented here.
According to Milne <cit.>, the general solution ψ of
the radial equation (<ref>) can be expressed in terms
of an amplitude y and a phase θ,

    ψ(R) = c  y(R) sin[θ(R)+θ_0],

where c and θ_0 are arbitrary constants.  The amplitude satisfies
the nonlinear equation

    y” = Uy + q^2/y^3,

with U=2μ(V_ eff-E), and the phase θ(R) is
constructed by integrating

    θ'=q/y^2.

In the equations above, q is arbitrary; the only restriction is
q^2>0 in the amplitude equation (<ref>).  We emphasize
that the amplitude and phase appearing in Milne's
parametrization (<ref>) are not unique; indeed, any
solution y of Eq. (<ref>) together with the associated
phase θ will give a valid representation of ψ.  This
undermines the advantage of Milne's method in numerical applications,
because the general solution y(R) of Milne's nonlinear equation has
an oscillatory behavior in classically allowed regions, and the unique
smooth amplitude is very difficult to find
<cit.>.  Despite this
difficulty, Milne's nonlinear equation (<ref>) has long
been used for computational work.  We remark that an equivalent
formulation based on a linear equation exists 
<cit.>, but it remained overlooked in the
physics community until
recently <cit.>.  As we show in
App. <ref>, the linear equation can be obtained by
simply replacing the amplitude with the envelope function
ρ,

    ρ(R)=y^2(R).

The envelope obeys a third order linear differential equation,

    ρ”' = 4Uρ' + 2 U'ρ.

In App. <ref> we present two
different derivations for the envelope equation (<ref>),
and in App. <ref> we discuss its equivalence with
Milne's nonlinear equation (<ref>).



Although Eq. (<ref>) is of third order, its linearity makes
it much more convenient than Milne's nonlinear
equation (<ref>).  However, finding the non-oscillatory
solution is still a difficult task.  To overcome this obstacle, we
devised a computational strategy for scattering problems (E>0) which
yields the smooth envelope in the asymptotic region, see
Sec. <ref>.






§  THE ENVELOPE EQUATION






We present here two different derivations of the envelope
equation (<ref>, <ref>).  The first
derivation is very brief; namely, we substitute y=√(ρ) in
equation (<ref>) and we find

    ρ” =  2Uρ + 1/ρ[1/2 (ρ')^2 + 2q^2].

Next, we multiply both sides  by ρ to obtain

    ρρ” = 2Uρ^2 + 1/2 (ρ')^2 + 2q^2,

which is still a nonlinear equation.  However, we now take the
derivatives of both sides,

    ρρ”' = 4Uρρ' + 2U'ρ^2,

and we divide by ρ  to finally obtain Eq. (<ref>).  



The second approach is similar to Milne's derivation <cit.> of
Eq. (<ref>).  Namely, we consider two solutions (ϕ and
χ) of the radial Schrödinger equation (<ref>),
and we try to find a differential equation for their product

    p ≡ϕχ.

Making use of Eq. (<ref>) for ϕ and χ, we obtain

    p” = 2Up + 2ϕ'χ',

and we now evaluate its derivative,

    p”' = 2U'p + 2Up' +2U(ϕχ'+ϕ'χ),

where we recognize p'=ϕχ'+ϕ'χ, and we find again the
envelope equation (<ref>),

    p”' =  4Up' + 2U'p.

As the product of any two solutions of the radial
equation obeys the envelope equation, ϕ^2 and χ^2 are also
valid solutions of Eq. (<ref>), as well as any linear
combination of ϕ^2, χ^2 and ϕχ
<cit.>.
In particular, ρ=ϕ^2+χ^2 is  a
valid solution, which corresponds to Milne's ansatz y=√(ϕ^2+χ^2).







§ MILNE'S AMPLITUDE EQUATION AS A CONSTRAINT FOR THE ENVELOPE EQUATION





The parameter q appears explicitly in Milne's nonlinear
equation (<ref>).  However, q is absent from the envelope
equation (<ref>), even though it is used when integrating
Eq. (<ref>) to obtain the phase θ.  This creates
some ambiguity, which stems from the fact Eq. (<ref>) is a
second order differential equation, while the
Eq. (<ref>) is of third order.  We now try to dispel the
ambiguity and show that the two equations are equivalent.  We first
remark that although q does not appear in Eq. (<ref>),
it should be assumed implicitly; indeed, if we recast
Eq. (<ref>) in the form

    1/2ρρ” - Uρ^2 - 1/4 (ρ')^2  =  q^2,

the expression on the left hand side can be interpreted as an
invariant of the envelope equation (<ref>), and any
solution ρ will also obey Eq. (<ref>) with a
particular value of q^2.  Recall that the equation above is
equivalent with Milne's equation (<ref>), which can thus be
regarded as a constraint for the envelope equation.  Indeed, as we
discuss below, Eq. (<ref>) should be used to enforce
the correct initial conditions for ρ, such that they correspond to
a fixed value for q.



To fully clarify the equivalence of Eqs. (<ref>) and
(<ref>), let us compare the sets of initial conditions
required in each case.  When we initialize ρ at R=R_0, we
consider given

    ρ(R_0),     ρ'(R_0),     ρ”(R_0),

which can be used in the constraint equation (<ref>)
evaluated at R=R_0 to obtain the value of q.  Conversely, if q
is considered given, we have

    ρ(R_0),     ρ'(R_0),       q =  fixed,

which we commonly employ in practice.
Equation (<ref>) is now used to obtain ρ”(R_0),
and thus initialize the solution of the envelope equation.
Equivalently, for Milne's amplitude equation we consider given

    y(R_0),      y'(R_0),      q =  fixed.

Moreover, the families of solutions for different values of q are
all equivalent.  Indeed, if y_1 is a solution of
Eq. (<ref>) for a given parameter q_1, then
y_2=(q_2/q_1)^1/2y_1 is a solution for q_2.  Similarly, we
have ρ_2/q_2=ρ_1/q_1, and θ'_2=θ'_1.  Varying the
parameter q is entirely redundant, as the phase θ remains
unchanged, therefore justifying the convenient choice q=k used
throughout this article.




Finally, we make use of the constraint in Eq. (<ref>)
to show that the choice q=k is consistent with the initial condition
ρ=1 at R=∞, which gives a convenient normalization for the
envelope.  Indeed, when R→∞, we have V_ eff(R)→0 and
thus U(R) ≈ U(∞) = -k^2 and U'≈ 0, while for the
envelope we have

    ρ(R)  ≈ρ(∞),     ρ'(R) ≈ 0,      ρ”(R)≈ 0.

Substituting  these asymptotic values in Eq. (<ref>)
we obtain

    -U(∞)ρ^2(∞) = q^2,

and using U(∞)=-k^2, we find the parameter q,

    q=kρ(∞).

Conversely, if one prefers to choose a certain value for q, the
equation above yields ρ(∞)=q/k.  However, as shown above,
the normalization constant ρ(∞)
is irrelevant; indeed, we have

    θ'(R)=q/ρ(R) = kρ(∞)/ρ(R),

which ensures

    θ'(R) ≈ k,      when R→∞,

and thus the  phase function suitable for scattering problems, as
defined in Sec. <ref>, has the desired  behavior:

    θ(R) ≈ kR,       when R→∞.










        