Integration scheme

For a detailed review of algorithms, please refer to W-SLDA Toolkit: A simulation platform for ultracold Fermi gases. Here we provide only basic information, and illustrate it with the AB4AM5 algorithm.

Time-dependent equations

The td-wslda code solves equations of the form:

i \dot{\varphi_n}(t) = \hat{H}(\{\varphi_n\},t)\varphi_n(t)

where the Hamiltonian is given by

\hat{H} = 
  \begin{pmatrix}
    h_{\uparrow}(\bm{r},t) & \Delta(\bm{r},t) \\
    \Delta^*(\bm{r},t)& -h_{\downarrow}^*(\bm{r},t) 
  \end{pmatrix}

and \varphi_n is a shorthand notation for quasiparticle wavefunctions (qpwfs):

\varphi_n(\bm{r},t)=
\begin{pmatrix}
    u_{n,\uparrow}(\bm{r},t) \\ 
    v_{n,\downarrow}(\bm{r},t)
  \end{pmatrix}

Note that all qpwfs fluctuate in time with a typical oscillating factor \exp(-i E_n t), where E_n is a quasi-particle energy. Since we evolve all states with quasi-particle energies from the interval [-E_c,E_c], the evolved states exhibit both very slow as well as very rapid oscillations in time. In order to improve the accuracy of the time integration, in each integration step, we subtract from the quasi-particle solution the typical frequency of oscillations, which in practice modifies the equation to the form:

i \dot{\varphi_n}(t) = (\hat{H}(\{\varphi_n\},t)-\langle H\rangle_n(t))\varphi_n(t)

where \langle H\rangle_n is instantaneous quasi-particle energy. Although this step incurs a high numerical cost in the calculations, it greatly improves the accuracy of the time integration.

ABM integration algorithm

We convert time-dependent equations into form:

 \dot{y}(t) = f(y,t)

where:

y(t) \equiv \varphi_n(\bm{r},t)

f(y,t) \equiv \dfrac{1}{i}(\hat{H}(\{\varphi_n\},t)-\langle H\rangle_n)\varphi_n(\bm{r},t)

Adams-Bashforth-Moulton (ABM) is designed as follows:

predictor, for example, of 4th order:

y_k^{(p)} = y_{k-1} + \dfrac{55}{24}\Delta t f_{k-1} - \dfrac{59}{24}\Delta t f_{k-2} + \dfrac{37}{24}\Delta t f_{k-3} - \dfrac{9}{24}\Delta t f_{k-4}

corrector, for example, of 5-th order:

y_k = y_{k-1} + \dfrac{251}{720}\Delta t f(y_k^{(p)},k\Delta t) +\dfrac{646}{720}\Delta t f_{k-1}-\dfrac{264}{720}\Delta t f_{k-2}+\dfrac{106}{720}\Delta t f_{k-3} - \dfrac{19}{720}\Delta t f_{k-4}

where f_{k}=f(y_k, k\Delta t). The accuracy of the method is set by the corrector.

The algorithm requires 5 (wave function) buffers called: y_{k-1}, f_{k-1}, f_{k-2}, f_{k-3}, f_{k-4}.

The algorithm is as follows (operations listed within a single step have to be executed simultaneously):

  \left\lbrace 
  \begin{array}{lll}
   y_{k-1} & \leftarrow & y_{k-1} + \dfrac{55}{24}\Delta t f_{k-1} - \dfrac{59}{24}\Delta t f_{k-2} + \dfrac{37}{24}\Delta t f_{k-3} - \dfrac{9}{24}\Delta t f_{k-4},\\
   f_{k-4} & \leftarrow & y_{k-1} +\dfrac{646}{720}\Delta t f_{k-1}-\dfrac{264}{720}\Delta t f_{k-2}+\dfrac{106}{720}\Delta t f_{k-3} - \dfrac{19}{720}\Delta t f_{k-4}
  \end{array}
  \right.

Compute densities from y_{k-1} and formulate potentials (like U(\bm{r}) and \Delta(\bm{r})) needed to formulate Hamiltonian. This step requires MPI communication. For densities and potentials, use additional buffers; note that these are small compared to wave-function buffers.
y_{k-1} \leftarrow f(y_{k-1},k\Delta t), this is the first application of Hamiltonian; use potentials formulated in step 2. This step can be done in batch mode. This is a numerically intensive step.
y_{k-1} \leftarrow \dfrac{251}{720}\Delta t y_{k-1} + f_{k-4}, this step can be merged with step 3. Now buffer y_{k-1} holds wave functions for t=k\Delta t.
Compute densities from y_{k-1} and formulate potentials (like U(\bm{r}) and \Delta(\bm{r})) needed to formulate Hamiltonian. This step requires MPI communication. For densities and potentials, use additional buffers; note that these are small compared to wave-function buffers.
Preparation of buffers for the next step. Note that the last operation corresponds to the second application of the Hamiltonian; use potentials formulated in step 5. This step can be done in batch mode. This is a numerically intensive step.

  \left\lbrace 
  \begin{array}{lll}
   f_{k-4} & \leftarrow & f_{k-3},\\
   f_{k-3} & \leftarrow & f_{k-2},\\
   f_{k-2} & \leftarrow & f_{k-1},\\
   f_{k-1} & \leftarrow & f(y_{k-1},k\Delta t).
  \end{array}
  \right.

Increase k by one and go to step 1.