Q31 - Circular well with linear potential

Consider a particle of mass $m$ in 2D, which is moving in a linear potential in a circular well of radius $R$.

$$ -\frac{\hbar^2}{2m} \nabla^2\Psi(x, y) + V (x, y)\Psi(x, y) = E\Psi(x, y) $$

where the potential is

$$ V (x, y) = \begin{cases} 0, & x^2 + y^2 > R^2 \\ -V_0 \left( 1 - \beta \frac{\sqrt{x^2+y^2}}{R} \right), & x^2 + y^2 \leq R^2 \end{cases} $$

This problem has an exact solution for $\beta = 0$, which you can look up and compare with your numerical results.

Consider values of potential $V_0 \frac{\hbar^2}{mR^2} = 0.1, 1.0, 10$ and $\beta = 0, 0.5, 1$ make tables for all possible combinations of these parameters:

• Calculate the ground state energy in units of $\frac{\hbar^2}{mR^2}$.
• Plot the ground state probability density $|\Psi(x, y)|^2$.
• Calculate the number of bound states.

For $\beta = 1$,

• What is the minimum $V_0$ that gives a bound state with $\langle L_z \rangle = \langle xp_y - yp_x \rangle = 2\hbar$?
• Plot $\frac{\langle \sqrt{x^2 + y^2}\rangle}{R}$ for the ground state as a function of $V_0 \frac{\hbar^2}{mR^2} = 0.1 \to 5$.

Solution

First of all, I will try to compare the methods I have and I will choose the optimum one for my purposes.

---

Finite Difference Method in Cartesian Coordinates (2D)

The first thing that comes to my mind is brute forcing the cartesian coordinates. If I want to use regular finite difference method, I would need a 3D Hamiltonian matrix. Instead, I will use kronecker sums and products to reduce the Hamiltonian in 2D. However, the wave function vector $\psi$ will now have $N^2$ elements instead of $N$. Let's visualize this procedure for $3 \times 3$ case.

A 2D wave function is a grid. For a simple $3 \times 3$ example, the grid looks like this:

$$ \Psi_{grid} = \begin{bmatrix} \psi_{11} & \psi_{12} & \psi_{13} \\ \psi_{21} & \psi_{22} & \psi_{23} \\ \psi_{31} & \psi_{32} & \psi_{33} \end{bmatrix} $$

This would require a 3D hamiltonian which I do not know how to construct it. Nevertheless, we can use a trick by converting this grid to a vector:

$$ \vec{\psi} = \begin{bmatrix} \psi_{11} \\ \psi_{21} \\ \psi_{31} \\ \mathbf{\psi_{12}} \\ \mathbf{\psi_{22}} \\ \mathbf{\psi_{32}} \\ \psi_{13} \\ \psi_{23} \\ \psi_{33} \end{bmatrix} $$

Now if I can define a Hamiltonian acting on this properly, I will find the eigenvalues so the bound states. Let's try to express the 2D laplacian in this form:

We define the basic 1D difference operator $D$ (ignoring the $1/(\Delta x)^2$ factor for now):

$$ D = \begin{bmatrix} -2 & 1 & 0 \\ 1 & -2 & 1 \\ 0 & 1 & -2 \end{bmatrix} $$

This assumes, at the boundaries $\psi(x)=0$. Now let's try to construct the $\frac{\partial^2}{\partial x^2}$ for this new $\vec{\psi}$ using the matrix $D$. In our vector $\vec{\psi}$, elements that are neighbors in the $x$-direction (e.g., $\psi_{11}$ and $\psi_{21}$) are adjacent to each other in the vector list.

This means the differentiation happens locally within the small blocks of the vector. Mathematically, this corresponds to the Kronecker Product $I \otimes D$:

$$ \frac{\partial^2}{\partial x^2} \rightarrow I \otimes D = \begin{bmatrix} \mathbf{D} & 0 & 0 \\ 0 & \mathbf{D} & 0 \\ 0 & 0 & \mathbf{D} \end{bmatrix} = \left[ \begin{array}{ccc|ccc|ccc} -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & -2 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 1 & -2 & 0 & 0 & 0 & 0 & 0 & 0 \\ \hline 0 & 0 & 0 & -2 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & -2 & 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & -2 & 0 & 0 & 0 \\ \hline 0 & 0 & 0 & 0 & 0 & 0 & -2 & 1 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 & 1 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & -2 \end{array} \right] $$

Similarly for $\frac{\partial^2}{\partial y^2}$, we can apply this procedure. Neighbors in the $y$-direction (e.g., $\psi_{11}$ and $\psi_{12}$) are far apart in our vector. $\psi_{11}$ is at index 0, but $\psi_{12}$ is at index 3.

This means we need a matrix that connects different blocks. This corresponds to $D \otimes I$:

$$ \frac{\partial^2}{\partial y^2} \rightarrow D \otimes I = \begin{bmatrix} -2I & I & 0 \\ I & -2I & I \\ 0 & I & -2I \end{bmatrix} = = \left[ \begin{array}{ccc|ccc|ccc} -2 & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 \\ 0 & -2 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & -2 & 0 & 0 & 1 & 0 & 0 & 0 \\ \hline 1 & 0 & 0 & -2 & 0 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 0 & -2 & 0 & 0 & 1 & 0 \\ 0 & 0 & 1 & 0 & 0 & -2 & 0 & 0 & 1 \\ \hline 0 & 0 & 0 & 1 & 0 & 0 & -2 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 0 & 0 & -2 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 & -2 \end{array} \right] $$

