{
"cells": [
{
"cell_type": "markdown",
"metadata": {
"id": "authorized-treasurer"
},
"source": [
"# Jost function method I: bound state"
]
},
{
"cell_type": "markdown",
"id": "valid-gothic",
"metadata": {},
"source": [
"---\n",
"This code is made available as part of the **FRIB-TA Summer School: A Practical Walk Through Formal Scattering Theory: Connecting Bound States, Resonances, and Scattering States in Exotic Nuclei and Beyond**, held virtually\n",
"August 4-6, 2021.\n",
"\n",
"https://fribtascattering.github.io/\n",
"\n",
"*Organizers/Lecturers:*\n",
"- Kévin Fossez (MSU/ANL)\n",
"- Sebastian König (NCSU)\n",
"- Heiko Hergert (MSU)\n",
"\n",
"*Author:* Kévin Fossez\n",
"\n",
"*Last update:* Aug 25, 2021\n",
"\n",
"---"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "L1YDR7xm7P2f"
},
"source": [
"The Jost function method will be introduced on the third day, which means that you will start implementing it without really knowing how it works. However, you already have all the tools to do it and understand the general idea.\n",
"\n",
"In short, given the momentum $k$ of any type of state (bound, scattering, resonant), it is a method to obtain the wave function using properties of the Jost functions. This is achieved by replacing the Schrödinger equation (2nd order ODE) by two first order ODEs for the \"pseudo-Jost functions\"."
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "BRABlb9ipf9o"
},
"source": [
"## Basic setup"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "i9xe0inXp3Om",
"outputId": "910fbe7f-57a6-4adf-aac5-2ec057322311"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"zsh:1: command not found: gdown\n",
"unzip: cannot find or open lib.zip, lib.zip.zip or lib.zip.ZIP.\n"
]
},
{
"name": "stderr",
"output_type": "stream",
"text": [
"/var/folders/ld/0k7b046s3t7dpld1w1b2n_zw0000gn/T/ipykernel_22011/4058438398.py:8: DeprecationWarning: `set_matplotlib_formats` is deprecated since IPython 7.23, directly use `matplotlib_inline.backend_inline.set_matplotlib_formats()`\n",
" set_matplotlib_formats('svg')\n"
]
}
],
"source": [
"# Download and install library helpers\n",
"!gdown --id 1sZT60tdWLtmBCFKzs870UhiMurVdSPRk\n",
"!unzip -o lib.zip\n",
"\n",
"# Set up pretty \n",
"from IPython.display import set_matplotlib_formats\n",
"%matplotlib inline\n",
"set_matplotlib_formats('svg')"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "aBftxc_esfTC"
},
"source": [
"## Riccati recursion"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "8b20fca9"
},
"source": [
"To use the Jost function method to obtain the radial wave functions of bound, scattering, and resonance states, we will need to do an analytical continuation of the spherical Bessel functions introduced in the lecture notes.\n",
"\n",
"Of course, libraries already exist, but it is instructive to know how it can be done (at least one way). Here, the idea is to use a recursion formula from the NIST handbook (10.51(i)) which gives the spherical Bessel functions at any order from the functions at the lowest orders, and the fact that the closed forms of the lowest orders are simple and can be analytically continued into the complex plane without any special care for $z \\neq 0$.\n",
"\n",
"The recursion formula is given by:\n",
"\n",
"\\begin{equation}\n",
"f_{n+1}(z) = \\frac{2n+1}{z} f_{n}(z) - f_{n-1}(z)\n",
"\\end{equation}\n",
"\n",
"This formula is actually not just for the spherical Bessel functions, so we will name the corresponding function accordingly.\n",
"\n",
"**- Implement the function using the recursion formula**"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"id": "0aeca466"
},
"outputs": [],
"source": [
"def ji_or_yi_or_h1_or_h2_recursion (n, z, f0, f1):\n",
" fi = 0.0\n",
" for i in range (2, n):\n",
" fi = (2.0*i-1.0)/z*f1 - f0\n",
"\n",
" f0 = f1\n",
" f1 = fi\n",
"\n",
" return fi"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "21ba0294"
},
"source": [
"Here, $n$ is the order, $z$ the argument, and $f_{0,1}$ are the values of the function in $z$ at the order 0 and 1 respectivelly. Tips: replace $n+1$ in the formula by $i$ in the code, so that $2n+1 = 2(i-1)+1 = 2i-1$. This is convenient to have $f_0$ on the right side of the formula at the first order. This will be, of course, a naive implementation and visible numerical differences will appear when compared against a library for $n > 6$.\n",
"\n",
"Then, the next step is to implement the actual spherical Bessel functions using the generic recursion formula above. We know that:\n",
"\n",
"\\begin{align}\n",
"& j_0(z) = \\frac{\\sin(z)}{z} \\\\\n",
"& j_1(z) = \\frac{\\sin(z)}{z^2} - \\frac{\\cos(z)}{z}\n",
"\\end{align}\n",
"and:\n",
"\n",
"\\begin{align}\n",
"& y_0(z) = -\\frac{\\cos(z)}{z} \\\\\n",
"& y_1(z) = -\\frac{\\cos(z)}{z^2} - \\frac{\\sin(z)}{z}\n",
"\\end{align}\n",
"\n",
"**- Implement the functions using the previous one**"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {
"id": "39ec1833"
},
"outputs": [],
"source": [
"import cmath\n",
"\n",
"def jn (l, z):\n",
" if z == 0.0:\n",
" return 0\n",
" print(\"jn: z must be nonzero\")\n",
" exit()\n",
"\n",
" f0 = cmath.sin(z)/z\n",
" if l == 0:\n",
" return f0\n",
"\n",
" f1 = cmath.sin(z)/(z*z) - cmath.cos(z)/z\n",
" if l == 1:\n",
" return f1\n",
"\n",
" return ji_or_yi_or_h1_or_h2_recursion(l, z, f0, f1)\n",
"\n",
"\n",
"def yn (l, z):\n",
" if z == 0.0:\n",
" return 0\n",
" print(\"yn: z must be nonzero\")\n",
" exit()\n",
"\n",
" f0 = - cmath.cos(z)/z\n",
" if l == 0:\n",
" return f0\n",
"\n",
" f1 = - cmath.cos(z)/(z*z) - cmath.sin(z)/z\n",
" if l == 1:\n",
" return f1\n",
"\n",
" return ji_or_yi_or_h1_or_h2_recursion(l, z, f0, f1)"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "d0654e34"
},
"source": [
"In addition, we can also define the incoming (-) and outgoing (+) Hankel functions in the special case where there is no Coulomb potential ($\\eta=0$):\n",
"\n",
"\\begin{equation}\n",
"H_{l}^\\pm(z) = z [ y_l(z) \\pm i j_l(z) ]\n",
"\\end{equation}"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {
"id": "cb985ce9"
},
"outputs": [],
"source": [
"def Hp_Hm (sign, l, z):\n",
" return z*(yn(l,z)+sign*1.0j*jn(l,z))"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "c58f4353"
},
"source": [
"**- Check the validity of what you implemented using the following python library**"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {
"id": "de967340"
},
"outputs": [],
"source": [
"import scipy.special as sp\n",
"l=0\n",
"kr=1.0\n",
"jl = sp.spherical_jn(l,kr)\n",
"yl = sp.spherical_yn(l,kr)"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "_zrthm9BsxPt"
},
"source": [
"# Potential"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "eb6a1326"
},
"source": [
"To solve a real nuclear physics case, we can also define a Woods-Saxon potential as follows:\n",
"\n",
"\\begin{equation}\n",
"V(r) = - V_0 f(r) - 4 \\vec{l}.\\vec{s} V_{so} \\frac{1}{r} \\frac{df(r)}{dr}\n",
"\\end{equation}\n",
"with:\n",
"\n",
"\\begin{equation}\n",
"f(r) = - \\frac{1}{ 1 + e^{\\frac{r-R_0}{d}} }\n",
"\\end{equation}\n",
"\n",
"We can also use the fact that (show it):\n",
"\n",
"\\begin{equation}\n",
"2 \\vec{l}.\\vec{s} = j(j+1) - l(l+1) - \\frac{3}{4}\n",
"\\end{equation}"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {
"id": "554f6a7e"
},
"outputs": [],
"source": [
"def pot_WS(r):\n",
"\n",
" d = 0.618\n",
" R0 = 2.162\n",
" Vo = 41.77\n",
" Vso = 6.991\n",
"\n",
" exp_r_minus_R0_over_d = np.exp ((r - R0)/d)\n",
" f = -1.0/(1.0 + exp_r_minus_R0_over_d)\n",
" Vnuclear = Vo * f;\n",
"\n",
" df = exp_r_minus_R0_over_d*f*f/d\n",
" two_l_scalar_s = j*(j+1) - l*(l+1) - 0.75\n",
"\n",
" Vspin_orbit = 0.0\n",
" if ((two_l_scalar_s != 0) and (Vso != 0.0) and r != 0.0):\n",
" Vspin_orbit = -2.0*two_l_scalar_s*Vso*df/r\n",
"\n",
" return Vnuclear + Vspin_orbit;"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "3OADwzZ_s93z"
},
"source": [
"# Wave function"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "a86ead50"
},
"source": [
"Now, we can implement the Jost function method naively, at first, and then improve it. We recall that the wave function can be obtained from the Jost functions as ($\\eta=0$ here):\n",
"\n",
"\\begin{equation}\n",
"u_l(r,k) = \\frac{i}{2} \\left( F_{l,\\eta}^+(k,r) H_{l,\\eta}^-(k,r) - F_{l,\\eta}^-(k,r) H_{l,\\eta}^+(k,r) \\right)\n",
"\\end{equation}\n",
"\n",
"**- Implement a function to return the wave function given the values of $k$, $r$, and $F^\\pm(r)$**"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {
"id": "434cfbe2"
},
"outputs": [],
"source": [
"def wf_u (k, r, Fp, Fm, l):\n",
" if (r == 0.0): \n",
" return 0.0\n",
"\n",
" return 1.0j*0.5*(Fp*Hp_Hm(-1.0,l, k*r) - Fm*Hp_Hm(+1.0,l, k*r))"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "kgAaJnRTtCxQ"
},
"source": [
"# Derivatives of the Jost functions"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "1ef1b8d8"
},
"source": [
"The critical point is then to calculate the derivatives of the Jost functions using the formula:\n",
"\n",
"\\begin{equation}\n",
"\\frac{\\partial}{\\partial r}F_{l,\\eta}^\\pm(k,r) = \\pm\\frac{1}{ik} \\frac{2m}{\\hbar^2} V(r)(\\mp i) H_{l,\\eta}^\\pm(k,r) u_{l}(r,k)\n",
"\\end{equation}\n",
"\n",
"**- Using the WS potential and the wave function, implement the derivative of the Jost functions and return [dFp,dFm]**"
]
},