Flavor mixing in ITensor

Hello Miles,

I recently run into a problem with flavor mixing. More specifically, I have fermion c_1 with charge q_1=3 hopping on odd sites and the other fermion c_2 with charge q_2=1 hopping on even sites. My U(1) invariant interaction will eventually mix the two flavors. However, at the benchmark stage, when there is no interaction and when the system size is relatively big (~20 unit cells), flavor mixing happens, meaning that three c_2 fermions disappear from the even sites and become one c_1 fermion on the odd sites, even if my initial state only consists of c_2 fermions. Even though this mixing still conserves total U(1) charge and fermion parity, ideally it’s not supposed to happen because the two flavors should be hopping independently before interaction is added.

Here is how I define the two site types:

function ITensors.space(
  ::SiteType"Odd";
  conserve_qns=false,
  conserve_Q=conserve_qns,
  conserve_nfparity=conserve_qns,
  qnname_Q="Q",
  qnname_nfparity="NfParity",
 )   
  if conserve_Q && conserve_nfparity
    return [
      # 0 particle
      QN((qnname_Q, 0),(qnname_nfparity,0,-2)) => 1,
      # 1 particle
      QN((qnname_Q, +3),(qnname_nfparity,1,-2)) => 1,
    ]
  end
  return 2
end

function ITensors.space(
  ::SiteType"Even";
  conserve_qns=false,
  conserve_Q=conserve_qns,
  conserve_nfparity=conserve_qns,
  qnname_Q="Q",
  qnname_nfparity="NfParity",
 )   
  if conserve_Q && conserve_nfparity
    return [
      # 0 particle
      QN((qnname_Q, 0),(qnname_nfparity,0,-2)) => 1,
      # 1 particle
      QN((qnname_Q, +1),(qnname_nfparity,1,-2)) => 1,
    ]
  end
  return 2
end

I have the following two questions:

  1. Is there a way to avoid the mixing at free fermion level?
  2. If the answer to the 1st question is yes, then would it prevent me from adding flavor-mixing interactions?

Thanks a lot for your time.
-Meng

Hi Meng,
Thanks for the question. Your space definitions look good to me.

Based on what you’re saying, the thing that is leading to the mixing is the numerical procedure you are doing afterward. But what is the procedure? Is it a DMRG calculation or a non-interacting calculation or something else? Thanks –

Sorry I should’ve made this clear. I did DMRG calculation using a simple hopping model for the benchmark stage: c_1 only hops on even sites and c_2 on odd sites, with the initial state only consisting of c_2 fermions. Thanks.

Here is a minimal version of my code:

using ITensors, HDF5, JLD, DelimitedFiles

function ITensors.space(
  ::SiteType"Odd";
  conserve_qns=false,
  conserve_Q=conserve_qns,
  conserve_nfparity=conserve_qns,
  qnname_Q="Q",
  qnname_nfparity="NfParity",
 )   
  if conserve_Q && conserve_nfparity
    return [
      # 0 particle
      QN((qnname_Q, 0),(qnname_nfparity,0,-2)) => 1,
      # 1 particle
      QN((qnname_Q, +3),(qnname_nfparity,1,-2)) => 1,
    ]
  end
  return 2
end

function ITensors.space(
  ::SiteType"Even";
  conserve_qns=false,
  conserve_Q=conserve_qns,
  conserve_nfparity=conserve_qns,
  qnname_Q="Q",
  qnname_nfparity="NfParity",
 )   
  if conserve_Q && conserve_nfparity
    return [
      # 0 particle
      QN((qnname_Q, 0),(qnname_nfparity,0,-2)) => 1,
      # 1 particle
      QN((qnname_Q, +1),(qnname_nfparity,1,-2)) => 1,
    ]
  end
  return 2
end


ITensors.val(::StateName"Emp",::SiteType"Odd" ) = 1
ITensors.val(::StateName"f1",::SiteType"Odd" ) = 2
ITensors.val(::StateName"Emp",::SiteType"Even" ) = 1
ITensors.val(::StateName"f2",::SiteType"Even" ) = 2

ITensors.state(::StateName"Emp",::SiteType"Odd" ) = [1.0 0.0]
ITensors.state(::StateName"f1",::SiteType"Odd" )  = [0.0 1.0]
ITensors.state(::StateName"Emp",::SiteType"Even" ) = [1.0 0.0]
ITensors.state(::StateName"f2",::SiteType"Even" )  = [0.0 1.0]