Thankfully, we have a kronecker sum property $A \oplus B = A \otimes I + I \otimes B$. Then,

$$ \nabla^2 \approx (D \oplus D) $$

Finally, we construct the full Hamiltonian $H$. Note that the potential $V$ must also be flattened into a diagonal matrix so it matches the shape of our vector.

$$ H = -\frac{\hbar^2}{2m(\Delta x)^2}(D \oplus D) + \text{diag}(V(x,y)) $$

For more details, you can check the video.

We can now feed this matrix $H$ into a solver to find the eigenvalues (Energy $E$) and eigenvectors (Wave functions $\vec{\psi}$). However, this would be a tough method. Because for $N \times N$ grid, I need $N^2 \times N^2$ Hamiltonian. Let's see how it goes:

# Libraries
import numpy as np
import matplotlib.pyplot as plt
from scipy.sparse.linalg import eigsh
from scipy.linalg import eigh_tridiagonal
from scipy.linalg import eig  
from scipy import sparse

# Defining the circular potential in cartesian coordinates
def potential_cart(x, y, beta, V0, R):
    r = np.sqrt(x**2 + y**2)
    V = np.zeros_like(r)
    indexes = r <= R
    V[indexes] = -V0 * (1 - beta * r[indexes] / R)
    return V

def cartesian_bound_state_finder(V0,beta):
    # Parameters, I did not scale them explicitly so I will assume them 1
    hbar = 1
    R = 1
    m = 1
    # Count of eigenvalues to find
    if V0 < 2:
        find = 3
    else:
        #I just find 7 eigenvalues because I tried many times and I know the 7th eigenvalue will be positive for V = 10.
        find = 7
    V0 = V0*hbar**2/m/R**2
    
    # Space grid
    space = 10          # I know this is a small space for shallow wells since wave function extends out of the well significantly. 
    N = 150             # But it cannot find anyways so I choose a small space for the time concerns.
    X, Y = np.meshgrid(np.linspace(-space, space, N),np.linspace(-space, space, N))
    dx = X[0,1]- X[0,0]
    
    # Constructing Hamiltonian
    diag = np.ones(N)
    diags = np.array([diag,-2*diag,diag])
    
    # D matrix that I defined above
    D = sparse.spdiags(diags, np.array([-1,0,1]),N,N)
    
    V = potential_cart(X, Y, beta, V0, R)
    # Kinetic Energy Part
    T = -1/2 * sparse.kronsum(D,D) / (m*dx**2)
    # Potential Energy Part
    U = sparse.diags(V.reshape(N**2), (0))
    # Hamiltonian
    H = T+U

    
    # Finding the bound state energies and wave vector,
    eigenvalues, eigenvectors = eigsh(H, k=find, which='SA')
    
    # Only the negative eigenvalue correspond to a bound state
    bound_states = eigenvalues[eigenvalues<0]

    # Displaying the result
    if len(bound_states) == 0:
        print(f'No bound states found for V0 = {V0:.2f} and beta = {beta:.2f}')
        print()
    else:
        print(f'Bound state energies for V0 = {V0:.2f} and beta = {beta:.2f} are:')
        for bound_state in bound_states:
            print(f'{bound_state:.4f}')
        print()

# Iteration over given conditions
V0s = [0.1,1.0,10.0]
betas = [0,0.5,1.0]
for V0 in V0s:
    for beta in betas:
        cartesian_bound_state_finder(V0,beta)

No bound states found for V0 = 0.10 and beta = 0.00

No bound states found for V0 = 0.10 and beta = 0.50

No bound states found for V0 = 0.10 and beta = 1.00

Bound state energies for V0 = 1.00 and beta = 0.00 are:
-0.1904

Bound state energies for V0 = 1.00 and beta = 0.50 are:
-0.0662

No bound states found for V0 = 1.00 and beta = 1.00

Bound state energies for V0 = 10.00 and beta = 0.00 are:
-8.0928
-5.2838
-5.2838
-2.0245
-1.6466
-1.0234

Bound state energies for V0 = 10.00 and beta = 0.50 are:
-5.7639
-2.4380
-2.4380

Bound state energies for V0 = 10.00 and beta = 1.00 are:
-3.6566
-0.0604
-0.0604

This method seemingly works. It yields some bound states for $V = 10.0$ and $V = 1.0$. However, for $V=0.1$, I do not get any bound states. This is unphysical since it is proven that there must be at least 1 bound state for any given circular potential well in 2D Reference. In the results above, some energy levels are identical or very close (for example, $-5.2838$). These correspond to degeneracies caused by angular momentum, which I will discuss in more detail later. To verify these findings, I will first solve the Schrödinger equation analytically for a constant potential well to compare it with my Cartesian results. After that, I will use other numerical methods to catch detect shallow bound states.

---

Polar Coordinates

Let's use to polar coordinates. We will have an ordinary differential equation since the potential is rotationally symmetric. We start with the time-independent Schrödinger equation in 2D Cartesian coordinates:

$$-\frac{\hbar^{2}}{2m}\nabla^{2}\Psi + V(r)\Psi = E\Psi$$

