-
Notifications
You must be signed in to change notification settings - Fork 0
Temporal_propagation
One can solve the time-dependent Schroedinger or Gross-Pitaevskii equation.
One starts with an arbitrary initial state, which can be anything: a Gaussian wavepacket, a plane wave, a state localized on a single site, or whatever the user builds.
There are two methods to solve the evolution equation: the slow ODE method and the fast Chebyshev method.
The discretized Schroedinger or Gross-Pitaevskii equation is a set of coupled differential equations between the various \psi_i(t), which can be solved using a standard Python package for ODEs. The advantage is that it works for arbitrary interaction strengths in the GPE. The drawback is that it is slow. Especially for the Schroedinger equation, the equations are linear, but this is not used by the algorithm.
In the non-interacting case, the evolution operator over time t is exp(-iHt) so that:
\psi(t) = \exp{-iHt}\psi(t=0)
The evolution operator is never explicitly built, but it can formally be written as an infinite series of Chebyshev polynomials of H, see Kernel Polynomial Method. The coefficients of theis series decay with the order, so that the series can be efficiently propagated. Thanks to the recurrence relations between Chebyshed polynomials, it is possible to build recursively \psi(t). If the series is truncated at order N, the calculation involves N matrix vector products of the type H|\phi>. Because of the sparse structure of H, this can be computed efficiently. The number of terms in the series is roughly range(H)*t, where range(H)=Emax_Emin is the spectral range of the Hamiltonian. This is valid when N is large, for smaller range(H)*t, N has to be taken slightly larger. As the CPU time needed scales like N, the cost of propagating over a long time (by chaining many elementary Chebyshev propagations) is almost independent of the time step, provided range(H)*t >> 1. A basic rule of thumb is not to take time steps making N smaller than 10.
The interacting case is slightly more complicated, as the non-linear term is time dependent and thus cannot be included in the Chebyshev expansion. The program uses a split-step algorithm by interleaving elementary propagation operators \exp{-iHdt} and \exp{-ig|\psi|^2dt}. Propagation with \exp{-iHdt} uses the previously described Chebyshev expansion while the non-linear propagation \exp{-i g|\psi|^2dt} is diagonal in configuration space. Of course, this is a good approximation only if g is sufficiently small. Hence, the larger g, the smaller the time step dt to be taken. This may be a drawback because the number N of terms in the Chebyshev expansion may become small, implying a loss of efficiency. The program records the maximum value of g|\psi|^2dt reached during the whole propagation (the so-called maximum nonlinear phase), and prints it. A rule of thumb for systems where the disorder is stronger than interaction is that it should be smaller than 0.1. Sometimes, it has been observed that smaller values, say 0.02-0.05, are needed for good convergence. When the interaction is larger than the disorder, the density |\psi|^2 is more uniform and larger values of the nonlinear phase can be used. In any case, it is always good to check the conservation of the total energy, setting "dispersion_energy" to True in the parameter file.