# 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

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

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