In two-dimensional Cartesian coordinates, the Laplacian is $\nabla^2 = \frac{\partial^2}{\partial x^2} + \frac{\partial^2}{\partial y^2}$. Converting to polar coordinates $(r, \theta)$, this becomes:

$$\nabla^2 = \frac{\partial^2}{\partial r^2} + \frac{1}{r}\frac{\partial}{\partial r} + \frac{1}{r^2}\frac{\partial^2}{\partial \theta^2}$$

Substituting this back into the Schrödinger equation:

$$-\frac{\hbar^{2}}{2m} \left[ \frac{\partial^2 \Psi}{\partial r^2} + \frac{1}{r}\frac{\partial \Psi}{\partial r} + \frac{1}{r^2}\frac{\partial^2 \Psi}{\partial \theta^2} \right] + V(r)\Psi = E\Psi$$

We assume the wavefunction $\Psi(r, \theta)$ can be separated into a radial part $R(r)$ and an angular part $\Theta(\theta)$:

$$\Psi(r, \theta) = R(r)\Theta(\theta)$$

Substituting this form into the equation:

$$-\frac{\hbar^{2}}{2m} \left[ \Theta \frac{d^2 R}{d r^2} + \frac{\Theta}{r} \frac{d R}{d r} + \frac{R}{r^2} \frac{d^2 \Theta}{d \theta^2} \right] + V(r)R\Theta = E R\Theta$$

To separate the terms, dividing the entire equation by $R\Theta$ and multiply by $-\frac{2mr^2}{\hbar^2}$:

$$r^2 \left( \frac{1}{R}\frac{d^2 R}{d r^2} + \frac{1}{rR}\frac{d R}{d r} \right) + \frac{1}{\Theta}\frac{d^2 \Theta}{d \theta^2} - \frac{2mr^2}{\hbar^2}(V(r) - E) = 0$$

Since there is no other angular dependance, let's say: $$\frac{d^2 \Theta}{d \theta^2} = -m_l^2 \Theta$$

The solution is $\Theta(\theta) = e^{i m_l \theta}$. For the wavefunction $\Psi(\theta) = \Psi(\theta + 2\pi)$, so $m_l$ must be an integer. Then, we left with the radial part:

$$r^2 \frac{1}{R} \left( \frac{d^2 R}{d r^2} + \frac{1}{r} \frac{d R}{d r} \right) - \frac{2mr^2}{\hbar^2}(V(r) - E) = m_l^2$$

Dividing by $r^2$ gives the final radial differential equation:

$$\frac{d^2 R}{d r^2} + \frac{1}{r}\frac{d R}{d r} + \left[ \frac{2m(E - V(r))}{\hbar^2} - \frac{m_l^2}{r^2} \right] R(r) = 0$$

Here, let me define dimensionless variables $\tilde{r} = \frac{r}{R}$ and $\tilde{E}-\tilde{V} = \frac{(E-V)}{\frac{\hbar^2}{mR^2}}$, and substitute them:

$$\frac{1}{R^2}\frac{d^2 R}{d \tilde{r}^2} + \frac{1}{\tilde{r} R^2}\frac{d R}{d \tilde{r}} + \left[ \frac{2(\tilde{E}-\tilde{V})}{R^2} - \frac{m_l^2}{\tilde{r}^2 R^2} \right] R(\tilde{r}) = 0$$

Simplifying: $$\frac{d^2 R}{d \tilde{r}^2} + \frac{1}{\tilde{r}}\frac{d R}{d \tilde{r}} + \left[ 2(\tilde{E}-\tilde{V}) - \frac{m_l^2}{\tilde{r}^2} \right] R(\tilde{r}) = 0$$

Finally, $$\frac{-1}{2}\frac{d^2 R}{d \tilde{r}^2} + \frac{-1}{\tilde{2r}}\frac{d R}{d \tilde{r}} + \left[ \tilde{V} + \frac{m_l^2}{2\tilde{r}^2} \right] R(\tilde{r}) = ER(\tilde{r})$$

---

Analytical Solution for Constant Potential Well

From now on, I will drop the tildes, let's rearrange the scaled radial Shrödinger's equation: $$-\frac{d^2 R}{d{r}^2} + -\frac{1}{r}\frac{d R}{d {r}} + \left[ 2(V(r)-E) + \frac{m_l^2}{{r}^2} \right] R({r}) = 0$$ $V = -V_0$ inside the well for $\beta = 0$, hence; $${r}^2\frac{d^2 R}{d{r}^2} +{r}\frac{d R}{d {r}} + \left[ 2(E+V_0){r}^2 - m_l^2 \right] R({r}) = 0$$

This differential equation has a very similar form to the Bessel Differential Equation: $$x^2 \frac{d^2y}{dx^2} + x \frac{dy}{dx} + (x^2 - \alpha^2)y = 0$$ which has the general solution: $$y(x) = C_1 J_\alpha(x) + C_2 Y_\alpha(x)$$

Let $k^2 = 2(E+V_0)$. Since we are looking for bound states within a deep well, $E+V_0 > 0$, making $k$ real. Substituting $x = kr$: $$x^2\frac{d^2 R}{d{x}^2} +x\frac{d R}{d{x}} + \left[ x^2 - m_l^2 \right] R({x}) = 0$$ The general solution is $R(r) = A J_{m_l}(kr) + B Y_{m_l}(kr)$. Because $Y_{m_l}(kr)$ diverges at the origin ($r \to 0$), it is unphysical. Thus: $$R_{in}(r) = A J_{m_l}(kr)$$

