I'm trying to find the energy eigenvalues of a Dirac delta potential: $$V(x)=-\alpha\delta(x)$$ with periodic boundary condition over some length $L$: $$\psi(x+L)=\psi(x)$$ and only even solutions: $$\psi(x)=\psi(-x).$$
I'm trying to solve it with a matrix method as follows:
- Consider the periodic boundary condition and even solution constraint for a zero potential. The solution basis is $\psi_n(x)=\sqrt{\frac{2}{L}}\cos(\frac{2\pi n}{L}x)$ with energy eigenvalues $E_n=\frac{2n^2\pi^2}{mL^2} ($since $k=\sqrt{2mE}$.
Write the Hamiltonian as $H=H_0+V$, where $H_0$ is the Hamiltonian for the zero-potential: $$H_0=\begin{bmatrix}E_1 & 0 & \dots \\ 0 & E_2 & \dots \\ \vdots & \vdots & \ddots\end{bmatrix}$$ and V is the potential operator expanded into this basis: $$V_{mn}=<\psi_m|V|\psi_n>=-\frac{2\alpha}{L}\int_{-\infty}^\infty\delta(x)\cos(\frac{2m\pi}{L}x)\cos(\frac{2n\pi}{L}x)dx$$
Numerically compute matrix entries for some arbitrarily chosen values of constants, up to a finite matrix size. Numerically compute the eigenvalues.
For a large enough finite basis, the lower-energy eigenvalues should approximately equal the true values. I coded up a solution in Python, and the eigenvalues clearly are converging as the matrix size increases. However, I think I'm doing something wrong and I'm unsure about whether a particular solution should be included in the basis.
Note that: $$\psi_0(x)=\sqrt{\frac{2}{L}}\cos(\frac{2\pi\cdot0}{L}x)=\sqrt{\frac{2}{L}}$$ is a valid solution that is constant everywhere, but has zero energy. I get drastically different solutions depending on whether or not I include it (if I do, the smallest eigenvalue is negative!). Is this physically meaningful, and is there any reason to discard it?
Additionally - I think I've done something wrong along the way, as the values I'm getting differ from a different solution method in both cases. Is there anything wrong with my approach?
EDIT 1:
So after some thought, I've made some changes. I'm now including the zero-energy constant solution in my basis, such that the basis is: $$ \psi_0(x) = \sqrt{\frac{2}{L}} $$ and for $n\in\mathbb{N}$: $$ \psi_n(x)=\sqrt{\frac{2}{L}}\cos(\frac{2n\pi}{L}x)$$ Note that all basis functions are equal to $\sqrt{2/L}$ at $x=0$, so: $$V_{mn}=<\psi_m|V|\psi_n>=\int_{-\infty}^\infty-\alpha\delta(x)\psi_m^*(x)\psi(x)dx=-\frac{2\alpha}{L}$$ so that all entries of $V$ are the same: $$V=\begin{bmatrix}-2\alpha/L & -2\alpha/L & \dots \\ -2\alpha/L & -2\alpha/L & \dots \\ \vdots & \vdots & \ddots\end{bmatrix}$$ and the bare Hamiltonian is just: $$H_0=\begin{bmatrix}E_0 & 0 & \dots \\ 0 & E_1 & \dots \\ \vdots & \vdots & \ddots\end{bmatrix}$$ where $E_n=\frac{2n^2\pi^2}{mL^2}$. I've implemented it and it does actually seem to converge in a reasonable manner, with the lowest eigenvalue negative (bound) and the subsequent ones positive. However, the bound state energy does not match the expected value and the higher eigenvalues also seem to be wrong.
What might I be doing wrong? Is this matrix formulation correct? Are my bare Hamiltonian values (allowed energies due to periodic boundary condition) wrong?
EDIT 2:
I've now noticed that the bound state energy does actually approach the continuum limit for larger box sizes. Perhaps it's the other method I'm comparing to that's wrong...
EDIT 3:
I calculated the bound state eigenvector numerically for a large matrix size, then used its entries as coefficients on the basis functions, and the resulting function is exactly what you'd expect for the perturbed bound state (note that it's smooth at the edges, but they're nonzero as you'd expect):
I think my implementation must be right, so I'm going to go ahead and consider this solved!
