How can this IBVP with regularity boundary conditions be solved?

Question

I have a radial Schrödinger equation for a particle in Coulomb potential:

$$i\partial_t f(r,t)=-\frac1{r^2}\partial_r\left(r^2\partial_r f(r,t)\right)-\frac2rf(r,t)$$

with initial condition

$$f(r,0)=e^{-r^2}$$

and boundary conditions

$$\begin{cases} |f(0,t)|<\infty\\ |f(\infty,t)|<\infty. \end{cases}$$

Trying to solve it, I couldn't come up to any analytical solution, so I tried to solve it numerically.

Simplest what I thought of was finite differences method. But while I can limit domain to make it finite, imposing zero Dirichlet condition at $r=A$ for some $A$, I would need to somehow impose the condition of regularity at $r=0$, where $f$ doesn't have to be neither zero, nor non-zero in general. I don't know how to do this with finite differences.

Another approach I tried was to expand the initial condition in eigenfunctions of the Hamiltonian operator (the RHS of the PDE) and then approximate the solution taking finite number of eigenstates into account. But the expansion appeared to converge extremely slowly, and after taking 200 bound eigenstates (simple because of closed-form solution in terms of Laguerre polynomials) and 70 eigenstates of continuous spectrum (taking those of them which vanish at $r=A$, so as to impose zero boundary condition at $A$; had $A=10$), I still got nothing similar to initial function.

So, the question is: is there any explicit solution for this IBVP (be it closed-form one or in some sort of series, but necessarily explicit and fast convergent)? If no, how can it be solved numerically?

Answer 1

Regarding the regularity condition when applying finite differences, some authors make the following approach:

Consider the RHS of the equation at $r \to 0 $ and compute the limit with the help of L'Hôpital's rule:

$$\lim_{r\to 0} \left(- \frac{1}{r^2} \partial_r (r^2 f_r) - \frac{2}{r} f \right) = \lim_{r\to 0} \left( - f_{rr} - \frac{2}{r} f_r - \frac{2}{r} f \right) = \lim_{r\to 0} \left( -f_{rr} - 2 f_{rr} - 2 f_r \right), $$

and hence you have:

$$ - \frac{1}{r^2} \partial_r (r^2 f_r) - \frac{2}{r} f \approx -3 f_{rr} - 2 f_r.$$

Now, you may discretize your equation at $i = 0$, if $\, r_i = r_0 + i \Delta r$, using this formula, avoiding the geometrical singularity.

Although this is a pretty "rough" way to approximate both the Laplacian and the source term, it's quite well known when applying a finite difference scheme.

I hope this helps you.

Cheers.

Answer 2

Another approach to set a regularity condition is to make a change of variables to switch the type of boundary condition. For this equation, if $f(r)$ is regular at $r=0$, then $\frac{\text{d}}{\text{d}\rho}f(\rho^2)=0$ at $\rho=0$.

We can switch to $g(\rho)=f(\rho^2)$ and rewrite the equation for $g(\rho)$:

$$i\partial_tg(\rho,t)=-\frac1{4\rho^2}\partial_{\rho\rho}g(\rho,t)-\frac3{4\rho^3}\partial_\rho g(\rho,t)-\frac2{\rho^2}g(\rho,t)$$

with boundary conditions:

\begin{cases} g^{(1,0)}(0,t)=0\\ |g(\infty,t)|<\infty, \end{cases}

i.e. we switched from regularity condition at $r=0$ to homogeneous Neumann boundary condition at $\rho=0$. Now this is trivial to setup for finite differences.

---

Yet simpler approach is to use $q(r,t)=rf(r,t)$. Then the equation transforms to:

$$i\partial_t q(r,t)=-\partial_{rr}q(r)-\frac2rq(r,t),$$

and the regularity condition becomes homogeneous Dirichlet condition (because finite $f$ is multiplied by zero $r$). Now we can't set good boundary condition at infinity though, but for numerical treatment this is not a problem. For finite domain $r\in[0,L]$ we have the following boundary conditions:

\begin{cases} q(0,t)=0\\ q(L,t)=0, \end{cases}

and this problem is also trivial to setup for finite differences.

Answer 3

This may be a wild goose chase, and at the very least it will be a mess, but here are my thoughts. The time-independent Schrodinger equation (TISE) with Coulomb potential can be solved exactly -- the negative-energy solutions are in any quantum mechanics book and the positive-energy solutions are confluent hypergeometric functions of some sort (I believe they are given in Schiff, for instance). So you can write a general solution of the time-dependent equation (TDSE) as a linear combination of all solutions of the TISE each multiplied by exp(-iEt); the coefficients are determined by the initial condition. (You only need the l=0 solutions, of course, since your initial condition is rotationally invariant.) I have not mentioned boundary conditions, and indeed I can only say that I think the above solutions of the TISE are well-behaved at the origin. The positive-energy solutions surely don't go to zero as r->infinity but I think the sum must do so.