Outside the well, $V(r) = 0$. For a bound state, the energy $E$ must be negative ($E < 0$). Let $\kappa^2 = -2E$, where $\kappa$ is real. The equation becomes the Modified Bessel Equation: $$r^2 \frac{d^2R}{dr^2} + r \frac{dR}{dr} - (\kappa^2 r^2 + m_l^2)R = 0$$ The general solution is: $$R_{out}(r) = C I_{m_l}(\kappa r) + D K_{m_l}(\kappa r)$$

$I_{m_l}(\kappa r)$ is the Modified Bessel Function of the First Kind. It grows exponentially as $r \to \infty$, so we must set $C = 0$ for a normalizable state. $K_{m_l}(\kappa r)$ is the Modified Bessel Function of the Second Kind. It decays exponentially as $r \to \infty$, which is the physical requirement for a bound state.

Thus: $$R_{out}(r) = D K_{m_l}(\kappa r)$$

From the smoothness of the wavefunction at $r = R$:

1. $A J_{m_l}(kR) = D K_{m_l}(\kappa R)$
2. $A k J_{m_l}'(kR) = D \kappa K_{m_l}'(\kappa R)$

Dividing the second equation by the first (and setting $R=1$) gives the discretization condition: $$\frac{k J_{m_l}'(k)}{J_{m_l}(k)} = \frac{\kappa K_{m_l}'(\kappa)}{K_{m_l}(\kappa)}$$

Rearranging into a root-finding form and using $\kappa = \sqrt{2V_0 - k^2}$: $$F(k) = k J_{m_l}'(k) K_{m_l}(\sqrt{2V_0-k^2}) - \sqrt{2V_0-k^2} K_{m_l}'(\sqrt{2V_0-k^2}) J_{m_l}(k)$$

I will find the roots of $F(k)$ using Brent's algorithm from scipy library.

from scipy.special import jv, jvp, kv, kvp
from scipy.optimize import brentq

def analytic_bound_state_finder(V0, m):
    # Let's first define the function to find its roots
    def function(k):
        kappa = np.sqrt(np.abs(2 * V0 - k**2))
        J, dJ = jv(m, k), jvp(m, k)
        K, dK = kv(m, kappa), kvp(m, kappa)
        return (k * dJ * K) - (kappa * dK * J)

    # evaluating the function for large intervals to use in the Brent's algorithm
    k_max = np.sqrt(2 * V0)  # maximum value of k when E = 0
    k_grid = np.linspace(0,k_max,int(1e4))
    k_grid.sort()
    
    y_vals = [function(k) for k in k_grid]

    
    roots = []
    # Brent's root finding algorithm
    for i in range(len(k_grid) - 1):
        # detecting sign change, if there is, finding the root between these intervals with higher precision using Brent's algorithm
        if y_vals[i] * y_vals[i+1] < 0:
            try:
                k_root = brentq(function, k_grid[i], k_grid[i+1], xtol=1e-12)
                E = 0.5 * k_root**2 - V0
                roots.append(E)
            except ValueError:
                pass
                
    return roots

# Potentials
V0s = [0.1, 1.0, 10.0]

# Finding and displaying the bound states
for V in V0s:
    print(f"\n--- V0 = {V} ---")
    for m in range(4):
        Es = analytic_bound_state_finder(V, m)
        if Es:
            # Scientific notation
            formatted_E = [f'{e:.4e}' for e in Es]
            print(f"m={m}: {len(Es)} bound state(s) -> E = {formatted_E}")
        else:
            print(f"m={m}: No states")
            break

--- V0 = 0.1 ---
m=0: 1 bound state(s) -> E = ['-2.1515e-09']
m=1: No states

--- V0 = 1.0 ---
m=0: 1 bound state(s) -> E = ['-1.9360e-01']
m=1: No states

--- V0 = 10.0 ---
m=0: 2 bound state(s) -> E = ['-8.0942e+00', '-9.4558e-01']
m=1: 1 bound state(s) -> E = ['-5.2645e+00']
m=2: 1 bound state(s) -> E = ['-1.7749e+00']
m=3: No states

It is clear that there is a bound state for $V_0 = 0.1$. Let's try to catch this value numerically using polar coordinates. We have the differential equation: $$\frac{-1}{2}\frac{d^2 R}{d {r}^2} + \frac{-1}{{2r}}\frac{d R}{d {r}} + \left[ {V} + \frac{m_l^2}{2{r}^2} \right] R({r}) = ER({r})$$

Let's try to solve in this form. Discretizations are:

$$ \frac{d^2R}{dr^2} \simeq \frac{R_{n+1}+R_{n-1}-2R_{n}}{(\Delta r)^2}$$ $$ \frac{dR}{dr} \simeq \frac{R_{n+1}-R_{n-1}}{2(\Delta r)}$$

For the boundary conditions, I assume $\psi_{N+1}=0$ and $\psi_{-1}=\psi_{1}$ since the probability density have a maximum at $r=0$. However, the yielding Hamiltonian will not be symmetric because first derivative term is not symmetric. But still I will try to solve in this form and take the real part of eigenvalues. I know this may not be a trustworthy method, but I will give it a try.

