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

ITensor Julia Questions
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