Time Evolving an MPS with Trotter Gates

A very accurate, efficient, and simple way to time evolve a matrix product state (MPS) is by using a Trotter decomposition of the time evolution operator.

Although the discussion below focuses on the case of a nearest-neighbor one-dimensional Hamiltonian, the method can be extended to Hamiltonians with arbitrary finite-range interactions by using swap gates to temporarily exchange sites. (For information about using swap gates, see New J. Phys. 12, 055026.)

If the Hamiltonian is a sum of local terms

$$ H = \sum_j h_{j,j+1} $$

where @@h_{j,j+1}@@ only acts non-trivially on sites j and (j+1), then a Trotter decomposition that is particularly well suited for use with MPS techniques is

$$ e^{-i \tau H} \approx e^{-i h_{1,2} \tau/2} e^{-i h_{2,3} \tau/2} \cdots e^{-i h_{N-1,N} \tau/2} e^{-i h_{N-1,N} \tau/2} e^{-i h_{N-2,N-1} \tau/2} \cdots e^{-i h_{1,2} \tau/2} + O(\tau^3) $$

Note the factors of two in each exponential. The error in the above decomposition is of order @@\tau^3@@ , so this will be the error accumulated per time step. Because of the time-step error, one takes @@\tau@@ to be small and then applies the above set of operators to an MPS as a single sweep, then does a number @@(t/\tau)@@ of sweeps to evolve for a total time @@t@@ . The total error will therefore scale as @@\tau^2@@ with this scheme, though other sources of error may dominate for long times, or very small @@\tau@@ , such as truncation errors.

The same decomposition can be used for imaginary time evolution just by replacing @@i \tau \rightarrow \tau@@ .

Below is a fully working code that applies the above ideas to time evolve an MPS which is initially a simple product state with the Heisenberg Hamiltonian.

The BondGate class has a constructor that accepts a local Hamiltonian term as a tensor, and automatically exponentiates it to make a Trotter gate with the specified time step. The gateTEvol function is a helper function that handles the logic of applying a container of gates to an MPS the right number of times and with the appropriate truncation of the MPS after each step.

#include "itensor/all.h"
#include "itensor/util/print_macro.h"

using namespace itensor;
using std::vector;

int main()
{
int N = 50; //number of sites
Real tstep = 0.02; //time step (smaller is generally more accurate)
Real ttotal = 1.0; //total time to evolve
Real cutoff = 1E-8; //truncation error cutoff when restoring MPS form

//Define a site set object "sites" which lets us
//easily obtain Site indices defining our Hilbert space
//and S=1/2 single-site operators
auto sites = SpinHalf(N);

//Make initial MPS psi to be in the Neel state
auto state = InitState(sites);
for(auto j : range1(N))
    {
    state.set(j,j%2==1?"Up":"Dn");
    }
auto psi = MPS(state);

//Create a std::vector (dynamically sizeable array)
//to hold the Trotter gates
auto gates = vector<BondGate>();

//Create the gates exp(-i*tstep/2*hterm)
//and add them to gates
for(int b = 1; b <= N-1; ++b)
    {
    auto hterm = op(sites,"Sz",b)*op(sites,"Sz",b+1);
    hterm += 0.5*op(sites,"S+",b)*op(sites,"S-",b+1);
    hterm += 0.5*op(sites,"S-",b)*op(sites,"S+",b+1);

    auto g = BondGate(sites,b,b+1,BondGate::tReal,tstep/2.,hterm);
    gates.push_back(g);
    }
//Create the gates exp(-i*tstep/2*hterm) in reverse
//order (to get a second order Trotter breakup which
//does a time step of "tstep") and add them to gates
for(int b = N-1; b >= 1; --b)
    {
    auto hterm = op(sites,"Sz",b)*op(sites,"Sz",b+1);
    hterm += 0.5*op(sites,"S+",b)*op(sites,"S-",b+1);
    hterm += 0.5*op(sites,"S-",b)*op(sites,"S+",b+1);

    auto g = BondGate(sites,b,b+1,BondGate::tReal,tstep/2.,hterm);
    gates.push_back(g);
    }

//Save initial state;
auto psi0 = psi;

//Time evolve, overwriting psi when done
gateTEvol(gates,ttotal,tstep,psi,{"Cutoff=",cutoff,"Verbose=",true});

printfln("Maximum MPS bond dimension after time evolution is %d",maxLinkDim(psi));

//Print overlap of final state with initial state
//(Will be complex so using innerC which can return complex);
Print(innerC(psi,psi0));

return 0;
}

Download the full example code

Back to Formulas

Back to Main