# Defining radial potential
def potential(rho, beta, V0, R):
    V = np.zeros_like(rho)
    inside = rho <= R
    V[inside] = -V0 * (1.0 - beta * rho[inside] / R)
    return V
    

def polar_bound_state_finder(R, beta, V0, space, N):
    # Defining the space variables
    r = np.linspace(0, space, N)
    dr = r[1] - r[0]

    
    results = []
    ml = 0

    # Iterating until no bound states left
    while True:

        # Physical potential
        V_r = potential(r, beta, V0, R)

        # Effective potential
        V_eff = V_r.copy()
        # To avoid singularity, I am ignoring the r[0] = 0 term.
        if ml != 0:
            V_eff[1:] += (ml**2) / (2 * r[1:]**2)

        # Building the Hamiltonian
        H = np.zeros((N, N))

        factor_2nd = -1.0 / (2.0 * dr**2)
        factor_1st = -1.0 / (4.0 * r[1:] * dr)

        idx = np.arange(1, N-1)

        # Diagonal
        H[idx, idx] = 1.0 / dr**2 + V_eff[1:-1]

        # Upper diagonal
        H[idx, idx+1] = factor_2nd + factor_1st[:-1]

        # Lower diagonal
        H[idx, idx-1] = factor_2nd - factor_1st[:-1]

        # r = 0 term, in order to avoid divison by zero, I am updating it manually.
        H[0, 0] = 2.0 / dr**2 + V_r[0] 
        H[0, 1] = -2.0 / dr**2   # This comes from boundary condition. I know the wave function is symmetric. 
                                 #Therefore, I assume psi_{-1}=psi_{1} so I multiply this term by 2.

        # Finding the eigenvalues
        vals, vecs = np.linalg.eig(H)
        vals = np.real(vals)

        # Bound states
        bound = vals[vals < 0]

        # Breaking the loop if there is no eigenvalues
        if len(bound) == 0:
            break   

        results.append((ml, np.sort(bound)))
        # next ml state
        ml += 1

    return results


# Parameters
R = 1
space = 15
N = 1200 

# Iteration over given conditions
V0s = [0.1,1.0,10.0]
betas = [0,0.5,1.0]
for V0 in V0s:
    print(f"\n{'='*50}")
    print(f"V0 = {V0}")
    print(f"{'='*50}")

    for beta in betas:
        print(f"\nβ = {beta}")
        print("-"*30)

        results = polar_bound_state_finder(R, beta, V0, space, N)

        if not results:
            print("  No bound states")
            continue

        for m, energies in results:
            energies_str = ", ".join(f"{E:.4e}" for E in energies)
            print(f"  m = {m:2d} | {len(energies)} state(s): {energies_str}")

==================================================
V0 = 0.1
==================================================

β = 0
------------------------------
  No bound states

β = 0.5
------------------------------
  No bound states

β = 1.0
------------------------------
  No bound states

==================================================
V0 = 1.0
==================================================

β = 0
------------------------------
  m =  0 | 1 state(s): -1.9042e-01

β = 0.5
------------------------------
  m =  0 | 1 state(s): -6.6343e-02

β = 1.0
------------------------------
  m =  0 | 1 state(s): -2.5368e-03

==================================================
V0 = 10.0
==================================================

β = 0
------------------------------
  m =  0 | 2 state(s): -8.0779e+00, -8.9108e-01
  m =  1 | 1 state(s): -5.2265e+00
  m =  2 | 1 state(s): -1.7130e+00

β = 0.5
------------------------------
  m =  0 | 1 state(s): -5.7521e+00
  m =  1 | 1 state(s): -2.4003e+00

β = 1.0
------------------------------
  m =  0 | 1 state(s): -3.6377e+00
  m =  1 | 1 state(s): -4.6029e-02

It seems that this is a valid method, causing no troubles. However, it is not better than the cartesian finder. Because this hamiltonian is not symmetric, it takes so much time to evaluate compared to symmetric ones. Nevertheless, it will help me to finding expectation value of $L_z$ and $r$. Let me try to use substitution to get rid of first derivative term. In this way I will have a symmetric hamiltonian which can be evaluated rapidly.

---

Substituted Polar Solution

I now substitute $R(r) = \frac{u(r)}{\sqrt(r)}$ to get rid of first derivative term \begin{align} \frac{dR(r)}{dr} &= \frac{-1}{2r^{3/2}}u(r)+\frac{1}{r^{1/2}}\frac{du(r)}{dr} \\ \frac{d^2R(r)}{dr^2} &= \frac{3}{4r^{5/2}}u(r)-\frac{1}{r^{3/2}}\frac{du(r)}{dr}+\frac{1}{r^{1/2}}\frac{du^2(r)}{dr^2} \\ \frac{3}{4r^{5/2}}u(r)-\frac{1}{r^{3/2}}\frac{du(r)}{dr}&+\frac{1}{r^{1/2}}\frac{du^2(r)}{dr^2} + \frac{-1}{2r^{5/2}}u(r)+\frac{1}{r^{3/2}}\frac{du(r)}{dr} + \left[ 2({E}-{V}) - \frac{m_l^2}{{r}^2} \right] \frac{1}{r^{1/2}}u(r) = 0 \\ \frac{1}{r^{1/2}}\frac{du^2(r)}{dr^2} &+ \frac{1}{4r^{5/2}}u(r) + \left[ 2({E}-{V}) - \frac{m_l^2}{{r}^2} \right] \frac{1}{r^{1/2}}u(r) = 0 \\ \frac{du^2(r)}{dr^2} &+ \left[ 2({E}-{V}) + \frac{1/4-m_l^2}{{r}^2} \right]u(r) = 0 \\ \end{align} Finally, the equation that I need to solve is: $$ \frac{-1}{2}\frac{du^2(r)}{dr^2} + \left[ {V(r)} + \frac{m_l^2-1/4}{{2r}^2} \right]u(r) = Eu(r) $$ Let's try to find bound states using this symmetric method, it is also tridiagonal so the eigenvalues can be found pretty fast.

