The relativistic J-matrix method - numerical example
Version: 2.0 / 24.01.2007 / R
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
Based on article of Pawel Horodecki, Phys. Rev. A 62 (2000) 052716

Collect input data

Properties of the projectile

l := 1                                      (* angular momentum of the incident electron,  
          should be non - negative integer *)
κ := 1                                        
ε:=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 := 
      5\)\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\( (*\ truncating\ parameter\ for\ the\ scattering\ potential, \ 
      should\ be\ integer\ and\ positive\ *) \)\)\)

Other parameters

(Subscript[v, light]):=137.0359996287515                  (* speed of light (in atomic units) *)

See also square-well parameters section below

Initializations

Energy rescaling and calculating the wave number

\(\(\(ϵ := \@\(ε\/\(ε + 2\ v\_light\^2\)\)

    k := \@\(ε \((ε + 2\ v\_light\^2)\)\)\/v\_light\)\(\ \ \ \)\)\)

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]]\)

\(ψ[n_, r_] := κ\/r\ ϕ[n, r] + 
        If[basis == "laguerre", 
		Exp[\(-\(\(λ\ r\)\/2\)\)] \((\ λ\^\(l + 1\)\ \((l + 1)\)\ r\^l\ \ LaguerreL[n, 
                    2\ l + 1, λ\ r] - \(\((λ\ r)\)\^\(l + 
                        1\)\) λ\/2\ LaguerreL[n, 
                    2\ l + 1, λ\ r] - \ 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \((λ\ r)\)\^\(l + 1\)\ \ \((LaguerreL[n, 2\ l + 2, λ\ r] - 
                      LaguerreL[n, 
                        2\ l + 
                          1, λ\ r])\))\), 
		Exp[\(-\(\(λ\^2\ r\^2\)\/2\)\)] \((\ λ\^\(l + 1\)\ \((l + 
                      1)\)\ r\^l\ \ LaguerreL[n, 
                    l + 1\/2, λ\^2\ r\^2] - \(\((λ\ r)\)\^\(l + 1\)\) λ\^2\ r\ LaguerreL[n, 
                    l + 1\/2, λ\^2\ r\^2] - 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \((λ\ r)\)\^\(l + 1\)\ \ \((LaguerreL[n, l + 3\/2, λ\^2\ r\^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

Alias["δ", KroneckerDelta];

\(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]

Jm1[n_, 
        m_] := -ε fifi[(Subscript[N, trunc]) - n - 1, 
          (Subscript[N, trunc]) - m - 1]
Jm2[n_, m_] := (Subscript[v, light]) psipsi[(Subscript[N, trunc]) - n - 1, m]
Jm3[n_, 
        m_] := -(ε + 2 ((Subscript[v, light])^2)) psipsi[n, 
          m]

Relativistic J-matrix

J1 = Array[Jm1, {(Subscript[N, trunc]), (Subscript[N, trunc])}, 0];
J2 = Array[Jm2, {(Subscript[N, trunc]), (Subscript[N, trunc])}, 0];
J3 = Array[Jm3, {(Subscript[N, trunc]), (Subscript[N, trunc])}, 0];
J4 = Transpose[Array[Jm2, {(Subscript[N, trunc]), (Subscript[N, trunc])}, 0]];

<<  LinearAlgebra`MatrixManipulation`

Jup = AppendRows[J1,  J2];

Jdown = AppendRows[J4,  J3];

J = AppendColumns[Jup,  Jdown];  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]) ,  
        If[n >= (Subscript[N, trunc ]) &&  m >= (Subscript[N, trunc]), 
          NIntegrate[ψ[n - (Subscript[N, trunc]), r] V[r] ψ[m - (Subscript[N, trunc]), 
                r], {r, a, b}], 0],  
        NIntegrate[ϕ[(Subscript[N, trunc]) - n - 1, r] V[r] ϕ[
              (Subscript[N, trunc]) - m - 1, r], {r, a, b}]]

(Subscript[V, N]) = Array[(Subscript[V, trunc]), {2 (Subscript[N, trunc]), 2 (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]);  
    Chop[MatrixForm[(Subscript[A, diag])], (10^(-7))]

Compute matrix approximating Green function

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

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

Compute the numerical phase shift

Main formula

(Subscript[g, N]) = g[[1, 1]];
(Subscript[J, N]) = Jm[(Subscript[N, trunc]), (Subscript[N, trunc]) - 1];
tmp = ((2 ϵ (Subscript[v, light]))/k);

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

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

The second, equivalent formula

(Subscript[g, N]) = g[[1, 2 (Subscript[N, trunc])]];
tmp = 2 (Subscript[v, light]);

tan(Subscript[δ, N]) = -((s[(Subscript[N, trunc]) - 1] + 
                tmp (Subscript[g, N]) (Subscript[J, N]) s[(Subscript[N, trunc])])/(c[(Subscript[N, trunc]) - 1] + 
                tmp (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

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

Energy rescaling and calculating wave number within the "potential area"

(Subscript[ϵ, prim]) := (Sqrt[((ε - (Subscript[V, 0]))/((ε - (Subscript[V, 0])) + 2 ((Subscript[v, light])^2)))])
(Subscript[k, prim]) := ((Sqrt[(ε - (Subscript[V, 0])) ((ε - (Subscript[V, 0])) + 2 ((Subscript[v, light])^2))])/(Subscript[v, light]))

Form the matrices

\(\(\(X[r_, k_, ϵ_] := 
      If[κ < 
          0, \ {\ {j[l, 
              k\ r], \(-n[l, k\ r]\)}, {\(-ϵ\)\ j[l + 1, 
                k\ r], ϵ\ n[l + 1, k\ r]}}, 
		{\ {j[l, 
              k\ r], \(-n[l, k\ r]\)}, {ϵ\ j[l - 1, 
                k\ r], \(-ϵ\)\ n[l - 1, k\ r]}}]\)\(
\)
    \)		\)
Y[r_, k_, ϵ_] := 
      If[κ < 
          0, { {ϵ n[l + 1, k r], 
            n[l, k r]}, {ϵ j[l + 1, k r], 
            j[l, k r]}}, { {-ϵ n[l - 1, k r], 
            n[l, k r]}, {-ϵ j[l - 1, k r], 
            j[l, k r]}}]

Y1 = Y[b, k, ϵ];

X1 = X[b, (Subscript[k, prim]), (Subscript[ϵ, prim])];

Y2 = Y[a, (Subscript[k, prim]), (Subscript[ϵ, prim])];

X2 = X[a, k, ϵ];

XY = Y1 . X1 . Y2 . X2 . {1, 0};

Tangent to the analytical phase shift, and the phase shift itself

tan(Subscript[δ, acc]) = XY[[2]]/XY[[1]];

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