# Question on TEBD for fermionic systems **URL:** \ **Category:** ITensor Julia Questions\ **Tags:** julia, fermions\ **Created:** [September 16, 2024, 2:30pm UTC](https://itensor.discourse.group/t/question-on-tebd-for-fermionic-systems/1991 "2024-09-16T14:30:49Z")\ **Posts on this page:** 2\ **Page:** 1
**Author:** ![wangfh5](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/wangfh5/32/619_2.png) [@wangfh5](https://itensor.discourse.group/u/wangfh5)\ **Post date:** [September 16, 2024, 2:30pm UTC](https://itensor.discourse.group/t/question-on-tebd-for-fermionic-systems/1991/1 "2024-09-16T14:30:49Z")
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](https://itensor.github.io/ITensors.jl/dev/tutorials/DMRG.html) and [TEBD](https://itensor.github.io/ITensors.jl/dev/tutorials/MPSTimeEvolution.html) 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](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/94fe9eb93f5734755eca758544ecb5ece775023e.jpeg) ```julia 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 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](https://itensor.discourse.group/t/proper-way-to-construct-fermion-gates/172/2), 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](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](https://itensor.discourse.group/t/proper-way-to-construct-fermion-gates/172/2), 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 ---
**Author:** ![wangfh5](https://yyz2.discourse-cdn.com/free1/user_avatar/itensor.discourse.group/wangfh5/32/619_2.png) [@wangfh5](https://itensor.discourse.group/u/wangfh5)\ **Post date:** [September 17, 2024, 6:11am UTC](https://itensor.discourse.group/t/question-on-tebd-for-fermionic-systems/1991/2 "2024-09-17T06:11:43Z")
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.