I have been following an excellent article about how to use the Crank-Nicolson method to solve the Schrodinger equation. In the article, it starts with a $V(x, y, t)$ but the potential seems to become constant later on. Am I able to recompute the potential at a new time value, take, the previous wave function, and plug it into a new set of matrices ($A$, and $M$ in the article) for the next step?

Additional context: I see in the post there are terms $$V(x, y, z)\psi(x, y, z) = \frac{1}{2}\left(V_{ij}^{n+1}\psi_{ij}^{n+1} + V_{ij}^{n}\psi_{ij}^{n}\right)$$ $$a_{ij} = \left(1 + 2r_x + 2r_y + i\frac{\Delta t}{2}V_{ij}^{n+1}\right)$$ $$b_{ij} = \left(1 + 2r_x + 2r_y + i\frac{\Delta t}{2}V_{ij}^{n}\right)$$ Where $b_{ij}$ and $a_{ij}$ are used each step to solve $\psi(x, y, t)$, this is for a 2D scenario in the post so I think $z$ in the first equation may have been intended to be $t$ as $n$ is the time index and the term $\psi(x, y, t)V(x, y, t)$ appears earlier in the post. I think $a$ and $b$ may not have a time index because the only term in them that depends on time is the potential, but if you don't have a time-varying potential, they would remain the same values throughout each step. So here, the author I think, excludes it because they use a constant potential ($V_{ij}^{n} = V_{ij}^{n+1}$). But if you did not have a constant potential, could you simply put $n$ it back onto $a$ and $b$ and recompute $a_{ij}^{n+1}$ and $b_{ij}^{n}$ each step? And if so, what are the limits, are there limits to how the potential can change from step to step?

  • 1
    $\begingroup$ The Schrodinger equation is just a time evolution PDE, any method can be used to solve it (in principle). Of course different methods will have different constraints on the time step for numerical stability and different accuracy consideration. $\endgroup$ Feb 1, 2023 at 0:36
  • $\begingroup$ @MaximUmansky Thanks for the reply, so I assume the answer is yes? The reason I ask is I have had some trouble finding good methods to simulate it with a time varying potential. It is difficult to simulate with something like Euler's method or Runge-Kutta because it has no first derivative in space, so it cant be integrated into the second derivative or be fed back into $\psi$. Most of the time people seem not to vary potential, when I look up Time Dependent theory, I get either the case where we don't vary the potential in time, so the equation separates out, or I get perturbation theory $\endgroup$
    – cgbsu
    Feb 1, 2023 at 1:03
  • $\begingroup$ @MaximUmansky I used numerov before but now I need a dynamic solution, I dont know how to adapt the numerov method for that. $\endgroup$
    – cgbsu
    Feb 1, 2023 at 1:05
  • 2
    $\begingroup$ I don't understand the point about Euler's method or Runge-Kutta. I suspect that you have fundamentally understood how these methods work, and you are poking at more complicated methods such as the C-N method because you have not understood how Euler's method works. My suggestion is to post a question in which you describe what you perceive the problems with Euler's method to be, and let people point out what your misconceptions are. If you understand how Euler's method works, you will see that there is no fundamental difference between Euler's method, Crank-Nicolson, and Runge-Kutta. $\endgroup$ Feb 1, 2023 at 2:33
  • $\begingroup$ @WolfgangBangerth You could be right. AFIK Euler takes an initial numerical value of a function, with an analytically-known first-derivative, it computes the next step in time based on this derivative and the previous value of the function $y_{n+1} = y_n + hf(t_n, y_n)$. So it integrates the first derivative (and the previous value of the function) into the next value of the function. Runge-Kutta (R-K) is similar, it takes an analytically-known first-derivative and integrates it into approximations for higher order derivative terms, then uses those to predict the next value of the function 1/2 $\endgroup$
    – cgbsu
    Feb 4, 2023 at 21:30

2 Answers 2


I suspect that section 3.3 of the following thesis may be what you're looking for: Numerical simulation of time-dependent quantum systems.

Edit: As noted by "davidhigh" in another answer the method discussed in the thesis that I linked is NOT norm conserving. I have implemented the method for my own purposes and I have noticed that the norm of the resulting time-dependent wavefunction (where the norm is taken at each time step) varies from 0.9999999899578632 to 1.0002454052762269. However, the variation in the norm (relative to 1.0) decreases with decreasing time-step size. So although the method is not norm-conserving, it does seem to be convergent.



