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?

  • $\begingroup$ Is using Mathematica an option? $\endgroup$
    – user21
    Jun 26, 2017 at 16:16
  • $\begingroup$ Due to type of potential, I don't think Mathematica is an appropriate method. The potential is kind of central potential like Yukawa and Lennard Jones. $\endgroup$
    – kelasmadin
    Jun 27, 2017 at 8:38
  • $\begingroup$ Can you add what k(t) and V(r) look like? $\endgroup$
    – user21
    Jun 27, 2017 at 20:11
  • $\begingroup$ For example, Yukawa potiential $V(r) = -V_0 \frac{e^{-\alpha r}}{r}$ and $k(t) = k_0 - F_0 \cos(\omega t)$ $\endgroup$
    – kelasmadin
    Jun 28, 2017 at 6:12
  • $\begingroup$ Could you provide more information regarding the number of particles you wish to use, domain, dimensions, etc? $\endgroup$
    – knl
    Jun 28, 2017 at 11:14

2 Answers 2


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).

  • 1
    $\begingroup$ I didn't completely catch your mean.As your mentioned, I inserted the expansion into the equation.\begin{align*}\sum_{\ell}\Big[ \frac{1}{R_{\ell}} \frac{d}{dr}(r^2 \frac{dR_{\ell}}{dr})-\frac{2mr^2}{\hbar^2}\big(\frac{1}{2}k(t)(x^2 + y^2)+V(r) \big)+\nonumber\\ \frac{2mr^2}{\hbar}\frac{1}{R_{\ell}} \frac{\partial R_{\ell}}{\partial t}\Big]+\sum_{\ell} \frac{1}{Y_{\ell 0}} \Big( \frac{1}{\sin\theta} \frac{\partial}{\partial \theta}\big(\sin \theta \frac{\partial Y_{\ell 0}}{\partial \theta}\big) \Big)=0.\end{align*}we have two terms separated.Should I solve numerically the first term? $\endgroup$
    – kelasmadin
    Jul 1, 2017 at 18:10
  • $\begingroup$ You project out the radial part, and by this you should notice that the spherical harmonics are eigenfunctions of the radial part of the Laplace operator. See here, for example. But my bad, instead of $R_l(r)$ use $\frac{R_l(r)}{r}$ -- it's more appropriate for the radial part of the Laplace operator. $\endgroup$
    – davidhigh
    Jul 1, 2017 at 20:48
  • $\begingroup$ thanks to you. I was wondering what were your mean of $x^2 + y^2 \sim r^2Y_{20}$? Can I insert it in the equation for numerically solving?($u=rR$) \begin{align*} \frac{d^2u(r,t)}{dr^2} - \frac{2mr^2}{\hbar^2}\Big(\frac{1}{2}k(t)r^2Y_{20} +\frac{\hbar^2 \ell (\ell +1)}{2mr^2}+V(r) \Big)u(r,t)+\frac{2imr^2}{\hbar}\frac{\partial u(r,t)}{\partial t} = 0. \end{align*} $\endgroup$
    – kelasmadin
    Jul 2, 2017 at 14:40
  • $\begingroup$ The idea is: express $x^2+y^2$ in terms of $r^2$ and some appropriate combination of spherical harmonics. For that, go here, and figure out which ones to take (my attempt above is wrong)... once you have it, insert it into the Schrödinger equation ... and only thereafter, you do the projection onto spherical harmonics. $\endgroup$
    – davidhigh
    Jul 2, 2017 at 15:46
  • $\begingroup$ All other terms remain as in the previous link, and the only difficulty is in evaluating the $x^2+y^2$ term. For that, you get terms like $\langle Y_{lm} | Y_{20} | Y_{\bar l \bar m} \rangle$ and similar (where the inner SH comes from your expansion of $x^2+y^2$. You need to evaluate these integrals (again, there closely related to the Wigner3j symbol). In your final equation, there should be no spherical harmonics anymore, it should depend on $\ell$ and $r$ only. $\endgroup$
    – davidhigh
    Jul 2, 2017 at 15:48

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

  • $\begingroup$ I will try to solve it with your suggestion and let you know if it works. $\endgroup$
    – kelasmadin
    Jun 29, 2017 at 12:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.