# Schrödinger equation with time dependent Hamiltonian

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?

• Is using Mathematica an option? – user21 Jun 26 '17 at 16:16
• 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. – kelasmadin Jun 27 '17 at 8:38
• Can you add what k(t) and V(r) look like? – user21 Jun 27 '17 at 20:11
• For example, Yukawa potiential $V(r) = -V_0 \frac{e^{-\alpha r}}{r}$ and $k(t) = k_0 - F_0 \cos(\omega t)$ – kelasmadin Jun 28 '17 at 6:12
• Could you provide more information regarding the number of particles you wish to use, domain, dimensions, etc? – knl Jun 28 '17 at 11:14

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

• 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? – kelasmadin Jul 1 '17 at 18:10
• 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. – davidhigh Jul 1 '17 at 20:48
• 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*} – kelasmadin Jul 2 '17 at 14:40
• 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. – davidhigh Jul 2 '17 at 15:46
• 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. – davidhigh Jul 2 '17 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.

• I will try to solve it with your suggestion and let you know if it works. – kelasmadin Jun 29 '17 at 12:36