Time evolution of local green function using two applyMPO for spinless-fermions.

Question

Hi,
I am trying to calculate the local green function Gii(t)= <\psi0|Cdagi(t)Ci|\psi0> at site i=N, using two applyMPO for one-dimensional Kitaev model [equation 1 of PRB 88, 161103(R) (2013)]. For this I am doing following steps:
(i) I am calculating |\phi>=Ci |\psi0> by applying local operator Ci to the ground state |psi0>.
(ii) Time evolution of |\phi> using applyMPO: |phi(t)>=exp(-iHt)|\phi>.
(iii) Application of local operator Cdagi to |\phi(t): Cdagi|phi(t)>.
(iv)Time evolution of |psi0> using applyMPO: |psi(t)>=exp(-iHt)|\psi_0>.
(v) Overlap of <psi(t)|phi(t)> using innerC.

After running the code, I am getting non-zero values for the odd time steps and zero for even number of time steps:
0.001 -0.000605309 -0.671164
0.002 0 -0
0.003 0.00300302 0.671158
0.004 0 -0
0.005 -0.00540069 -0.671143

Here is the code:

`//-------------Creating |phi> = C_i1|Psig>--------

 auto i1=N;  // for site i1=N. 
 psig.position(i1);
 auto newpsi = noPrime(psig(i1)*op(sites,"A",i1));
 psig.set(N, newpsi);

//--------------Jordan-Wigner string----------------------

 for(int k = i1-1; k >=1; k--)
{
 psig.position(k);
auto newpsi1 = noPrime(psig(k)*op(sites,"F",k)); 
   psig.set(k, newpsi1);
}
  psig.noPrime().normalize();

auto tau=0.001;
auto ii = Complexi;
auto args = Args("Method=","DensityMatrix","Cutoff=",1E-14,"MaxDim=",7000);
auto expH = toExpH(ampo,tstep*Cplxi);
//-------------Time Evolution--------------------------

auto nt = int(ttotal/tau+(1e-9*(ttotal/tau)));
for(int n = 1; n <= nt;++n)
{
           psig = applyMPO(expH,psig,args);     
           psig.noPrime().normalize();

//--Jordan-Wigner string

 for(int k = 1; k <i1; ++k)
  {  
       psig.position(k);
       auto newpsi2 = noPrime(psig(k)*op(sites,"F",k)); 
       psig.set(k, newpsi2);
}

// ----------------- Cdag|phi(t)>--------------------

    psig.position(i1);
    auto newpsi3 = noPrime(psig(i1)*op(sites,"Adag",i1));
    psig.set(i1, newpsi3);
    psig.noPrime().normalize();

//------------------psi(t)>=exp(-iHt)|psig>-----------------------

      psi = applyMPO(expH,psi,args);
      psi.noPrime().normalize();
      auto result2 = -ii*innerC(psi,psig);
      file2<<float(n*tstep)<<' '<< result2.real()<<' ' <<result2.imag() <<std::endl;

Answer 1

Hi, thanks for the good question and your patience with a slow reply.

While I'm not sure why you getting the value of zero on every other time step, I think there are a few major improvements to your approach which you should consider anyway, and after you make them this problem may disappear in the process of changing your code.

Here are the improvements or fixes:
1. since one of the exp(-i Ht) operators acts on the ground state of H, you can rigorously replace it by a scalar exp(-i E0 t) where E0 is the ground state energy <psi0|H|psi0> and so skip that step

1. a much better choice for doing the time evolution would be the TDVP algorithm, for which a high-quality implementation in ITensor is now available here:
   https://github.com/ITensor/TDVP
   including a helpful sample code

Especially #2 above might eliminate the zeros because I and others have found that unfortunately the MPO method (based on toExpH) has not turned out to be as accurate as originally hoped (I mean as an algorithm - not that there is anything buggy about how it is implemented by ITensor). TDVP is basically just a superior method in almost every way, with one caveat that one has to be careful about growing or adapting the MPS basis within TDVP. The 2-site TDVP method does this automatically but is somewhat less accurate than the 1-site TDVP method. For the 1-site TDVP method, the basis adaptation isn't automatic but can be achieved through a basis extension method (the addBasis function) shown in the TDVP code example.

Or, last but very much not least, you can just use the Trotter gate or TEBD method for your time evolution. For your model which is short ranged and 1D, TEBD is a great method. It is efficient and automatically adapts the bond dimension in a very 'foolproof' and reliable way. The only advantage of TDVP over Trotter/TEBD for you model would be that TDVP allows larger time steps, basically. If I were in your shoes, I would implement both and use one as a check on the other, and not use the toExpH method.

Best regards,
Miles

Comments

Hi Miles,  thank you for your detailed answer. Yes, the code is working fine with exp(-i E0 t)  without giving zero in alternative steps. My primary motivation is to calculate the  local green function for a non-equilibrium state (where \psi(t) is not a ground state of H :
G(t,t')=〈ψ(t)|c†_i exp(iHt′)c_j exp(−iHt′)|ψ(t)〉and |ψ(t)〉=exp(−iHt)|ψ(0)). So that we have time-dependent local-density of state, after Fourier transform with respect to t′ of the correlator G(t,t'). For that reason, I need two-time evolution in one time step. Thank you for your suggestion to use TEBD/TDVP.  I will try with TEBD and TDVP for two-time evolution process.
With regards
Bradraj Pandey