def substituted_polar_ground_state_finder(R, beta, V0, space, N):
    # parameters
    dr = space / N
    r_axis = np.arange(dr, space + dr, dr)
    ml = 0
    results = []

    # Iterates until no bound states left
    while True:
        # Potential term
        V_eff = potential(r_axis, beta, V0, R) \
                + (ml**2 - 0.25) / (2 * r_axis**2)

        # Creating Hamiltonian diagonals
        main_diag = (1.0 / dr**2) + V_eff
        off_diag  = -1.0 / (2 * dr**2) * np.ones(N - 1)

        # Finding the eigenvalues and eigenvectors of tridiagonal matrix.
        eigenvalues, _ = eigh_tridiagonal(
            main_diag, off_diag, select='i', select_range=(0, 9)
        )
        # Bound states, energies less than 0.
        bound = eigenvalues[eigenvalues < 0]

        if bound.size == 0:
            break

        results.append((ml, bound))
        ml += 1

    return results



R = 1
space = 100
N = 100000

# Iteration over given conditions
V0s = [0.1,1.0,10.0]
betas = [0,0.5,1.0]
for V0 in V0s:
    print(f"\n{'='*50}")
    print(f"V0 = {V0}")
    print(f"{'='*50}")

    for beta in betas:
        print(f"\nβ = {beta}")
        print("-"*30)

        results = substituted_polar_ground_state_finder(R, beta, V0, space, N)

        if not results:
            print("  No bound states")
            continue

        for m, energies in results:
            energies_str = ", ".join(f"{E:.4e}" for E in energies)
            print(f"  m = {m:2d} | {len(energies)} state(s): {energies_str}")

==================================================
V0 = 0.1
==================================================

β = 0
------------------------------
  No bound states

β = 0.5
------------------------------
  No bound states

β = 1.0
------------------------------
  No bound states

==================================================
V0 = 1.0
==================================================

β = 0
------------------------------
  m =  0 | 1 state(s): -1.2161e-01

β = 0.5
------------------------------
  m =  0 | 1 state(s): -2.6438e-02

β = 1.0
------------------------------
  m =  0 | 1 state(s): -2.8301e-05

==================================================
V0 = 10.0
==================================================

β = 0
------------------------------
  m =  0 | 2 state(s): -7.8224e+00, -5.1858e-01
  m =  1 | 1 state(s): -5.2682e+00
  m =  2 | 1 state(s): -1.7807e+00

β = 0.5
------------------------------
  m =  0 | 1 state(s): -5.3677e+00
  m =  1 | 1 state(s): -2.4211e+00

β = 1.0
------------------------------
  m =  0 | 1 state(s): -3.1110e+00
  m =  1 | 1 state(s): -4.5401e-02

This method is much faster, but I could not figure out how to handle the singularity at $r=0$. I tried using the same boundary conditions as the initial polar function, but they did not work this time. Although the energy values I get are close to the bound states, the singularity at the origin prevents the results from reaching the exact analytical values. Increasing the grid size does not solve this; the results still will approach these obtained eigenvalues. It seems that polar coordinate solution without substitution is the most suitable among these. I will continue with that algorithm. However, all of them cannot catch the bound state energies for $V = 0.1$ because the potential difference is pretty small. However, there should be at least 1 bound state for any attractive potential Reference. That's why I will count them as they have 1 bound state even if I cannot explicitly find their values.

---

Part a)

Bound states are:

# Parameters
R = 1
space = 20
N = 2000 

# Iteration over given conditions
V0s = [0.1,1.0,10.0]
betas = [0,0.5,1.0]
for V0 in V0s:
    print(f"\n{'='*50}")
    print(f"V0 = {V0}")
    print(f"{'='*50}")

    for beta in betas:
        print(f"\nβ = {beta}")
        print("-"*30)

        results = polar_bound_state_finder(R, beta, V0, space, N)

        if not results:
            print("  No bound states")
            continue

        for m, energies in results:
            energies_str = ", ".join(f"{E:.4e}" for E in energies)
            print(f"  m = {m:2d} | {len(energies)} state(s): {energies_str}")

==================================================
V0 = 0.1
==================================================

β = 0
------------------------------
  No bound states

β = 0.5
------------------------------
  No bound states

β = 1.0
------------------------------
  No bound states

==================================================
V0 = 1.0
==================================================

β = 0
------------------------------
  m =  0 | 1 state(s): -1.9096e-01

β = 0.5
------------------------------
  m =  0 | 1 state(s): -6.6546e-02

β = 1.0
------------------------------
  m =  0 | 1 state(s): -3.9062e-03

==================================================
V0 = 10.0
==================================================

