Numerov wavefunction for continuum about arc-alkali-rydberg-calculator

[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

[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.