Schrödinger equation with time dependent Hamiltonian

Question

I need to solve the Schrödinger equation with a time dependent Hamiltonian

$$i\hbar \frac{\partial}{\partial t} \Psi = \left[-\frac{\hbar^2}{2m}\nabla^2 +\frac{1}{2} k(t)(x^2+y^2) + V(r)\right]\Psi $$

Can anybody recommend me an efficient numerical method or software package for solving the problem.

EDIT: This is a suitable form of the equation for numerically solving(which we discussed in the one of the below answers with @davidhigh). \begin{align*} i\frac{\partial}{\partial t}u_{\ell}(r,t) = \Bigg(-\frac{1}{2} \frac{\partial ^2}{\partial r^2} + \frac{\ell(\ell+1)}{2r^2} + V(r)\Bigg)u_{\ell}(r,t)\quad \quad \quad \quad \quad \quad \quad \nonumber\\ + \frac{2}{3}k(t)r^2 \sum_{\ell' = {max}(\ell-2,0)}^{{min}(\ell+2,L_{{max}})} \Bigg(\delta_{\ell,\ell'} - \sqrt{ \frac{4\pi}{5}} \alpha(\ell, \ell') \Bigg) u_{\ell'}(r,t) \end{align*} where the coefficient $\alpha(\ell, \ell')$ is given by \begin{align*} \alpha(\ell, \ell') = \int Y^*_{\ell0}(\Omega)Y_{20}(\Omega)Y_{\ell'0}(\Omega) d\Omega \end{align*}

which the integral can be solved by Wigner-3j coefficients (see here ).

Can anyone suggest an efficient numerically method for solving the equation?

Answer 1

I'd suggest to try it on your own. Do an expansion of your wavefunction in terms of spherical harmonics,

$$ \psi(\mathbf r) \ = \ \sum_{\ell} R_\ell(r,t) \, Y_{\ell 0} (\theta,\phi)\,. $$

Note that I've set the index $m$ in $Y_{\ell m}$ to zero, in order to account for the symmetry of your Hamiltonian with respect to rotations around the $x-y$ plane. This makes your problem essentially two-dimensional. You can also write the spherical harmonic in terms of Legendre polynomials, $$ Y_{\ell0}(\theta,\varphi) = \sqrt{\frac{2\ell+1}{4\pi}} P_{\ell}(\cos\theta). $$

Insertion into your Schrödinger equation and project on the spherical hamonics, and you'll end up with a coupled equation for the radial functions $R_l(r,t)$. Here, "coupled" means coupled in $\ell$, i.e. the function $R_\ell(r,t)$ depends on all the other quantum numbers $0\leq \ell \leq L_{\max}$. Solve that with standard finite differences, it's not that hard.

In the insertion-and-projection step above, the only problem appears in evaluating the matrix elements with $x^2+y^2$. If you write it as $r^2 Y_{10}$, it boils down to an integral over three spherical harmonics, which is related to the Clebsch-Gordan or Wigner-3j coefficients. But It's an easy one, for which analytical formulae exist (--just google for the buzzwords in the previous sentence).

If you arrived at the working formula and need further assistance, let me know.

---

EDIT: summarizing our lengthy discussion in the comment section, here is the final equation which is about to be solved numerically.

$$ i\frac{\partial}{\partial t} u_\ell(r,t) = \left(-\frac{1}{2} \frac{\partial^2}{\partial r^2} + \frac{\ell(\ell+1)}{2r^2} + V(r)\right) u_\ell(r,t) \\ \qquad \qquad\qquad\qquad\qquad\quad+ \frac{2}3 \,k(t)\, r^2 \; \sum_{\ell^\prime=\max(\ell-2,0)}^{\min(\ell+2,L_\max)} \left( \delta_{\ell,\ell^\prime} - \sqrt{\frac{4\pi}5} \alpha(\ell,\ell^\prime) \right)\,u_{\ell^\prime}(r,t) $$

Here the coefficient $\alpha(\ell,\ell^\prime)$ which you introduced is given by $$ \alpha(\ell,\ell^\prime) \ = \ \int Y^\ast_{\ell 0}(\Omega) Y_{20}(\Omega) Y_{\ell^\prime 0}(\Omega)\ d\Omega $$ (you can also express that more in standard terms such as Wigner3j symbols, see e.g. here).

Note the restriction in the summation indices which come from the fact that $0\leq \ell^\prime \leq L_\max$ (where $L_\max$ is the maximum angular quantum number chosen in the numerical representation).

Answer 2

FEniCS users have solved this problem before, but keep in mind that FEniCS does not natively support complex numbers right now in its code. Therefore you have to make a workaround.

See: https://fenicsproject.org/qa/9209/how-to-use-complex-numbers-iterative-solvers https://fenicsproject.org/qa/10671/mixed-function-spaces-for-complex-valued-problems-in-2016-0