Numerov wavefunction for continuum

[zzpwahaha]

Hi
I am trying to compute the continuum radial wavefuntion of Rb87 with the parameteric model potential https://journals.aps.org/pra/pdf/10.1103/PhysRevA.49.982. I looked up the cpp code in arc_c_extensions.c and have couple of questions:

1. Does the radialWavefunction function support continuum energy? I try with positive stateEnergy and it spit out some resonable thing but see 3.
2. In the code, you re-parametrize the radius with x = sqrt(r). Is that for a better percision? since it is not the regular way we write the radial wavefunction.
3. I try a Numerov method with the radial wf u =r * R, like this with the corresponding paramatric potential taken from you code. Below is my code. What I found is that the solution of the two methods sometimes agree but mostly do not. For example,
   En = 1.5, Rl = np.linspace(1e-7,100,10000)
   does not, but
   En = 0.5, Rl = np.linspace(1e-7,100,10000)
   agree. I am wondering what could cause that, would it be the way of re-parametrization?

        from scipy import *
        import numpy as np
        from scipy import integrate
        from scipy import optimize
        import matplotlib.pyplot as plt
        from scipy import constants as cc
        from arc import *
        atom = Rubidium87()
        def effRadialPotential(l, s, j, r, E):
            # used as f in Numerov(f, x0, dx, dh)
            mu = (atom.mass - cc.m_e) / atom.mass #reduce mass of electron
            return 2 * ( mu * (atom.potential(l,s,j,r) - E) + l * (l + 1)/(2 * r**2) )

        def Numerov(f, x0, dx, dh):
            """Given precomputed function f(x), solves for x(t), which satisfies:
                  x''(t) = f(t) x(t)
            """
            x = np.zeros(len(f))
            x[0] = x0
            x[1] = x0+dh*dx

            h2 = dh**2
            h12 = h2/12.

            w0=x0*(1-h12*f[0])
            w1=x[1]*(1-h12*f[1])
            xi = x[1]
            fi = f[1]
            for i in range(2,len(f)):
                w2 = 2*w1-w0+h2*fi*xi  # here fi=f1
                fi = f[i]  # fi=f2
                xi = w2/(1-h12*fi)
                x[i]=xi
                w0 = w1
                w1 = w2
            return x

        Rl = np.linspace(1e-7,10,10000)
        l=0
        En=0.5
        feff = np.array([effRadialPotential(l,0.5,l+.5, i, En) for i in Rl])
        ur = Numerov(feff,0.0,1e-1,Rl[1]-Rl[0])
        norm = integrate.simps(ur**2,x=Rl)
        ur *= 1/np.sqrt(abs(norm))
        arcr, arcRad = atom.radialWavefunction(l,0.5,l+0.5,En,Rl[0],Rl[-1],Rl[1])

        plt.plot(Rl,ur,'b')
        plt.plot(arcr,arcRad,'r')
        # plt.ylim(-1,1)
        norm

[nikolasibalic]

Hi @zzpwahaha

Thank you for your question.

1. & 3. It was made for integration of bound states. Could it be adopted in similar form for free states? Probably. However, note that for bound states integration starts not from 0 (as in your case) but from large radius R inwards. This is because for bound states large R value of wave function is well defined (we know wave function is 0 at large R), whereas precise potential for small R is not actually that well known. That's why radialWavefunction would normally stop with integration before reaching 0 [at self.alphaC**(1 / 3.0)] for alkali atoms. In any case I guess that this difference in initial condition and direction of integration (you start with 0 wave function, and integrate from core outwards in your example) is making difference between two codes.

2. This is to adjust mesh so that we have denser mesh where radial wave function changes more quickly (close to core, for R->0) and vice versa. One can use also other meshes, for example have a constant step, but then number of integration steps required would be larger and it would take bit more time to get results.

Hope this clarifies all. I would have to think how to extend the existing method correctly for free states.

[zzpwahaha]

Hi Nikola
Thanks for your reply.
I have played around with the free state wf. Here is the way I did it

%%cython -a
import numpy as np
cdef extern from "math.h":
    double fabs(double x)
def Numerovc(double [:] f, double x0, double dx, double dh):
    cdef int divergentPoint = 0
    cdef double maxValue = 1e16
    cdef double [:] x = np.zeros(len(f))
    x[0]=x0
    x[1]=x0 + dh*dx
    cdef double h2 = dh*dh
    cdef double h12 = h2/12.
    
    cdef double w0 = x[0]*(1-h12*f[0])
    cdef double w1 = x[1]*(1-h12*f[1])
    cdef double w2
    cdef double xi = x[1]
    cdef double fi = f[1]
    for i in range(2,len(f)):
        w2 = 2*w1 - w0 + h2*fi*xi  # here fi=f1
        fi = f[i]                     # fi=f2
        xi = w2/(1-h12*fi)
        x[i] = xi
        w0 = w1
        w1 = w2
        # check convergence to avoid instability near origin 
        
    return x

def getRadialWf_Contin(Enconti,Rlsqrt,l,j):
    feffsqrt = - np.array([effRadialPotentialsqrt(l,0.5,j, i, Enconti) for i in Rlsqrt])
    ursqrt = Numerovc(feffsqrt,
                      0.0,
                      -1e-15,
                      Rlsqrt[1]-Rlsqrt[0]) * np.sqrt(Rlsqrt)
    normsqrt = integrate.simps(ursqrt**2,x=Rlsqrt**2)
    return ursqrt*1/np.sqrt(abs(normsqrt))

def effRadialPotentialsqrt(l, s, j, x, stateEnergy):
    r = x * x
    mu = (atom.mass - cc.m_e) / atom.mass
    return -3. / (4. * r) + 4. * r * (
        2. * mu * (stateEnergy - atom.potential(l, s, j, r))
        - l * (l + 1) / (r**2)
        )

where the method is the same as in ARC in the sense that we use the same sqrt(r) discretization as well as the parametric model potential. The start point of integration is at origin. The problem with that is in some case, the wavefunction's amplitude is too large and the normalization returns a Nan.

Hope that will help you
Best,
Zhenpu