function ITensors.op!(Op::ITensor, ::OpName"C1", ::SiteType"Odd", s::Index)
  return Op[s' => 1, s => 2] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"C1d", ::SiteType"Odd", s::Index)
  return Op[s' => 2, s => 1] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"C2", ::SiteType"Even", s::Index)
  return Op[s' => 1, s => 2] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"C2d", ::SiteType"Even", s::Index)
  return Op[s' => 2, s => 1] = 1.0
end

function ITensors.op!(Op::ITensor, ::OpName"F", ::SiteType"Odd", s::Index)
  Op[s'=>1,s=>1] = +1.0
  Op[s'=>2,s=>2] = -1.0
end

function ITensors.op!(Op::ITensor, ::OpName"F", ::SiteType"Even", s::Index)
  Op[s'=>1,s=>1] = +1.0
  Op[s'=>2,s=>2] = -1.0
end


ITensors.has_fermion_string(::OpName"C1", ::SiteType"Odd") = true
ITensors.has_fermion_string(::OpName"C1d", ::SiteType"Odd") = true
ITensors.has_fermion_string(::OpName"C2", ::SiteType"Even") = true
ITensors.has_fermion_string(::OpName"C2d", ::SiteType"Even") = true

let
        N = 12
        t1 = 1
        t2 = -1
        sweeps = Sweeps(10)
        maxdim!(sweeps,10,100)
        cutoff!(sweeps,1e-10)
        noise!(sweeps,1E-5)
        
	function site_type(n::Int)
       		if mod(n,2)==1
         		return "Odd"
         	end
         	if mod(n,2)==0
         		return "Even"
         	end
	end

	sites = siteinds(n -> site_type(n), 2*N; conserve_qns=true)

        states = ["Emp" for i=1:2*N]       
        for i = 1:floor(Int,3*N/4)
           states[2*i-1] = "Emp"
           states[2*i] = "f2"
        end

        ampo = AutoMPO()
        
        for i = 1:(N-1)
            ########## f1 ##########################
           ampo += t1, "C1", 2*i-1, "C1d", 2*i+1
           ampo += t1, "C1", 2*i+1, "C1d", 2*i-1
            ########## f2 ##########################
           ampo += t2, "C2", 2*i, "C2d", 2*i+2
           ampo += t2, "C2", 2*i+2, "C2d", 2*i
        end
     
        H = MPO(ampo,sites)
        using HDF5,JLD,DelimitedFiles
        psi0 = productMPS(sites,states)
        energy0,psi = dmrg(H,psi0,sweeps)
        
        ######Check GS particle occupation on each site        
        for i = 1:N
        	numb1 = AutoMPO()
        	numb1 += "C1d",2*i-1, "C1", 2*i-1
        	N1 = MPO(numb1,sites)
        	println(inner(psi,N1,psi))
        end
        
        for i = 1:N
        	numb2 = AutoMPO()
        	numb2 += "C2d",2*i, "C2", 2*i
        	N2 = MPO(numb2,sites)
        	println(inner(psi,N2,psi))
        end
        return
end          

I see, thanks. At a high level, there’s no reason to expect that DMRG will respect the seperate particle numbers unless there something that strictly prevents this, such as separately conserving the quantum numbers. Otherwise, you can think of DMRG as producing “virtual” processes which aren’t in the Hamiltonian but just come from small, random perturbations of the wavefunction which “nucleate” a small density of particles on the other sites. Then if it’s beneficial to put more particles on those sites to lower the energy, then DMRG will rapidly do that.

I think the only solution to this issue besides enforcing the separate particle number conservation (which you say you can’t do for your case) would be to put an artificial energy penalty, such as an on-site potential, to prevent DMRG from putting particles on the sites you don’t want.

Thank you Miles. I suspected the same thing. I feel that once the flavor-mixing interaction is added, no matter how small, the final result will be correct, i.e. this is only an issue for the free case without explicit mixing terms. That said, I’m current trying to perturbatively calculate the small interaction case to see if my intuition is correct or not. I’ll update you once I figure this out.

Hi Miles, sorry for the late response. Instead of doing perturbations, I directly did ED, and the results agree in the presence of flavor-mixing interactions. So we can conclude that the above problem only exists when there are no flavor-mixing interactions.

Thanks for the update. So it sounds you were able to confirm the code is correct now?

Yes. Thanks again.