
{
"cell_type": "code",
"execution_count": 8,
"metadata": {
"id": "084448bb"
},
"outputs": [],
"source": [
"def dF(F, r, k, two_amu_over_hbar2):\n",
" Fp = F[0]\n",
" Fm = F[1]\n",
"\n",
" u = wf_u(k,r,Fp,Fm, l)\n",
"\n",
" dFp = 0.0\n",
" dFm = 0.0\n",
"\n",
" if r==0.0:\n",
" return [0.0,0.0]\n",
"\n",
" Hp = Hp_Hm (+1.0,l, k*r)\n",
" Hm = Hp_Hm (-1.0,l, k*r)\n",
" dFp = +(1.0/(1.0j*k))*two_amu_over_hbar2*pot_WS(r)*(-1.0j*Hp)*u\n",
" dFm = -(1.0/(1.0j*k))*two_amu_over_hbar2*pot_WS(r)*(+1.0j*Hm)*u\n",
"\n",
" return [dFp,dFm]"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "ALUBXsrwtHPu"
},
"source": [
"# Numerical integration"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "c0fdca24"
},
"source": [
"We are now almost ready to solve the Schrödinger equation and obtain the wave function by integrating the ODEs for the Jost functions. Due to the way the standard ODE solver used in the example presented before (simple 2nd order ODE) is implemented, it cannot immediately handle complex numbers and pass arguments trivially. To use it, we would need to implement an interface, but this would be a little convoluted for this exercise, even though it is fairly standard. We will thus use a different package.\n",
"\n",
"In the following function, we initialize the Jost functions as they should, and then solve the ODEs using the derivatives defined above. Then, we extract the Jost functions obtained and calculate the wave function.\n",
"\n",
"We remind that the Jost function in zero must equal 1, or `1.0+0.0j` to avoid issues with the ODE solver.\n",
"\n",
"**- Implement the function to calculate the wave function by solving the ODEs for the Jost functions first**"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {
"colab": {
"base_uri": "https://localhost:8080/"
},
"id": "mMhQ5JTszcbA",
"outputId": "6b837a06-ef27-4ffa-9bbe-21cc138a6005"
},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"\u001b[33mDEPRECATION: Configuring installation scheme with distutils config files is deprecated and will no longer work in the near future. If you are using a Homebrew or Linuxbrew Python, please see discussion at https://github.com/Homebrew/homebrew-core/issues/76621\u001b[0m\n",
"Requirement already satisfied: odeintw in /usr/local/lib/python3.9/site-packages (0.1.1)\n",
"Requirement already satisfied: scipy in /usr/local/lib/python3.9/site-packages (from odeintw) (1.6.0)\n",
"Requirement already satisfied: numpy>=1.16.5 in /usr/local/lib/python3.9/site-packages (from scipy->odeintw) (1.19.4)\n"
]
}
],
"source": [
"!pip3 install odeintw\n",
"\n",
"from odeintw import odeintw"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {
"id": "b164dd9d"
},
"outputs": [],
"source": [
"def integrate_wf (r_range, k, l, two_amu_over_hbar2):\n",
"\n",
" u_range=np.array([])\n",
" pot_range = np.array([])\n",
"\n",
" # initial values of the r-dependent Jost function at r=0\n",
" # they must be declared as complex (even if real) for the ODE solver\n",
" #\n",
" Fp0 = 1.0+0.0j\n",
" Fm0 = 1.0+0.0j\n",
" F0 = np.array([Fp0, Fm0])\n",
"\n",
" F = odeintw(dF,F0,r_range, args=(k, two_amu_over_hbar2,))\n",
"\n",
" # extract the outgoing and incoming r-dependent Jost functions\n",
" Fp_range = np.array([F[i][0] for i in range(Nr)])\n",
" Fm_range = np.array([F[i][1] for i in range(Nr)])\n",
"\n",
"\n",
" # compute the wave function using the Jost functions obtained\n",
" for i in range(Nr):\n",
" Fp = Fp_range[i]\n",
" Fm = Fm_range[i]\n",
" r = r_range[i]\n",
" u = wf_u(k,r,Fp,Fm, l)\n",
" u_range = np.append(u_range, u)\n",
"\n",
" pot_range = np.append(pot_range, pot_WS(r))\n",
"\n",
" return u_range, pot_range, Fp_range, Fm_range\n"
]
},
{
"cell_type": "markdown",
"metadata": {
"id": "884f2a8a"
},
"source": [
"and we run the program for a given value of the momentum corresponding to a bound state (given here, but how could you obtain it?):"
]
},
{
"cell_type": "code",
"execution_count": 11,
"metadata": {
"id": "30f84d9f",
"outputId": "2f829b4c-7ef0-47e8-98e3-af1330514f69"
},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
"