Cooper-Zare-Theory - Photoionization simulated in Python

My concrete problem is that I want to calculate cross-sections for photoionization transitions in argon in Python. For example from the excited state 3S^23p^56d into the continuum. For this, I found Cooper-Zare theory in several sources, including "Introduction to Photoelectron Angular Distributions". The formula for the cross section is, according to the above-mentioned book(p.87): (i am doing some reverse engineering here) $\sigma_{pd} = \frac{16\pi^2}{3(2l_0 + 1)} [l_0\Re^2_{l_0-1}+(l_0 + 1)\Re^2_{l_0+1}]$ where $l_0$ is the initial azimuthal quantum number (in the upper case: $l_0=2$ , because of the d-state) and $\Re_{l_0\pm1} = \int_0^{\infty} r^3R_{nl_0}(r) G_{k,l_0\pm1}(r) \, dr$ where $n$ is the principal quantum number (here $n=6$ ). And $G_{kl}(r) = \frac{e^{i\delta_l} \sin(kr - l\pi/2)+ e^{-i(kr + l\pi/2)} e^{-2i\delta_l}\sin(\delta_l)}{kr}$ which originates from $| f \rangle \propto 4\pi\sum_{l, m} i^l e^{-i\delta_l} Y_{lm}^* (\hat{k})Y_{lm}(\hat{r}) G_{kl}(r)$ where $\hat{k}$ defines a unit vector in the direction of the photo-electron momentum vector. And this comes from: $| f \rangle \propto \sum_{l=0}^{\infty} (2l + 1)i^l e^{-i\delta_l} P_l(\cos\theta_s) \left[\frac{e^{i\delta_l} \sin(kr - l\pi/2) \left(\frac{1}{2} + e^{-i kr + l\pi/2}\right) e^{-2i\delta_l} \sin(\delta_l)}{kr}\right]$ which is assumed to be the final state for the dipole approximation $\frac{d\sigma_{pd}}{d\Omega} = |\langle f |O| i \rangle|^2$ with the operator for linearly polarized light: $O = r^4 \sqrt{\frac{4\pi}{3}} Y_{1 0}(\mathbf{r})$ My problem is that I don't know where to get the phase shifts $\delta_l$ of the partial wave, which are used in the formula for $G$ . I would be very grateful for any tips. import scipy.special as spe import scipy.integrate as integrate

#radial wave function def psi_R(r,n=1,l=0): rr=r coeff = np.sqrt((2.0/n)**3 * spe.factorial(n-l-1) /(2.0nspe.factorial(n+l))) laguerre = spe.assoc_laguerre(2.0rr/n,n-l-1,2l+1) return coeff * np.exp(-rr/n) * (2.0*rr/n)**l * laguerre

the radial portion of the final-state wave function divided by kr

def G_kl (r,k,l,d_l): gg=(np.exp(1jd_l)np.sin(kr-lnp.pi/2)+np.exp(-1j*(kr+lnp.pi/2))np.exp(-2jd_l)np.sin(d_l))/(kr) return gg

e radial dipole integral

def R_l0(n,l0,l,k,d_l): integ= integrate.quad(lambda r: np.real(r**3*psi_R(r,n,l0)*G_kl(r,k,l,d_l)), 0, np.inf,limit=100000, epsabs=np.inf) return integ[0]

asymmetry parameter β

def beta(n,l,k,dl_p1,dl_m1): R_p1=R_l0(n,l,l+1,k,dl_p1) if l-1>=0: R_m1=R_l0(n,l,l-1,k,dl_m1) else: R_m1=0 chi= l*(l-1)R_m1**2+(l+1)(l+2)R_p1**2-6l*(l+1)R_m1R_p1np.cos(dl_p1-dl_m1) epsi=(2l+1)(lR_m12+(l+1)*R_p12)

return chi/epsi

def cross_sec(n,l0,k,dl_p1,dl_m1): r_lp1=R_l0(n,l0,l0+1,k,dl_p1) r_lm1=R_l0(n,l0,l0-1,k,dl_m1) return 16np.pi**2/(3(2l0+1))(l0*r_lm12+(l0+1)*r_lp12)