Question on TEBD for fermionic systems

ITensor Julia Questions
,
1

Let’s consider a quantum quench scenario. At t=0, system is in the ground state of some Hamiltonian H(\lambda_0) with parameter \lambda_0. And then it evolves under another parameter H(\lambda), |\psi(t)\rangle=e^{-\mathrm{i} tH(\lambda)}|\psi_0\rangle. This problem can be studied using DMRG and TEBD in ITensors.jl.

When \lambda = \lambda_0, system should stay in the initial state up to a phase as time evolves, |\psi(t)\rangle=e^{\mathrm{i}\theta(t)}|\psi_0\rangle, and all the physical observables should be unchanged as time evolves. I use this fact to do a basic benchmark of the TEBD results. Indeed, many spin models past this test, including transversed filed Ising model (with OBC or PBC) and even more complicated models with next-nearest-neighboring and three-body interaction (e.g., S_i^xS_{i+1}^xS_{i+2}^z). For small L, the TEBD results also converge to ED results as the Trotter decomposition time slice getting smaller \Delta_t\rightarrow0.

However, when I start to play with fermionic system, this simply test can not be pasted. For free fermions with only NN hopping t and chemical potential \mu, this test is trivial, since the ground state of \mu/t is also an eigenstate of arbitrary \mu^\prime/t. And the state is exactly unchanged (in general cases, the results are only almost unchanged due to Trotter error). So we consider the Kitaev’s chain model with OBC,

\begin{equation} \hat{H}=-\sum_{j=1}^{L-1}\left[tc_{j}^{\dagger}c_{j+1}+tc_{j+1}^{\dagger}c_{j}+\Delta c_{j}c_{j+1}+\Delta c_{j+1}^{\dagger}c_{j}^{\dagger}\right]-\mu\sum_{j=1}^L\left(c_{j}^{\dagger}c_{j}-\frac{1}{2}\right).\tag{1} \end{equation}

I choose t=\Delta=1.0,\mu=0.5 (which is not a gapless point). However, the results are not expected. The average density n(t)=\langle \psi(t)|\sum_j c^\dagger_j c_j|\psi(t)\rangle/L is not unchanged, as shown in the figure below. And the code can be found below.

image
image1356×998 98.9 KB

using ITensors, ITensorMPS

function Kitaev_chain_MPO(L::Int64, t::Float64, μ::Float64, Δ::Float64)
    sites = siteinds("Fermion",L)
    os = OpSum()
    for j ∈ 1:L
        os += -μ, "n", j
    end
    for j ∈ 1:L-1
        # hopping
        os += -t, "c†", j, "c", j+1
        os += -t, "c†", j+1, "c", j
        # pairing
        os += -Δ, "c", j, "c", j+1
        os += -Δ, "c†", j+1, "c†", j
    end
    H = MPO(os,sites)
    # initialize the MPS
    psi0 = random_mps(sites;linkdims=10)
    return sites, H, psi0
end
function ITensors.op(::OpName"expitKC", ::SiteType"Fermion", s1::Index, s2::Index; Δt,t,μ,Δ,bc_c=1.0)
    h = (-μ) * op("n", s1) * op("I", s2) # chemical potential
    # hopping
    # h +=   (-t) * product(op("c†", s1),op("F",s1)) * op("c", s2)
    #      - (-t) * product(op("c", s1),op("F",s1)) * op("c†", s2)
    h +=   (-t) * op("c†", s1) * op("c", s2)
         + (-t) * op("c†", s2) * op("c", s1)
    # pairing
    # h +=   (-Δ) * product(op("c", s1),op("F",s1)) * op("c", s2)
    #      - (-Δ) * product(op("c†", s1),op("F",s1)) * op("c†", s2)
    h += - (-Δ) * op("c", s1) * op("c", s2)
         - (-Δ) * op("c†", s2) * op("c†", s1)
    return exp(-im * Δt / 2 * h)
end
function average_density(L::Int64, psi::MPS)
    ns = expect(psi,"n")
    ave_n = sum(ns)/L
    return ave_n, ns
end


## evolution using TEBD: 
## use the same initial Hamiltonian and quench Hamiltonian, 
## so the state is unchanged. 
begin # OBC
    L = 8
    t = 1.0
    μ = 0.5
    Δ = 1.0
    # initial ground state
    sites, H, psi0 = Kitaev_chain_MPO(L, t, μ, Δ) # OBC
    nsweeps = 25
    maxdims = [2,2,8,8,20,20,50,50,100,100,200,200,400]
    cutoff1 = [1E-11]
    noise1 = [1E-4, 1E-4, 1E-5, 1E-6, 1E-7, 1E-8, 0.0]
    Eg, psi = dmrg(H,psi0;nsweeps,maxdim=maxdims,cutoff=cutoff1, noise=noise1)
    n0 = average_density(L,psi)[1]
    @show n0
