The non-relativistic J-matrix method - numerical example
Version: 2.0 / 24.01.2007 / NR
Author: Pawel Syty, Gdansk University of Technology, Gdansk, Poland, e-mail: sylas@mif.pg.gda.pl
Fortran 95 code is also available, see http://jmatrix.sylas.info

Collect input data

Properties of the projectile

l := 1                                      (* angular momentum of the incident electron, should be non - negative integer *)
ε:=3.0                                   (* 
        kinetic energy of the electron (in atomic units),  
        should be positive *)

Basis set, its parameters and size

\(\(\(basis := "gauss"\ \ \ \ \ \ \ \ \  (*\ 
      basis\ set, \ "laguerre"\ or\ "gauss"\ *) 

    λ := 
      1.0\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  (*\ scaling\ parameter, \ 
      should\ not\ be\ changed\ in\ the\ present\ version\ *) 

    N\_trunc := 
      10\)\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\( (*\ truncating\ parameter\ for\ the\ scattering\ potential, \ 
      should\ be\ integer\ and\ positive\ *) \)\)\)

See also square-well parameters section below

Initializations

Wave number

k := (Sqrt[2 ε])

Basis functions

\(ϕ[n_, r_] := 
      If[basis == "laguerre", 
		\((λ\ r)\)\^\(l + 1\)\ Exp[\(-\(\(λ\ r\)\/2\)\)]\ LaguerreL[n, 
            2\ l + 1, λ\ r], 
		\(\((λ\ r)\)\^\(l + 1\)\) 
          Exp[\(-\(\(λ\^2\ r\^2\)\/2\)\)]\ LaguerreL[n, 
            l + 1\/2, λ\^2\ r\^2]]\)

Expansion coefficients of the sine- and cosine-like solutions

sint := (((k/λ))/(((k^2)/(λ^2)) + (1/4)));       cost := ((((k^2)/(λ^2)) - (1/4))/(((k^2)/(λ^2)) + (1/4)));        sin2t2 := ((λ^2)/(4 (k^2) + (λ^2)));     η := 
      (k/λ)

\(s[n_] := 
      If[basis == "laguerre", 
		\(2\^l\ Gamma[l + 1]\ \(n!\)\ sint\^\(l + 1\)\)\/Gamma[n + 2\ l + 2]\ GegenbauerC[n, l + 1, 
            cost], 
		\(\@\(2\ π\)\ \(n!\)\ \((\(-1\))\)\^n\)\/Gamma[n + l + 3\/2]\ Exp[\(-η\^2\)\/2]\ η\^\(l + 1\)\ LaguerreL[n, 
            l + 1\/2, η\^2]]\)

\(c[n_] := 
      If[basis == "laguerre", 
		\(\(-2\^l\)\ Gamma[l + 1\/2]\ \(n!\)\)\/\(\@π\ Gamma[n + 2\ l + 2]\ sint\^l\)\ Hypergeometric2F1[\(-n\) - 2\ l - 1, n + 1, 1\/2 - l, 
            sin2t2], 
		\(\(\@\(2\/π\)\) Gamma[l + 1\/2]\ \((\(-1\))\)\^n\ \(n!\)\)\/\(\(\ \)\(Gamma[n + \ l + 3\/2]\)\)\ Exp[\(-η\^2\)\/2]\ η\^\(-l\)\ \ Hypergeometric1F1[\(-n\) - l - 1\/2, 
            1\/2 - l, η\^2]]\)

Form the J-matrix

J-matrix elements

δ[a_, b_]:=If[a == b, 1.0, 0.0]

\(fifi[n_, m_] := 
      If[basis == "laguerre", 
		Gamma[n + 2\ l + 2]\/\(λ\ \(n!\)\)\ \((2\ \((n + l + 1)\)\ δ[m, n] - 
              n\ δ[m, n - 1] - \((n + 2\ l + 2)\)\ δ[m, 
                  n + 1])\), 
		Gamma[n + l + 3\/2]\/\(2\ \(n!\)\)\ δ[m, n]]\)

\(psipsi[n_, m_] := 
      If[basis == "laguerre", 
		Gamma[n + 2\ l + 2]\/\(4\ \(n!\)\)\ \((2\ \((n + l + 1)\)\ \((2\ λ - 1)\)\ δ[m, n] + 
              n\ δ[m, n - 1] + \((n + 2\ l + 2)\)\ δ[m, 
                  n + 1])\), 
		\(\(λ\^2\ Gamma[
                  n + l + 3\/2]\)\/\(2\ \(n!\)\)\) \((\ \((2\ n + l + 
                    3\/2)\)\ \ δ[m, n] + 
              n\ δ[m, n - 1] + \((n + l + 3\/2)\)\ δ[m, 
                  n + 1])\)]\)

Jm[n_, m_] := (1/2) psipsi[n, m] - ((k^2)/2) fifi[n, m]

Non-relativistic J-matrix

J = Array[Jm, {(Subscript[N, trunc]), (Subscript[N, trunc])}, 0];  MatrixForm[J]

Scattering potential

Square-well parameters

\(\(\(a := 
      0.8\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  (*\ 
      left\ bound\ *) 

    b := 1.0\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \  (*\ 
      right\ bound\ *) 

    V\_0 := \(-1.0\)\)\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\( (*\ depth\ *) \)\)\)

V[r_] := If[r > a  &&  r < b, (Subscript[V, 0]), 0]