β = 0
------------------------------
  m =  0 | 2 state(s): -8.0806e+00, -9.0000e-01
  m =  1 | 1 state(s): -5.2326e+00
  m =  2 | 1 state(s): -1.7233e+00

β = 0.5
------------------------------
  m =  0 | 1 state(s): -5.7531e+00
  m =  1 | 1 state(s): -2.4033e+00

β = 1.0
------------------------------
  m =  0 | 1 state(s): -3.6375e+00
  m =  1 | 1 state(s): -4.5805e-02

These are the bound state energies. I will now plot the ground state $\|\psi(x,y)\|^2$s using cartesian solver because it is easier.

def plot_ground_state_grid(V0s, betas):
    # Fixed Parameters
    hbar, R, m = 1, 1, 1
    space, N = 15, 200
    
    x = np.linspace(-space, space, N)
    y = np.linspace(-space, space, N)
    X, Y = np.meshgrid(x, y)
    dx = x[1] - x[0]
    
    # Figure
    fig, axes = plt.subplots(len(V0s), len(betas), figsize=(15, 13))
    fig.suptitle("Ground State Probability Density $|\psi|^2$ in $xy$-Plane", fontsize=20, y=0.95)

    # Kinetic energy part
    diag = np.ones(N)
    diags = np.array([diag, -2*diag, diag])
    D = sparse.spdiags(diags, np.array([-1, 0, 1]), N, N)
    T = -0.5 * sparse.kronsum(D, D) / (m * dx**2)

    # Solving and plotting the probability density
    for i, V0_val in enumerate(V0s):
        # potential
        V0_scaled = V0_val * hbar**2 / (m * R**2)
        
        for j, beta in enumerate(betas):
            ax = axes[i, j]
            
            # Potential and Hamiltonian
            V = potential_cart(X, Y, beta, V0_scaled, R)
            U = sparse.diags(V.reshape(N**2), 0)
            H = T + U
            
            # Solving for the lowest state
            try:
                vals, vecs = eigsh(H, k=1, which='SA')
                energy = vals[0]
                
                if energy < 0:
                    # Ground state density
                    psi = vecs[:, 0].reshape(N, N)
                    prob_density = np.abs(psi)**2
                    
                    # Plotting
                    im = ax.pcolormesh(X, Y, prob_density, cmap='magma', shading='gouraud')
                    ax.set_title(f"$V_0={V0_val}, \\beta={beta}$\n$E={energy:.4f}$", fontsize=10)
                    
                    # Potential well boundaries
                    circle = plt.Circle((0, 0), R, color='cyan', fill=False, lw=1.5, linestyle='--')
                    ax.add_artist(circle)
                else:
                    ax.text(0.5, 0.5, "No Bound State", ha='center', transform=ax.transAxes)
            except:
                ax.text(0.5, 0.5, "Solver Error", ha='center', transform=ax.transAxes)

            # Setting the limits
            ax.set_aspect('equal')
            ax.set_xlim(-R*4, R*4) # Zoom into the well area
            ax.set_ylim(-R*4, R*4)
            if i == 2: ax.set_xlabel("x")
            if j == 0: ax.set_ylabel("y")

    plt.tight_layout(rect=[0, 0.03, 1, 0.95])
    plt.show()

# Calling the function
plot_ground_state_grid([0.1, 1.0, 10.0], [0, 0.5, 1.0])

It seems that my space boundary is good enough to find the bound states because the probability densitites do not exceed the well too much. I have plotted all the bound states that I have found. In addition to these, there are degeneracies for $m_l \neq 0$ like $m_l = \pm 1$ have same energies for both cases. This can be seen in cartesian solution clearly. There occurs some consecutive energies which have pretty close values. These correspond to degeneracies. Therefore, I will multiply each $ml \neq 0$ bound state with 2 to find all bound states. Hence there are,

Potential Depth ($V_0$)	($\beta$)	Total Bound States
$V_0 = 0.1$	$\beta = 0.0$	1
	$\beta = 0.5$	1
	$\beta = 1.0$	1
$V_0 = 1.0$	$\beta = 0.0$	1
	$\beta = 0.5$	1
	$\beta = 1.0$	1
$V_0 = 10.0$	$\beta = 0.0$	6
	$\beta = 0.5$	3
	$\beta = 1.0$	3

---

Part b)

$$ \hat L_z = (\hat x \hat P_y - \hat y \hat P_x) = -i\hbar\frac{\partial}{\partial\theta}. $$

The eigenfunctions can be written in polar coordinates using separation of variables: $$ \Psi(r,\phi)=\psi_{m_l}(r)e^{im_l\theta}, \qquad m_l\in\mathbb{Z}. $$

Acting with the operator $\hat L_z$ yields the eigenvalue equation: $$ \hat L_z \Psi(r,\phi) = -i\hbar\frac{\partial \Psi(r,\theta)}{\partial\phi} = m_l \hbar \Psi(r,\theta) $$

Therefore, in order to find the minimum $V_0$ that gives a bound state with $\langle L_z \rangle = 2 \hbar$, I all need to do is to fix $m_l = 2$ in my code, thanks to the polar coordinate system. Then iterate over potential values.

