Hi,

I'm running tDMRG for a simple anderson impurity coupled to two semi-infinite chain at different chemical potential, but get different dynamics for different spin states. Details are following:

First in the DMRG part, I decouple left chain, impurity, right chain to calculate the ground state of the decoupled system. Left chain and right chain are both half filled with same number of spin-up and spin-down electrons, and impurity is empty. My DMRG results are correct compared with analytic results (for a finite chain with even number of sites, the ground state of half filling has exactly local occupation of 0.5 for both spin-up and spin-down states.)

Secondly, I quench the whole system by coupling the impurity site and chains, and measure the dynamics of local occupation of spin-up and spin-down states. But these two occupations deviate very fast, which is unphysical considering that the Hamiltonian is symmetric under spin inversion). See the following results

t=0.01, Nup = 0.000945655, Ndn = 0.000945656

t=0.02, Nup = 0.00370455, Ndn = 0.0037134

t=0.03, Nup = 0.00812654, Ndn = 0.00818666

t=0.04, Nup = 0.0140243, Ndn = 0.0142351

t=0.05, Nup = 0.0211824, Ndn = 0.0217184

t=0.06, Nup = 0.0293665, Ndn = 0.03049

t=0.07, Nup = 0.0383324, Ndn = 0.0404003

t=0.08, Nup = 0.0478346, Ndn = 0.0512999

t=0.09, Nup = 0.057634, Ndn = 0.0630419

t=0.10, Nup = 0.0675047, Ndn = 0.0754844

My code is:

#include"itensor/all.h"

#include <math.h>

#include

using namespace itensor;

```
int main(int argc, char *argv[])
{
float U=4;
float e0=-2;
float tb=10;
float td=3.1623;
float V=4;
int L=101;
double tau=0.01;
float t_tot=3;
auto Nc=L/2+1;//center site, impurity
double el=V/2;
double er=-V/2;
//Read ground-state
Electron sites;
readFromFile("site_file",sites);
MPS psi(sites);
readFromFile("ground_state",psi);
auto ampo = AutoMPO(sites);
//Right Chain
for(int j=Nc+1;j <=L-1; j+=1)
{
ampo += -tb,"Cdagup",j,"Cup",j+1;
ampo += -tb,"Cdagup",j+1,"Cup",j;
ampo += -tb,"Cdagdn",j,"Cdn",j+1;
ampo += -tb,"Cdagdn",j+1,"Cdn",j;
ampo += er,"Ntot",j;
}
ampo+=er,"Ntot",L;
//Left chain
for(int j=1;j < Nc-1; j+=1)
{
ampo += -tb,"Cdagup",j,"Cup",j+1;
ampo += -tb,"Cdagup",j+1,"Cup",j;
ampo += -tb,"Cdagdn",j,"Cdn",j+1;
ampo += -tb,"Cdagdn",j+1,"Cdn",j;
ampo += el,"Ntot",j;
}
ampo+=el,"Ntot",Nc-1;
//Impurity site
ampo+=e0,"Ntot",Nc;
ampo+=U,"Nup",Nc,"Ndn",Nc;
//Coupling
ampo+= -td,"Cdagup",Nc-1,"Cup",Nc;
ampo+= -td,"Cdagup",Nc,"Cup",Nc-1;
ampo+= -td,"Cdagdn",Nc-1,"Cdn",Nc;
ampo+= -td,"Cdagdn",Nc,"Cdn",Nc-1;
ampo+= -td,"Cdagup",Nc+1,"Cup",Nc;
ampo+= -td,"Cdagup",Nc,"Cup",Nc+1;
ampo+= -td,"Cdagdn",Nc+1,"Cdn",Nc;
ampo+= -td,"Cdagdn",Nc,"Cdn",Nc+1;
auto nt=round(t_tot/tau);
auto expH=toExpH(ampo,tau*Cplx_i);
auto args=Args("Method=","DensityMatrix","Cutoff=",1E-9,"MaxDim=",1000);
Real Ild, Ilu, Ird, Iru, nu, nd;
for(int j=1;j<=nt;j++){
psi=applyMPO(expH,psi,args);
psi.noPrime().normalize();
printf("t=%g\n",tau*j);
psi.position(Nc);
//Occu
nu=std::real(eltC(psi(Nc)*Nu*dag(prime(psi(Nc),"Site"))));
nd=std::real(eltC(psi(Nc)*Nd*dag(prime(psi(Nc),"Site"))));
printf("Nu=%g, Nd=%g, ",nu,nd);
return 0;
}
```

Can anyone help? THANK YOU!