end
begin
    # quench parameters of Hamiltonian
    t_q = t
    μ_q = μ
    Δ_q = Δ
    # TEBD
    cutoff2 = 1E-10
    Δt = 0.01
    Ttot = 2.0
    # observables
    nt = zeros(ComplexF64, Int(Ttot/Δt)+1)
    # Make gates (1,2),(2,3),(3,4),...
    gates = ITensor[]
    for j in 1:L-1
        s1 = sites[j]
        s2 = sites[j+1]
        Gj = op("expitKC", s1, s2; Δt=Δt, t=t_q, μ=μ_q, Δ=Δ_q)
        push!(gates, Gj)
    end
    # Include gates in reverse order too
    # (N-1,N),(N-2,N-1),(N-3,N-2)...
    append!(gates, reverse(gates))
    # Compute <n> at each time step
    # then apply the gates to go to the next time
    i = 1
    for t in 0.0:Δt:Ttot
        nt[i] = average_density(L,psi)[1]
        t≈Ttot && break

        psi = apply(gates, psi; cutoff=cutoff2)
        normalize!(psi)
        i += 1
    end
    @show nt[1:10]
end

Did I make some mistakes in the understanding of TEBD for fermionic systems?

The following are my own understanding about how to deal with the Jordan-Wigner strings and I suspect some of them are incorrect.

  1. DMRG: as @miles pointed out in https://itensor.discourse.group/t/proper-way-to-construct-fermion-gates/172/2?u=wangfh5, the Opsum system would deal with the JW strings properly. So I just write the Hamiltonian according to Eq.(1). And the results (ground-state energy and \langle n(0)\rangle) are consistent with ED.
  2. TEBD: there are two steps in TEBD,
    1. Construct local gates: according to some related posts and practical experience, the local fermionic gate should consider the JW strings by hand. Following the doc in https://itensor.org/docs.cgi?page=tutorials/fermions, I figure out two equivalent approaches for NN gates (as shown in the function ITensors.op(::OpName"expitKC", ... in the code) and they are formulated below.
    \begin{aligned} c_{i}^{\dagger}c_{i+1}&=\hat{S}_{i}^{+}F_{i}\hat{S}_{i+1}^{-}=\texttt{product(op("cdag",i),op("F",i))*op("c",i+1)}\\c_{i+1}^{\dagger}c_{i}=-c_{i}c_{i+1}^{\dagger}&=-\hat{S}_{i}^{-}F_{i}\hat{S}_{i+1}^{+}=-\texttt{product(op("c",i),op("F",i))*op("cdag",i+1)}\\c_{i}c_{i+1}&=\hat{S}_{i}^{-}F_{i}\hat{S}_{i+1}^{-}=\texttt{product(op("c",i),op("F",i))*op("c",i+1)}\\c_{i+1}^{\dagger}c_{i}^{\dagger}=-c_{i}^{\dagger}c_{i+1}^{\dagger}&=-\hat{S}_{i}^{+}F_{i}\hat{S}_{i+1}^{+}=-\texttt{product(op("cdag",i),op("F",i))*op("cdag",i+1)} \end{aligned}
    \begin{aligned} c_{i}^{\dagger}c_{i+1}&=\hat{S}_{i}^{+}\hat{S}_{i+1}^{-}=\texttt{op("c†",i)*op("c",i+1)}\\c_{i+1}^{\dagger}c_{i}=-c_{i}c_{i+1}^{\dagger}&=\hat{S}_{i}^{-}\hat{S}_{i+1}^{+}=\texttt{op("c",i)*op("c†",i+1)}=\texttt{op("c†",i+1)*op("c",i)}\\c_{i}c_{i+1}&=-\hat{S}_{i}^{-}\hat{S}_{i+1}^{-}=-\texttt{op("c",i)*op("c",i+1)}\\c_{i+1}^{\dagger}c_{i}^{\dagger}=-c_{i}^{\dagger}c_{i+1}^{\dagger}&=-\hat{S}_{i}^{+}\hat{S}_{i+1}^{+}=-\texttt{op("c†",i)*op("c†",i+1)}=-\texttt{op("c†",i+1)*op("c†",i)} \end{aligned}
    Both of them give the same results.
    2. Apply the local gates to the MPS: again, as @miles pointed out in https://itensor.discourse.group/t/proper-way-to-construct-fermion-gates/172/2?u=wangfh5, the apply system should deal with the JW strings properly. I am not sure how it work.

Any idea and advice is appreciated!
Fo-Hong

2

I found that TDVP give me excellent results in all cases today. So in practical, the question in this post is not essential to me, but I am still curious about the right way to perform TEBD of fermionic systems.

1 Like