\documentclass[12pt]{article}
\usepackage{epsfig}
%\topmargin=+0.0in
%\oddsidemargin=0.25in
\textwidth=6.0in
%\textheight=8in
%\footskip=6ex
%\footheight=2ex
\begin{document}
\begin{center}
{\bf Lecture 2} \\
Finding roots of a function, $\Lambda$ hypernuclei
\end{center}
We start the course with a standard numerical problem: finding the
roots of a function. This exercise is meant to get you started
with programming in gcc or ROOT, and hopefully also in using
latex. We will apply a method to finding the energy
levels of a $\Lambda$ particle bound in a nucleus. First we
discuss an algorithm, then the physics.\\
\centerline{\bf Finding the roots of a Function}
Often in Physics one has to solve an algebraic equation that does
not have a simple form. One can usually cast the equation into the
form
\begin{equation}
f(x) = 0
\end{equation}
\noindent The goal is to find the value(s) of $x$ that satisfy this
equation, i.e. the roots of $f(x)$. If there is no analytic solution,
then numerical methods can be used to find a solution to the desired
accuracy. In lecture I will describe a method we call the
"Bracket and Half". In the appendix I also describe two more
standard algorithms: the Bisection Method, and Newton's method. In each method,
an algorithm is developed for producing a series of $x's$, $x_i$. The
algorithm is designed such that the series $x_i$ converges to a root
of $f(x)$.\\
\centerline{\bf Bracket and Half Method}
The bracket and half method is best explained by describing the algorithm.
One choses an initial starting value, $x_1$, and an initial step size $\delta$.
The next value of $x$, $x_2 = x_1 + \delta$. If $f(x_2)f(x_1) > 0$, then we
have not passed a root. This is because if the product is greater than
zero, the sign of $f(x_2)$ and $f(x_1)$ is the same. We keep the same value of $\delta$
and continue
with $x_3=x_2+\delta$. We keep increasing $x_i$ by $\delta$
until $f(x_{i+1})f(x_i) < 0$. When this occurs, we have steped over the root,
because now the signs of $f(x_i)$ and $f(x_{i+1})$ are different.
Now we halve the step size and reverse its direction by
changing $\delta$ to $-\delta /2$. The process is continued
until we have reached the desired accuracy.\\
\newpage
Below is an example code for the bracket and half method:\\
\noindent a=starting value \\
del=initial step size \\
for $(i=1; ixtol)$ \\
$\{$ \\
b = a + del; \\
if $(f(b)*f(a)<0)$ then del=-del/2; \\
a = b; \\
$\}$ \\
The bracket and half method is useful for finding the zeros above (or below)
a certain value of $x$. It is also useful in finding the bound state
energy levels of the non-relativistic Schroedinger equation, which we
will consider in the next assignment. For this
application, the algorithm is refered as the "shooting method".\\
The particular equation that we will solve is
\begin{equation}
tan(\sqrt{\alpha (x-1)} \; ) + \sqrt{x-1} = 0
\end{equation}
\noindent where $x>1$. Since $\alpha$ will be greater than zero, the solution
for $x$ will occur when the tan function is negative, i.e. when the argument of the tangent
is in second quadrant. You should start your initial value at $x_1 > 1$. Let me discuss
the physics behind this equation.
\newpage
\centerline{\bf Lambda Hypernuclei}
When particles are confined to a region of space, the energies that they can
have are quantized. For our first example, we will examine the energies
that a lambda ($\Lambda$) particle can have when it is "trapped" inside of
a nucleus. The $\Lambda$ particle is similar to a neutron and proton, and
is a member of their "baryon octet":
\begin{center}
\begin{tabular}{c|c|c|c|c|c}
baryon & mass($MeV/c^2$) & charge & lifetime (sec) & quarks & spin \\
\hline
proton & 938.3 & +e & $\infty$ & uud & 1/2 \\
neutron & 939.7 & 0 & 880 & udd & 1/2 \\
$\Lambda$ & 1115.7 & 0 & $2.63 \times 10^{-10}$ & uds & 1/2 \\
$\Sigma^+$ & 1189.4 & +e & $8.02 \times 10^{-11}$ & uus & 1/2 \\
$\Sigma^0$ & 1192.6 & 0 & $7.4 \times 10^{-20}$ & uds & 1/2 \\
$\Sigma^-$ & 1197.4 & -e & $1.48 \times 10^{-10}$ & dds & 1/2 \\
$\Xi^0$ & 1314.9 & 0 & $2.9 \times 10^{-10}$ & uss & 1/2 \\
$\Xi^-$ & 1321.7 & -e & $1.6 \times 10^{-10}$ & dss & 1/2 \\
\end{tabular}
\end{center}
Although the $\Lambda$ particle has a short lifetime, it can exist
long enough to be "trapped" within a nucleus. Nuclei that have a
$\Lambda$ particle trapped inside are called Lambda-hypernuclei.
The $\Lambda$ is neutral,
so it does not experience the electrostatic force. However, it is
attacted to neutrons and protons via the strong force. The $\Lambda$
"lives" long enough for experimentalists to measure the quantized energy
levels of the $\Lambda$ when it is in the nucleus. If one wanted to
put a $\Lambda$ particle in say a $^{12}C$ nucleus, one would collide
$K^-$ particles on a $^{13}C$ nucleus. One of the neutrons in the
$^{13}C$ carbon nucleus would be converted into a $\Lambda$ and a pion
would be emitted:
\begin{equation}
K^- + n \rightarrow \Lambda + \pi^-
\end{equation}
\noindent By measuring the energy of the $\pi^-$ one can determine the
binding energy of the $\Lambda$ within the $^{12}C$ nucleus.
\newpage
Theoreticians can calculate the binding energies of the $\Lambda$ by
describing the interaction of the $\Lambda$ with the rest of the
nucleus with a potential $V(\vec{r})$. The potential function $V(\vec{r})$
is the potential energy due to the strong interaction between the
$\Lambda$ and the rest of the nucleus, where $\vec{r}$ is the position
vector of the $\Lambda$ from the center of the nucleus. There are some properties that
we believe $V(\vec{r})$ should have:
\begin{enumerate}
\item $V(\vec{r})$ should be spherically symmetric. That is,
$V(\vec{r}) = V(r)$ where $r=|\vec{r}|$.
\item Since the range of the interaction between the $\Lambda$ and
a nucleon is very short, $V(r)$ should go to zero just outside the
nucleus.
\item Since the $\Lambda -nucleon$ interaction is very short range,
the $\Lambda$ is only affected by its nearest neighbor nucleons.
Thus, $V(r)$ should be fairly constant for $rR
\end{eqnarray*}
\noindent where $R$ is the radius of the nucleus, and $V_0$ is the "depth"
of the potential.
To solve for the energy levels that the $\Lambda$ can have, one uses
$V(r)$ as the potential in the Schroedinger equation. Solving the
Schroedinger equation for the spherical square well potential (we will show this
later), one obtains for $l=0$:
\begin{equation}
tan \Big( \sqrt{ {{2m(V_0-|E|)} \over {\hbar^2}} } \; R \Big) = -\sqrt{{{V_0-|E|} \over {|E|}}}
\end{equation}
\noindent where $m$ is the mass of the $\Lambda$, and the binding energy is $|E|$.
To solve this equation for $|E|$, we need to choose the appropriate units and cast
the equation into its simplest form.
\newpage
In particle physics a common unit for mass is $MeV/c^2$, or $GeV/c^2$. One ususally
leaves off the $c^2$. A common unit for length is the Fermi, $fm$, where
$1 \; fm = 10^{-15}m$. The Fermi is roughly the radius of a neutron or proton.
Another important constant that enters in quantum mechanics is Plank's constant, $h$,
or $\hbar = h/(2 \pi)$. Usually $\hbar$ can be coupled with $c$, as $\hbar c$.
$\hbar c$ has units of (energy)(distance). In particle physics units,
$\hbar c \approx 197.33 \; MeV-fm$.
To solve for the potential strength $V_0$, it is easiest to express the equation in
unitless quantities. This can be done as follows:
\begin{eqnarray*}
tan \Big( \sqrt{ {{2m(V_0-|E|)} \over {\hbar^2}} } \; R \Big) & = &-\sqrt{{{V_0-|E|} \over {|E|}}}\\
tan \Big( \sqrt{ {{2mc^2(V_0/|E|-1)} \over {\hbar^2c^2}} } \; R \Big) & = &-\sqrt{V_0/|E|-1}
\end{eqnarray*}
\noindent If we let $\alpha = 2mc^2|E|R^2/(\hbar c)^2$ and $x=V_0/|E|$, the equation becomes
\begin{equation}
tan(\sqrt{\alpha (x-1)} \; ) = - \sqrt{x-1}
\end{equation}
\noindent or
\begin{equation}
tan(\sqrt{\alpha (x-1)} \; ) + \sqrt{x-1} = 0
\end{equation}
\newpage
\begin{center}
{\bf Solving the Non-Relativistic Schroedinger Equation \\
for a spherically symmetric potential}
\end{center}
\bigskip
If the energy of a particle is non-relativistic, and its interaction
is described by a potential energy function, the "physics" is described
by solutions to the the time independent Schr\"{o}dinger equation:
\begin{equation}
-{{\hbar^2}\over{2m}} \nabla^2 \Psi + V(r) \Psi = E \Psi
\end{equation}
\noindent Whether one is performing a scattering experiment or measuring the
bound state energies, will determine the boundary conditions of the
solution $\Psi (\vec{r})$.
\noindent For bound state solutions, the wavefunction $\Psi$ and the integral
$\int \Psi^* \Psi dV$ over all space must be finite.
Thus the "boundary conditions" at infinity are:
$r \rightarrow \infty$, $\Psi$ must approach zero
faster than $1/r$. Since the potential is spherically symmetric,
the angular dependence can be separated from the radial. Writing
$\Psi = R(r) Y_{lm}(\theta, \phi)$ as a product of a radial part times a
{\it spherical harmonic} ($Y_{lm}(\theta, \Phi)$), the above equation reduces to
\begin{equation}
-{{\hbar^2}\over{2m}} ({1 \over r^2} {{d} \over {d r}}
(r^2 {{d R} \over {d r}}) - {{l(l+1)} \over {r^2}} R(r))
+ V(r) R(r) = E R(r)
\end{equation}
\noindent The integer $l$ is related to the particles {\it orbital
angular momentum}.
A further simplification is obtained by writing $R(r)$ as
$u(r)=R(r)/r$. The radial part of the Schr\"{o}dinger equation finally
becomes the somewhat simple form:
\begin{equation}
-{{\hbar^2}\over{2m}} ( {{d^2 u(r)} \over {d r^2}}
- {{l(l+1)} \over {r^2}} u(r) )
+ V(r) u(r) = E u(r)
\end{equation}
\noindent For $\Psi$ to be finite, $u(0)$ equals $0$, and for bound
states, $u(r)$ goes to zero as $r \rightarrow \infty$.
\newpage
For $l=0$, no orbital angular momentum, the equation further simplifies to
\begin{equation}
-{{\hbar^2}\over{2m}} {{d^2 u(r)} \over {d r^2}}+ V(r) u(r) = E u(r)
\end{equation}
For the spherical square well potential, we can solve the Schr\"{o}dinger equation
exactly for $rR$. Then, we can require that $u(r)$ and it's
derivative are continuous at $r=R$.
For $rR$, $V(r)=0$, and we have
\begin{eqnarray*}
{{d^2 u(r)} \over {d r^2}} & = & {{2m|E|} \over \hbar^2} u(r) \\
{{d^2 u(r)} \over {d r^2}} & = & k^2 u(r)
\end{eqnarray*}
\noindent where $k=\sqrt{2m|E|/\hbar^2}$. The solution to this equation that
has $u(r \rightarrow \infty) \rightarrow 0$ is
\begin{equation}
u(r) = Be^{-kr}
\end{equation}
\newpage
\noindent Requiring $u(r)$ to be continuous at $R$ gives:
\begin{equation}
A sin(k'R) = Be^{-kR}
\end{equation}
\noindent and requiring $u'(r)$ to be continuous at $R$ gives:
\begin{equation}
Ak'cos(k'R) = -Bke^{-k'R}
\end{equation}
\noindent Dividing these two equations yields:
\begin{equation}
{{tan(k'R)} \over k'} = - {1 \over k}
\end{equation}
\noindent Expressing this equation with our original variables gives
\begin{eqnarray*}
tan(k'R) & = & - {k' \over k} \\
tan \Big( \sqrt{ {{2m(V_0-|E|)} \over {\hbar^2}} } \; R \Big) & = & -\sqrt{{{V_0-|E|} \over {|E|}}}
\end{eqnarray*}
\noindent which is the equation that we are solving numerically.
\newpage
\centerline{\bf Appendix}
\centerline{\bf Bisection Method}
The bisection method finds roots in the interval between the values
$a$ and $b$, where $a** tol)$ \\
$\{$ \\
c=(a+b)/2; \\
if $(f(c)*f(a)<0)$ then b = c; \\
if $(f(c)*f(b)<0)$ then a = c; \\
$\}$ \\
\noindent where "tol" can be determined from the accuracy desired or the machine
tolerance. Instead of $f(c)$ being the condition, the while loop could continue
till $|b-a|$ is less than a tolerance.
One limitation of the bisection method is that one needs to choose $a$ and $b$
such that there is a root between them. The next method avoids this requirement,
and only needs a starting value and step size.\\
\centerline{\bf Newton's Method}
Newton's method can be used if the first derivative, $f'(x)$, is known.
We will just state the algorithm here. For a derivation as to its
validity, refer to a text on numerical methods.
One starts with an initial guess for the root of $f(x)$, $x_1$. The
algorthym is an equation to obtain $x_{i+1}$ from $x_i$:
\begin{equation}
x_{i+1} = x_i - {{f(x_i)} \over {f'(x_i)}}
\end{equation}
\noindent Under certain conditions, the series $x_i$ will converge to
a root of $f(x)$. Usually in physics applications one does not know
the analytic form of $f'(x)$, so Newton's method is not commonly
used.
\end{document}
**