{
"cells": [
{
"cell_type": "markdown",
"id": "decimal-appliance",
"metadata": {},
"source": [
"# Radial Schrödinger equation with local potential"
]
},
{
"cell_type": "markdown",
"id": "annual-quick",
"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:* Sebastian König\n",
"\n",
"*Last update:* Aug 11, 2021\n",
"\n",
"---"
]
},
{
"cell_type": "markdown",
"id": "actual-superintendent",
"metadata": {},
"source": [
"We consider here the radial Schrödinger equation written as\n",
"\n",
"\\begin{equation}\n",
" \\left[\\frac{\\mathrm{d}^2}{\\mathrm{d}r^2} - \\frac{l(l+1)}{r^2} - 2\\mu V(r) + k^2\\right] u(r) = 0 \\,.\n",
"\\end{equation}"
]
},
{
"cell_type": "markdown",
"id": "industrial-barcelona",
"metadata": {},
"source": [
"Like we did for the Lippmann-Schwinger equation in momentum space, we start by defining our system. The default sets $m = 2\\mu = 1$ and also uses $\\hbar = 1$. For the potential, we also pick the same as in the Lippmann-Schwinger example, an attractive Gaussian with $V_0=-4.0$ and $R=2.0$."
]
},
{
"cell_type": "code",
"execution_count": 1,
"id": "suffering-knife",
"metadata": {},
"outputs": [],
"source": [
"from lib.system import *\n",
"from lib.potential import *\n",
"\n",
"sys = System()\n",
"V = V_Gauss(sys, -4.0, 2.0)"
]
},
{
"cell_type": "markdown",
"id": "emerging-minute",
"metadata": {},
"source": [
"To solve the **ordinary differential equation (ODE)**, we use `odeint` from `scipy.integrate`:"
]
},
{
"cell_type": "code",
"execution_count": 2,
"id": "ready-barrel",
"metadata": {},
"outputs": [],
"source": [
"from scipy.integrate import odeint"
]
},
{
"cell_type": "markdown",
"id": "collectible-metallic",
"metadata": {},
"source": [
"The radial Schrödinger equation is a **second order** ODE. To handle this with the numerical solver, we need to tranform it into a system of first-order ODEs. To that end, let us first write\n",
"\n",
"\\begin{equation}\n",
" u''(r) = \\left[\\frac{l(l+1)}{r^2} + 2\\mu V(r) - k^2\\right] u(r) \\,,\n",
"\\end{equation}\n",
"\n",
"and set $v(r) = u'(r)$. Then we have two coupled first-order equations:\n",
"\n",
"\\begin{align}\n",
" u'(r) &= v(r) \\\\\n",
" v'(r) &= \\left[\\frac{l(l+1)}{r^2} + 2\\mu V(r) - k^2\\right] u(r)\n",
"\\end{align}"
]
},
{
"cell_type": "markdown",
"id": "changed-hardwood",
"metadata": {},
"source": [
"This we can now directly turn into an operator defintion as expected by `odeint`:"
]
},
{
"cell_type": "code",
"execution_count": 3,
"id": "mobile-season",
"metadata": {},
"outputs": [],
"source": [
"def radseq(uv, r, k, l):\n",
" u, v = uv\n",
" return [v, (l * (l + 1) / r**2 + 2.0 * sys.mu * V(r) - k**2) * u]"
]
},
{
"cell_type": "markdown",
"id": "manufactured-january",
"metadata": {},
"source": [
"With this, it is straightforward to solve the ODE. We use here a grid starting essentially at the origin (but avoiding the singularity exactly at $r=0$). As initial conditions we use $u(0)=0$ and $u'(0)=1$, and we solve for $l=0$ and $k = 1.0$:"
]
},
{
"cell_type": "code",
"execution_count": 4,
"id": "popular-bulgarian",
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"\n",
"rs = np.linspace(0.0001, 20.0, 1000)\n",
"\n",
"uv0 = [0.0, 1.0]\n",
"\n",
"l = 0\n",
"k = 1.0\n",
"\n",
"sol = odeint(radseq, uv0, rs, args=(k, l))"
]
},
{
"cell_type": "markdown",
"id": "overhead-riverside",
"metadata": {},
"source": [
"The solution is now a pair of vectors (for $u(r)$ and $v(r)=u'(r)$), stored together in `sol`. We extract the vector for $u(r)$ and plot it, together with the free solution for comparison:"
]
},
{
"cell_type": "code",
"execution_count": 5,
"id": "available-organic",
"metadata": {},
"outputs": [
{
"data": {
"image/svg+xml": [
"\n",
"\n",
"\n",
"\n"
],
"text/plain": [
"