Plot[V[r], {r, a - 1, b + 1}, 
      PlotStyle -> {Thickness[ .01], Hue[ .7]}]

Truncate the potential in the selected basis

(Subscript[V, trunc])[n_, m_] := 
      If[n >= (Subscript[N, trunc])  ||  m >= (Subscript[N, trunc]) ,  0, 
        NIntegrate[ϕ[n, r] V[r] ϕ[m, r], {r, a, b}]]

(Subscript[V, N]) = Array[(Subscript[V, trunc]), {(Subscript[N, trunc]), (Subscript[N, trunc])}, 0];  
    MatrixForm[(Subscript[V, N])]

The matrix A=J+V

Form the matrix

A = J + (Subscript[V, N]);  MatrixForm[A]

Diagonalize the matrix

Γ = Eigenvectors[A];

(Subscript[Γ, tr]) = Transpose[Γ];

(Subscript[A, diag]) = Γ . A . (Subscript[Γ, tr]);  
    MatrixForm[(Subscript[A, diag])] // Chop

Compute matrix approximating Green function

gm[n_, m_] := 
      Sum[((Subscript[Γ, tr])[[n, j]] (Subscript[Γ, tr])[[m, 
                  j]])/(Subscript[A, diag])[[j, j]], {j, 1, (Subscript[N, trunc])}]

g = Array[gm, {(Subscript[N, trunc]), (Subscript[N, trunc])}];

Compute the numerical phase shift

(Subscript[g, N]) = g[[(Subscript[N, trunc]), (Subscript[N, trunc])]];
(Subscript[J, N]) = Jm[(Subscript[N, trunc]), (Subscript[N, trunc]) - 1];

tan(Subscript[δ, N]) = -((s[(Subscript[N, trunc]) - 1] + 
                (Subscript[g, N]) (Subscript[J, N]) s[(Subscript[N, trunc])])/(c[(Subscript[N, trunc]) - 1] + 
                (Subscript[g, N]) (Subscript[J, N]) c[(Subscript[N, trunc])]));

(Subscript[δ, N]) = ArcTan[tan(Subscript[δ, N])]

Analytical formula for square-well

Spherical Bessel and Neumann functions and their derivatives

j[v_, z_] :=  (Sqrt[(π/(2 z))]) BesselJ[v + (1/2), z]
n[v_, z_] := (Sqrt[(π/(2 z))]) BesselY[v + (1/2), z]

jp[v_, 
        z_] := (Sqrt[(π/8)]) (z^(-(3/2))) (z BesselJ[-(1/2) + v, z] - BesselJ[(1/2) + v, z] - z BesselJ[(3/2) + v, z])

np[v_, 
        z_] := (Sqrt[(π/8)]) (z^(-(3/2))) (z BesselY[-(1/2) + v, z] - BesselY[(1/2) + v, z] - z BesselY[(3/2) + v, z])

Calculate wave number within the "potential area"

(Subscript[k, 1]) = (Sqrt[2 (ε - (Subscript[V, 0]))]);

Tangent to the analytical phase shift

\(\(tanδ\_acc\  = \ 
\((k\ jp[l, k\ b]\ n[l, 
                  k\_1\ b]\ \((k\_1\ j[l, k\ a]\ jp[l, k\_1\ a] - 
                    k\ j[l, k\_1\ a]\ jp[l, k\ a])\) + \ 
              k\ j[l, k\_1\ b]\ jp[l, 
                  k\ b]\ \((jp[l, k\ a]\ k\ n[l, k\_1\ a] - 
                    k\_1\ j[l, k\ a]\ np[l, k\_1\ a])\) + \ 
              k\_1\ j[l, 
                  k\ b]\ \((j[l, k\ a]\ jp[l, k\_1\ b]\ k\_1\ np[l, 
                        k\_1\ a] - 
                    k\ jp[l, k\_1\ b]\ jp[l, k\ a]\ n[l, \ 
                        k\_1\ a] + 
k\ j[l, k\_1\ a]\ jp[l, 
                        k\ a]\ np[l, k\_1\ b] - 
                    k\_1\ j[l, k\ a]\ jp[l, k\_1\ a]\ np[l, 
                        k\_1\ b])\))\)/
\((k\ np[l, 
                  k\ b]\ \((k\ j[l, k\_1\ b]\ jp[l, k\ a]\ n[l, k\_1\ a] - 
                    k\ j[l, k\_1\ a]\ jp[l, \ k\ a]\ n[l, 
                        k\_1\ b] + 
\ 
                    k\_1\ j[l, k\ a]\ jp[l, k\_1\ a]\ n[l, k\_1\ b] - 
                    k\_1\ j[l, k\_1\ b]*
                      j[l, k\ a]\ np[l, k\_1\ a])\) - 
\ 
              k\_1\ n[l, 
                  k\ b]\ \((k\ jp[l, k\_1\ b]\ jp[l, k\ a]\ n[l, k\_1\ a] - 
                    k\_1\ j[l, k\ a]\ jp[l, k\_1\ b]\ np[l, 
                        k\_1\ a] - 
\ 
                    k\ j[l, k\_1\ a]\ jp[l, k\ a]\ np[l, k\_1\ b] + 
                    k\_1\ j[l, k\ a]\ jp[l, k\_1\ a]\ np[l, 
                        k\_1\ b])\))\);\)\)

Analytical phase shift

(Subscript[δ, acc]) = ArcTan[tan(Subscript[δ, acc])]