Quantum-mechanical time-evolution assigns to a given wavefunction at time $t_0$ a new one at time $t$ under preservation of the norm. Its action is thus expressed through a unitary time-evolution operator $\hat U(t,t_0)$ defined by the relation \begin{align} |{\Psi(t)}\rangle \ = \ \hat U(t,t_0) \, |\Psi(t_0)\rangle\,. \end{align} The equation of motion and the initial condition follow directly from the time-dependent Schrödinger equation, \begin{align} \label{eq:eom_time_evolution} i\partial_t \, \hat U(t,t_0) \ &= \ \hat H(t) \, \hat U(t,t_0)\,,\\ \label{eq:ic_time_evolution} \hat U(t_0,t_0) \ &= \ 1\,. \end{align} By iterating the integral of the previous equation, one can state the formal solution \begin{align} \label{eq:operator_time_evolution} \hat U(t,t_0) \ &= \ \sum_{n=0}^\infty U_n(t,t_0)\,,\\ U_n(t,t_0) \ &= \ (-i)^n \int_{t_0}^t dt_1 \int_{t_0}^{t_1} dt_2 \cdots \int_{t_0}^{t_n} dt_n \; \hat H(t_1) \cdots \hat H(t_n)\,, \end{align} which is known as the Dyson series. It may alternatively be written as \begin{align} \hat U(t,t_0) \ = \ \hat{\mathcal{ T}} \ \exp \left\{-i\int_{t_0}^t \hat H(\bar t)\, d\bar t\,\right\}\;. \end{align} where the time-ordering operator $\hat {\mathcal{T}}$ orders later times to the left.

The direct use of the Dyson series to determine $U(t,t_0)$ is tedious because truncation at a given order does not lead to a unitary operator.

There are basically two ways to handle the Dyson series -- either account for or neglect the time-ordering. An example of the first approach is the Magnus expansion. On the other hand, there are the more frequently encountered methods which neglect the time-ordering. Their foundation is given by the following piecewise-in-time approximation to the Hamiltonian, \begin{align} \hat U(t,t_0) \ &= \ \hat{\mathcal{ T}} \ \exp \left\{-i\int_{t_0}^t \hat H(\bar t)\, d\bar t\,\right\}\\ &\approx \ \hat{\mathcal{ T}} \ \exp \left\{-i \sum_{j=1}^{N_t}\hat H(t_j) \Delta t \,\right\}\\ \label{eq:time_evol_approx_3} &\approx \ \hat{\mathcal{ T}} \ \exp \left\{-i \hat H(\bar t_{N_t})\Delta t\right\} \cdots \exp \left\{-i \hat H(\bar t_{1})\Delta t\right\}\\[0.5em] &= \ \exp \left\{-i \hat H(\bar t_{N_t})\Delta t\right\} \cdots \exp \left\{-i \hat H(\bar t_{1})\Delta t\right\}\,. \end{align} Here, the integral is evaluated via the trapezoidal rule using $N_t$ intermediate time points, and subsequently the exponential is factorized under neglect of commutators of higher order in $\Delta t$. As the resulting term is already time-ordered, $\hat{\mathcal{ T}}$ can be trivially applied.

This also answers the question in the title: Any method for the time-propagation of stanionary Hamiltonians can also be applied for the solution of time-dependent Hamiltonians, but this always involves an approximation.

Crank-Nicolson propagator

There seem to be two variants of the Crank-Nicolson propagator.

  1. The one which is common to me takes the temporal separation as shown before as a basis, and then uses an approximation of the time-evolution operator in terms of the Caley operator: \begin{align} \exp \left\{-i \hat H(t_{j})\Delta t\right\} \approx \frac{1+\frac{i}2 H(t_{j})\Delta t}{1-\frac{i}2 H(t_{j})\Delta t} \end{align} The good thing about this approximation is that it is unitary. The not-so-good thing is that it doesn't account for time-ordering (as shown before).

    Using this aproximated form of the time-evolution operator, the CN-method then proceeds to solve the linear system \begin{align} \bigg(1-\frac{i}2 H(t_{j})\Delta t\bigg)|\Psi(t_{j+1})\rangle = \bigg(1+\frac{i}2 H(t_{j})\Delta t\bigg)|\Psi(t_{j})\rangle\,, \end{align} where the right-hand side is known (resp. can be evaluated at time $t_j$).

  2. The second type of CN-propagator is used in the link of the other answer. It basically uses the Hamiltonian at the later time $t_{j} + \Delta t$ for back-propagation: \begin{align} \frac{1+\frac{i}2 H(t_{j})\Delta t}{1-\frac{i}2 H(t_{j} + \Delta t)\Delta t} \end{align} I don't know whether this methods is theoretically justified. It might be a bit better regarding time-ordering, however, in general it is not unitary.

    EDIT: Wikipedia lists this second alternative as the standard CN-method, and justifies it via a Forward and Backward Euler step. As the Hamiltonian is evaluated at the right-boundary time, it is further possible to interpret the CN-method as an implicit Runge-Kutta method.

Summarizing: In quantum-mechanical applications, I would use the CN-propagator of alternative 1., by which you obtain norm-conservatuon. However, one has to accept that the time-ordering is only approximately incorporated, and, as usual, reduce the time-step to attain greater accuracy.


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.