Hi,

I am trying to calculate the local green function G*ii(t)= <\psi*0|Cdag*i(t)C*i|\psi*0> 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>=C*i |\psi

*0> by applying local operator C*i to the ground state |psi

*0>.*

(ii) Time evolution of |\phi> using applyMPO: |phi(t)>=exp(-iHt)|\phi>.

(iii) Application of local operator Cdagi to |\phi(t): Cdag

(ii) Time evolution of |\phi> using applyMPO: |phi(t)>=exp(-iHt)|\phi>.

(iii) Application of local operator Cdag

*i|phi(t)>.*

(iv)Time evolution of |psi0> using applyMPO: |psi(t)>=exp(-iHt)|\psi_0>.

(iv)Time evolution of |psi

(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 = Complex*i;
auto args = Args("Method=","DensityMatrix","Cutoff=",1E-14,"MaxDim=",7000);
auto expH = toExpH(ampo,tstep*Cplx*i);

//-------------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;
```