def ml_bound_state_finder(R, beta, V0, space, N, ml):
    # Defining the space variables
    r = np.linspace(0, space, N)
    dr = r[1] - r[0]

    results = []

    # Physical potential
    V_r = potential(r, beta, V0, R)

    # Effective potential
    V_eff = V_r.copy()
    # To avoid singularity, I am ignoring the r[0] = 0 term.
    if ml != 0:
        V_eff[1:] += (ml**2) / (2 * r[1:]**2)

    # Building the Hamiltonian
    H = np.zeros((N, N))

    factor_2nd = -1.0 / (2.0 * dr**2)
    factor_1st = -1.0 / (4.0 * r[1:] * dr)

    idx = np.arange(1, N-1)

    # Diagonal
    H[idx, idx] = 1.0 / dr**2 + V_eff[1:-1]

    # Upper diagonal
    H[idx, idx+1] = factor_2nd + factor_1st[:-1]

    # Lower diagonal
    H[idx, idx-1] = factor_2nd - factor_1st[:-1]

    # r = 0 boundary
    H[0, 0] = 2.0 / dr**2 + V_r[0]
    H[0, 1] = -2.0 / dr**2

    # Finding the eigenvalues
    vals, vecs = np.linalg.eig(H)
    vals = np.real(vals)

    # Bound states
    bound = vals[vals < 0]
    return bound

R = 1
space = 10
N = 1500

beta = 1
ml = 2

V0_vals = np.linspace(26.0, 30.0, 200) # Starting from 26 for your time

V0_min = None

for V0 in V0_vals:
    bound = ml_bound_state_finder(R, beta, V0, space, N, ml)
    
    if len(bound) > 0:
        V0_min = V0
        print(f"First bound state for m={ml} at V0 = {V0:.3f}")
        break
if V0_min is None:
    print("No bound state found in range")

First bound state for m=2 at V0 = 26.060

This seems valid. Let me also find the $\langle \sqrt{x^2 + y^2} \rangle = \langle r \rangle $ which is a lot easier thanks to polar coordinates.

def ground_state_r_expectation(R, beta, V0, space, N):
    r = np.linspace(0, space, N)
    dr = r[1] - r[0]

    # m = 0 ground state
    V_r = potential(r, beta, V0, R)

    H = np.zeros((N, N))

    factor_2nd = -1.0 / (2.0 * dr**2)
    factor_1st = -1.0 / (4.0 * r[1:] * dr)
    idx = np.arange(1, N-1)

    H[idx, idx] = 1.0 / dr**2 + V_r[1:-1]
    H[idx, idx+1] = factor_2nd + factor_1st[:-1]
    H[idx, idx-1] = factor_2nd - factor_1st[:-1]

    H[0, 0] = 2.0 / dr**2 + V_r[0]
    H[0, 1] = -2.0 / dr**2
    H[-1, -1] = 1.0

    vals, vecs = np.linalg.eig(H)

    # ground state = lowest eigenvalue
    psi = vecs[:, np.argmin(vals)]

    # normalize in 2D
    psi /= np.sqrt(np.trapezoid(np.abs(psi)**2 * r, r))

    # expectation value of r
    r_mean = np.trapezoid(r**2 * np.abs(psi)**2, r)

    return r_mean

R = 1
space = 30
N = 600
beta = 1

V0_vals = np.arange(0.1,5.1,0.1)
r_mean_vals = []

for V0 in V0_vals:
    # use your solver to check ground state exists
    results = polar_bound_state_finder(R, beta, V0, space, N)

    if len(results) == 0 or results[0][0] != 0:
        r_mean_vals.append(np.nan)
        continue

    r_mean = ground_state_r_expectation(R, beta, V0, space, N)
    r_mean_vals.append(r_mean)

# Plotting
plt.figure(figsize=(8, 5))
plt.plot(V0_vals,r_mean_vals, 'o-', label=r'$\langle r \rangle$')
plt.axhline(R, color='r', linestyle='--', label='Well Radius $R$')
plt.title(r'$\langle \sqrt{x^2 + y^2} \rangle$ for the ground state as a function of $V_0 \frac{\hbar^2}{mR^2}$', fontsize=12)
plt.xlabel(r'$V_0$ $\frac{\hbar^2}{mR^2}$')
plt.ylabel(r'Average Radius $\langle r \rangle$')
plt.grid(True, alpha=0.3)
plt.legend()
plt.show()

It seems grounded as the depth increases, the particle gets more confined inside the well.

Conclusion

Finding bound state energies for quantum particles often requires solving eigenvalue problems on a computer, especially when the potential has no analytic solution. While effective, this numerical approach has limitations. For instance, obtaining precision beyond the second significant digit can be difficult because eigenvalue solvers are often slow. The method generally works better for ground state energies than for excited states, and it performs most reliably with deep potential wells. In shallow potentials, finding energy values becomes much harder. Consequently, more negative bound state energies obtained from eigensolvers tend to be more trustworthy than less ones. To improve performance, it is important to use symmetries, as they significantly speed up the code. However, choosing the right coordinate system involves a trade-off. A Cartesian solver avoids mathematical singularities but requires a long time to reach high precision. In contrast, a polar solver is much faster but introduces singularities when symmetry is applied. While I managed to resolve one singularity by implementing the correct boundary conditions, I could not succeed with the substituted case and need to experiment with different boundary settings. A final note on cartesian solver: it helped me to see the degeneracies which at first I thought it was a numerical error. And it made plotting wavefunctions on an $(x,y)$ grid much simpler. Therefore, all methods have their own benefits.