Skip to content

time_evolve_expokit

Computes the real-time evolution,

\[\vert \psi(t) \rangle = e^{-iHt} \vert \psi_0\rangle,\]

of a State \(\vert \psi_0 \rangle\) and a Hermitian operator \(H\) using the iterative algorithm implemented by Expokit

Expokit: A Software Package for Computing Matrix Exponentials
Roger B. Sidje
ACM Trans. Math. Softw., 24(1):130-156, 1998. (1998)
DOI: 10.1145/285861.285868

The algorithm features automatic stepsize control and computes approximate solutions with high precision according to our tests. Yet, the evolve_lanczos implementation is currently faster and more memory efficient.

Sources
time_evolve_expokit.hpp
time_evolve_expokit.cpp
time_evolve_expokit.jl

Definition

The method is provided in two variants:

  1. Returning a new state while the input state remains untouched. This variant is safe to use and simple to code.

    TimeEvolveExpokitResult time_evolve_expokit(
    OpSum const &ops, State state, double time, double precision = 1e-12,
    int64_t m = 30, double anorm = 0., int64_t nnorm = 2);

    time_evolve_expokit(ops::OpSum, state::State, time::Float64;
                    precision::Float64=1e-12, m::Int64 = 30, 
                    anorm::Float64 = 0.0, nnorm::Int64 = 2)

  2. An inplace variant time_evolve_expokit_inplace, where the input state is overwritten and contains the time evolved state upon exit. This version is more memory efficient than time_evolve_expokit.

    TimeEvolveExpokitInplaceResult time_evolve_expokit_inplace(
    OpSum const &ops, State &state, double time, double precision = 1e-12,
    int64_t m = 30, double anorm = 0., int64_t nnorm = 2);

    time_evolve_expokit_inplace(ops::OpSum, state::State, time::Float64;
                            precision::Float64=1e-12, m::Int64 = 30, 
                            anorm::Float64 = 0.0, nnorm::Int64 = 2)

Parameters

Name Description Default
H OpSum defining the hermitian operator \(H\) for time evolution
psi0 initial State \(\vert \psi_0 \rangle\) of the time evolution
time time \(t\) until which the state is evolved
precision accuracy of the computed time evolved state 1e-12
m dimension of used Krylov space, main memory requirement 30
anorm 1-norm estimate of the operator \(H\), if unknown default 0. computes it fresh 0.
nnorm number of random samples to estimate 1-norm, usually not more than 2 required 2

Returns

A struct with the following entries

Entry Description
error the computed error estimate during evolution
hump the "hump" as defined in Expokit 10.1145/285861.285868
state time-evolved State \(\vert \psi(t)\rangle\) (not defined for inplace variant)

Convergence criterion

The algorithm is estimating the following error,

\[ \varepsilon = \parallel \vert \tilde{\psi}(t)\rangle - e^{z(H - \delta)} \vert \psi_0\rangle \parallel_2, \]

where \(\vert \tilde{\psi}(t) \rangle\) denotes the approximation computed during the algorithm. As the exact solution is not available this error is estimated using the method described by Algorithm 2 in

Expokit: A Software Package for Computing Matrix Exponentials
Roger B. Sidje
ACM Trans. Math. Softw., 24(1):130-156, 1998. (1998)
DOI: 10.1145/285861.285868

Usage Example

int N = 8;
int nup = N / 2;
auto block = Spinhalf(N, nup);

// Define the nearest-neighbor Heisenberg model
auto ops = OpSum();
for (int i=0; i<N; ++i) {
  ops += Op("SdotS", {i, (i+1) % N});
}

auto psi0 = product_state(block, {"Up", "Dn", "Up", "Dn", "Up", "Dn", "Up", "Dn"});
double time = 1.0;
double precision = 1e-8;
auto res1 = time_evolve_expokit(ops, psi0, time, precision);
auto res2 = time_evolve_expokit_inplace(ops, psi0, time, precision);
XDIAG_SHOW(isapprox(psi0, res1.state));
XDIAG_SHOW(res1.error);
XDIAG_SHOW(res1.hump);

let
    N = 8
    nup = N รท 2
    block = Spinhalf(N, nup)

    # Define the nearest-neighbor Heisenberg model
    ops = OpSum()
    for i in 1:N
        ops += Op("SdotS", [i, mod1(i+1, N)])
    end

    psi0 = product_state(block, ["Up", "Dn", "Up", "Dn", "Up", "Dn", "Up", "Dn"])
    time = 1.0
    res1 = time_evolve_expokit(ops, psi0, time, precision=1e-8)
    res2 = time_evolve_expokit_inplace(ops, psi0, time, precision=1e-8)
    @show isapprox(psi0, res1.state)
    @show res1.error
    @show res1.hump
end