Manual construction of nearest neighbor hopping of spinless Fermions on a 1D chain

ITensor Julia Questions
1

Hello everyone,

I am currently trying to construct a hopping term between nearest neighbors for spinless Fermions on a one dimensional chain with L sites manually, i.e. without using OpSum (I want to learn the manual construction since my actual Hamiltonian contains exponentials of operators, which, to my understanding, do not work in OpSum). This is my Hamiltonian:

H = -t \sum_{j=1}^{L-1} \left( c^\dagger_j c_{j+1} + c^\dagger_{j+1} c_j \right)

I follow the first half of the paper by Crosswhite and Bacon (https://arxiv.org/pdf/0708.1221) for constructing MPOs by using finite automata. To apply this formalism, I need the terms to be site ordered (at least to my understanding), which is why I rewrite the Hamiltonian

H = -t \sum_{j=1}^{L-1} \left( c^\dagger_j c_{j+1} - c_j c^\dagger_{j+1} \right)

and then construct the boundary matrices W_0 and W_L as well as the bulk matrices W_bulk by writing down the finite state automata and the corresponding matrix product diagram. Below is the code for three functions, one for constructing the MPO with OpSum, one for constructing it manually and one for constructing it manually and neglecting the sign change:

using ITensors
using ITensorMPS
using Random

function hopping_opsum(sites::Vector{<:Index}, t::Real=1.0)::MPO
    os = OpSum()
    for j in 1:(length(sites)-1)
        os += -t, "c†", j, "c", j + 1
        os += -t, "c†", j + 1, "c", j # or +t, "c", j, "c†", j + 1 due to anticommutation
    end
    return MPO(os, sites)
end

function hopping_manual(sites::Vector{<:Index}, t::Real=1.0)::MPO
    # Local operators
    id = ComplexF64[1 0; 0 1]
    n = ComplexF64[0 0; 0 1]
    c_dag = ComplexF64[0 0; 1 0]
    c = ComplexF64[0 1; 0 0]

    # Virtual bond indices (dimension 4), excluding boundaries
    L = length(sites)
    links = [Index(4, "link,l=$i") for i in 1:(L - 1)]
    W = Vector{ITensor}(undef, L) # undef: do not initialize yet

    # First site
    W[1] = ITensor(links[1], prime(sites[1]), sites[1])
    for i in 1:2, j in 1:2
        W[1][1, i, j] = id[i, j]
        W[1][2, i, j] = -t * c_dag[i, j]
        W[1][3, i, j] = +t * c[i, j] # !! sign change due to anticommutation
    end

    # Bulk sites
    for l in 2:(L-1) # n is already used
        W[l] = ITensor(links[l-1], links[l], prime(sites[l]), sites[l])
        for i in 1:2, j in 1:2
            W[l][1, 1, i, j] = id[i, j]
            W[l][1, 2, i, j] = -t * c_dag[i, j]
            W[l][1, 3, i, j] = +t * c[i, j] # !! sign change due to anticommutation

            W[l][2, 4, i, j] = c[i, j]
            W[l][3, 4, i, j] = c_dag[i, j]
            W[l][4, 4, i, j] = id[i, j]
        end
    end

    # last site
    W[L] = ITensor(links[L-1], prime(sites[L]), sites[L])
    for i in 1:2, j in 1:2
        W[L][2, i, j] = c[i, j]
        W[L][3, i, j] = c_dag[i, j]
        W[L][4, i, j] = id[i, j]
    end
    return MPO(W)
end

function hopping_manual_changed_sign(sites::Vector{<:Index}, t::Real=1.0)::MPO
    # Local operators
    id = ComplexF64[1 0; 0 1]
    n = ComplexF64[0 0; 0 1]
    c_dag = ComplexF64[0 0; 1 0]
    c = ComplexF64[0 1; 0 0]

    # Virtual bond indices (dimension 4), excluding boundaries
    L = length(sites)
    links = [Index(4, "link,l=$i") for i in 1:(L - 1)]
    W = Vector{ITensor}(undef, L) # undef: do not initialize yet

    # First site
    W[1] = ITensor(links[1], prime(sites[1]), sites[1])
    for i in 1:2, j in 1:2
        W[1][1, i, j] = id[i, j]
        W[1][2, i, j] = -t * c_dag[i, j]
        W[1][3, i, j] = -t * c[i, j] # !!
    end

    # Bulk sites
    for l in 2:(L-1) # n is already used
        W[l] = ITensor(links[l-1], links[l], prime(sites[l]), sites[l])
        for i in 1:2, j in 1:2
            W[l][1, 1, i, j] = id[i, j]
            W[l][1, 2, i, j] = -t * c_dag[i, j]
            W[l][1, 3, i, j] = -t * c[i, j] # !!

            W[l][2, 4, i, j] = c[i, j]
            W[l][3, 4, i, j] = c_dag[i, j]
            W[l][4, 4, i, j] = id[i, j]
        end
    end

    # last site
    W[L] = ITensor(links[L-1], prime(sites[L]), sites[L])
    for i in 1:2, j in 1:2
        W[L][2, i, j] = c[i, j]
        W[L][3, i, j] = c_dag[i, j]
        W[L][4, i, j] = id[i, j]
    end
    return MPO(W)
end

let
    # setup for simple system
    Random.seed!(42)
    L = 3
    sites = siteinds("Fermion", L)
    H_opsum = hopping_opsum(sites)
    H_manual = hopping_manual(sites)
    H_manual_changed = hopping_manual_changed_sign(sites)

    # compute energies for a random s
    s = randomMPS(sites, 2)
    @show energy_opsum = inner(s', H_opsum, s)
    @show energy_manual = inner(s', H_manual, s)
    @show energy_manual_changed = inner(s', H_manual_changed, s)

    # compute action of Hamiltonians on s
    s_opsum = apply(H_opsum, s)
    s_manual = apply(H_manual, s)
    s_manual_changed = apply(H_manual_changed, s)

    # compare results with opsum as reference
    @show norm(s_opsum - s_manual)
    @show norm(s_opsum - s_manual_changed)

end

nothing

The lines affected by the sign change are marked with # !!. The let block is for quickly checking, if the MPOs encode the same operator, i.e. if they act the same on a random MPS. Executing this code within the REPL yields:
```
julia> include(“hopping.jl”)
energy_opsum = inner(s’, H_opsum, s) = -0.07088577164847563
energy_manual = inner(s’, H_manual, s) = -1.629062980549775e-17 + 0.0im
energy_manual_changed = inner(s’, H_manual_changed, s) = -0.07088577164847563 + 0.
0im
norm(s_opsum - s_manual) = 1.085567382692007
norm(s_opsum - s_manual_changed) = 1.7775626050888642e-16
```
which brings me to the following questions:

  1. Is the general syntax that I use to manually construct the MPO correct? I have googled a lot during the last weeks on how to achieve this in ITensors, but have found no nice reference/example.

  2. Why does the function with the ‘wrong’ sign convention work correctly, i.e. match the result from OpSum? It feels like I have not yet fully understood how to construct any Hamiltonian including only nearest neighbor interactions with the finite automata formalism.

  3. Does the formalism change extensively when implementing periodic boundary conditions and Abelian symmetries?

I am very new to Julia and ITensor, so please forgive all the silly mistakes I probably did. I look forward to hearing from you!

2
  1. Seems fine to me. What is the Hamiltonian you are trying to construct? Why do you expect the OpSum stuff not to work? I’m familiar with the algorithm OpSum uses but much less so with the finite automata approach, so I’m interested to know why that is a better choice.
  2. This is because you are not taking the Jordan Wigner string into account. Consider following four site example
c^\dagger_2 \rightarrow Z \otimes \sigma_- \otimes I \otimes I \\ c_3 \rightarrow Z \otimes Z \otimes \sigma_+ \otimes I \ .

Then

c^\dagger_2 c_3 \rightarrow I \otimes (\sigma_- Z) \otimes \sigma_+ \otimes I \\ = I \otimes \sigma_- \otimes \sigma_+ \otimes I \ .

One could use c^\dagger_3 c_2 = (c^\dagger_2 c_3)^\dagger to arrive at c^\dagger_3 c_2 \rightarrow I \otimes \sigma_+ \otimes \sigma_- \otimes I. But to be explicit

c^\dagger_3 c_2 \rightarrow (Z \otimes Z \otimes \sigma_- \otimes I) (Z \otimes \sigma_+ \otimes I \otimes I) \\ = I \otimes (Z \sigma_+) \otimes \sigma_- \otimes I \\ = I \otimes \sigma_+ \otimes \sigma_- \otimes I \ .

You could also do the same thing by first transforming c^\dagger_3 c_2 = - c_2 c^\dagger_3 but then you wind up with the term \sigma_+ Z = -\sigma_+ cancelling out the minus sign introduced by the anticommutation.

  1. This would work with almost no technical changes for Abelian symmetries. As for periodic boundary conditions you’ll need to add in new bond dimensions to carry the interaction between the ends. This will require some technical changes.
2 Likes
3

First of all thank you for your quick and detailed reply!

  1. I later want to include a (constant) photon mode which dresses the hopping terms via the quantized Peierls phase \exp(\pm \frac{ig}{\sqrt{L}} \left( a^\dagger + a \right)) . After reading your reply, I realized that you are right and that this is possible via OpSum. (If I understand correctly: Determine the matrix representation of the operator on a truncated Hilbert space, setting some appropriate cutoff for the photon number. Use the ITensor.op function to define the operator such that one can use it in OpSum.)

    What algorithm does OpSum actually use? I didn’t manage to find a reference within the documentation and I was quite overwhelmed by the source code. Does ITensors OpSum work something like that:

    i. Construct MPOs for all of the single site operators c^\dagger_i , \, c_i for i = 1,\dots, L
    ii. Take the desired products of single site MPOs to obtain the hopping terms c^\dagger_{i+1} c_i , \, c^\dagger_{i} c_{i+1}
    iii. Sum up the resulting MPOs to obtain the hopping Hamiltonian
    iv. Compress the MPO via SVD or some other algorithm

    as outlined in this paper (https://arxiv.org/pdf/1611.02498)?

    The reason why I started to learn the FSM approach was, that I saw it in a short review by Schollwöck (https://royalsocietypublishing.org/rsta/article/369/1946/2643/114494/The-density-matrix-renormalization-group-a-short) and I also wanted to be able to do some pen and paper derivation to test my code (and my understanding). The FSM approach seemed especially nice for the pen and paper derivation, since it doesn’t really depend on the number of sites and many actual calculations as summing two MPOs etc.

  2. What I don’t understand about the Jordan-Wigner transformation (JW) for my case is, how it actually affects the manual construction of the Hamiltonian. I understand, that I could use it to map my system to one of spins. But I want the resulting Hamiltonian to be defined on a chain of Fermions and not on one of spins. I can’t quite figure out, what the JW means for my fermionic system and the corresponding operator constructed via FSM.

  3. Okay, that is good to hear thank you!

4

1

Yes, that approach should work great. The OpSum algorithm will work best (it will produce an optimal bond dimension in many cases) if the single site operators you use to expand the Hamiltonian are linearly independent. For example if on site j you have terms involving I_j, Z_j and \exp(i t Z_j) you might achieve a smaller bond dimension by rewriting the OpSum to eliminate one of those terms.

ITensors uses the algorithm described in section IV A of https://doi.org/10.1063/1.4955108. The core of the algorithm takes a bipartite operator

\hat{O} = \sum_{a = 1}^{N_A} \sum_{b = 1}^{N_B} \gamma_{ab} \hat{A}_a \otimes \hat{B}_b \ ,

and turns it into a two-site MPO of bond dimension w

\hat{O} = \sum_{\chi = 1}^w \left( \sum_{a = 1}^{N_A} \alpha_{\chi a} \hat{A}_a \right) \otimes \left( \sum_{b = 1}^{N_B} \beta_{\chi b} \hat{B}_b \right) \ ,

where \alpha = U and \beta = S V^\dagger obtained from the SVD of \gamma, and w is the rank.

2

At the end of the day the fermion sites are simply qubit sites with some fancy magic behind the scenes to account for the anti-commutivity of the operators (done by the Jordan-Wigner mapping). Most of the time as a user you don’t notice the difference from a “true fermion” site (whatever that is), but if you are implementing something, say MPO construction, you need to add the magic yourself. This can lead to misleading translations from pen-and-paper MPOs to code.

For example consider

H = c^\dagger_1 c_2 + c^\dagger_2 c_1 \ ,

which could be written as a product of operator valued matrices as

H = \begin{pmatrix} c & c^\dagger \end{pmatrix} \begin{pmatrix} -c^\dagger \\ c \end{pmatrix} = -c_1 c^\dagger_2 + c^\dagger_1 c_2 \ .

This is almost an MPO, but since the fermionic operators are non-local it is not technically an MPO. Your mistake is that you used this non-local construction but then used local definitions of the fermionic operators (see c_dag = ComplexF64[0 0; 1 0]).

5

1. Okay, thank you for the reference and the brief summary!

2.

Yes, now I understand. My manual code works without the added sign change due to the Jordan-Wigner mapping. As you implied, the site type is not important for manual constructions. I constructed a minimal example to show the problem and how to solve it. One can still manually construct operators via finite state machines, which I find pedagogically useful.

Lets examine \hat{H} = \sum_{j=1}^L \hat{c}^\dagger_j. The fermionic creation operator can also be written as \hat{c}^\dagger_j = \hat{Z}_1 \dots \hat{Z}_{j-1} \hat{a}^\dagger_j, with the bosonic creation operator \hat{a}^\dagger_j and the Jordan-Wigner strings \hat{Z}_j = 1 - 2\hat{n}_j. One can use this to rewrite the Hamiltonian as \hat{H} = \sum_{j=0}^L \hat{Z}_1 \dots \hat{Z}_{j-1} \hat{a}^\dagger_j with the benefit that the anticommutation relation is embedded directly in the string operators.

For L = 2, we have for example H |10\rangle = - |11\rangle, using the convention |11\rangle \equiv c^\dagger_2 c^\dagger_1 |00 \rangle (= - c^\dagger_1 c^\dagger_2 |00 \rangle). The code below shows the “naive” approach, which does not implement the anticommutation relation correctly and the Jordan-Wigner approach:

using ITensorMPS
using ITensors

function fermionic_creator_sum_wrong_manual(sites::Vector{<:Index})::MPO
    _check_spinless_fermions_sites(sites)

    # local operators
    c_dag = ComplexF64[0 0; 1 0]
    id = ComplexF64[1 0; 0 1]

    # Virtual bond indices (dimension 2), excluding boundaries
    L = length(sites)
    links = [Index(2, "link,l=$i") for i in 1:(L - 1)]

    # allocate vector of length L to hold site tensors
    W = Vector{ITensor}(undef, L) # undef: do not initialize yet

    # First site
    W[1] = ITensor(links[1], prime(sites[1]), sites[1])
    for i in 1:2, j in 1:2
        W[1][1, i, j] = c_dag[i, j]
        W[1][2, i, j] = id[i, j]
    end

    # Bulk sites
    for n in 2:(L - 1)
        W[n] = ITensor(links[n - 1], links[n], prime(sites[n]), sites[n])
        for i in 1:2, j in 1:2
            W[n][1, 1, i, j] = id[i, j]
            W[n][2, 1, i, j] = c_dag[i, j]
            W[n][2, 2, i, j] = id[i, j]
        end
    end

    # Last site
    W[L] = ITensor(links[L - 1], prime(sites[L]), sites[L])
    for i in 1:2, j in 1:2
        W[L][1, i, j] = id[i, j]
        W[L][2, i, j] = c_dag[i, j]
    end

    return MPO(W)
end

function fermionic_creator_sum_manual(sites::Vector{<:Index})::MPO
    _check_spinless_fermions_sites(sites)

    # local operators
    a_dag = ComplexF64[0 0; 1 0]
    z = ComplexF64[1 0; 0 -1] # Jordan-Wigner string operator
    id = ComplexF64[1 0; 0 1]

    # Virtual bond indices (dimension 2), excluding boundaries
    L = length(sites)
    links = [Index(2, "link,l=$i") for i in 1:(L - 1)]

    # allocate vector of length L to hold site tensors
    W = Vector{ITensor}(undef, L) # undef: do not initialize yet

    # First site
    W[1] = ITensor(links[1], prime(sites[1]), sites[1])
    for i in 1:2, j in 1:2
        W[1][1, i, j] = a_dag[i, j]
        W[1][2, i, j] = z[i, j]
    end

    # Bulk sites
    for n in 2:(L - 1)
        W[n] = ITensor(links[n - 1], links[n], prime(sites[n]), sites[n])
        for i in 1:2, j in 1:2
            W[n][1, 1, i, j] = z[i, j]
            W[n][2, 1, i, j] = a_dag[i, j]
            W[n][2, 2, i, j] = id[i, j]
        end
    end

    # Last site
    W[L] = ITensor(links[L - 1], prime(sites[L]), sites[L])
    for i in 1:2, j in 1:2
        W[L][1, i, j] = id[i, j]
        W[L][2, i, j] = a_dag[i, j]
    end

    return MPO(W)
end

let
    L = 2
    sites = siteinds("Fermion", L)
    H_wrong = fermionic_creator_sum_wrong_manual(sites)
    H = fermionic_creator_sum_manual(sites)

    initial_state = ["1", "0"]

    final_state = ["1", "1"]

    psi_initial = MPS(sites, initial_state)
    psi_final = MPS(sites, final_state)

    @show inner(psi_final', H_wrong, psi_initial) # should be -1 but is 1
    @show inner(psi_final', H, psi_initial) # should be -1
end

nothing

The difference between the two approaches, becomes clear when computing for example \langle 11 | H | 01 \rangle = \langle 11 | (-1 | 11 \rangle ) = -1. The above code evaluates to

```
julia> include(“example.jl”)
inner(psi_final’, H_wrong, psi_initial) = 1.0 + 0.0im
inner(psi_final’, H, psi_initial) = -1.0 + 0.0im
```

which shows, that the Jordan-Wigner mapping takes care of the anticommutations. This example is I think even an example for something that OpSum cannot do yet because “Parity-odd fermionic terms [are] not yet supported by OpSum to MPO conversion”.

Thank you @corbett5 for all your help!

1 Like
6

This topic was automatically closed 10 days after the last reply. New replies are no longer allowed.