Different results for spin up and down

Question

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!

Comments

It looks like the occupations match up to a certain time. Perhaps this is due to truncation or timestep error? Have you tried decreasing the timestep or the cutoff?