Incorrect fermionic anticommutation relationships

[acroscarrillo]

I'm fairly new to the package but it seems to me like the fermionic anticommutation relations

$$\{a_i, a^\dagger_j\}\doteq a_i a^\dagger_j + a^\dagger_j a_i = \delta_{ij} \quad \text{and} \quad \{a_i^\dagger, a^\dagger_j\} = \{a_i, a_j\} = 0,$$

have been incorrectly implemented:

julia> b = NLevelBasis(2)
NLevel(N=2)

julia> b_mb = ManyBodyBasis(b, fermionstates(b, [0,1,2]))
ManyBody(onebodybasis=NLevel(N=2), states:4)

julia> (create(b_mb,2)*create(b_mb,1)+create(b_mb,1)*create(b_mb,2))*basisstate(b_mb, [0,0])
Ket(dim=4)
  basis: ManyBody(onebodybasis=NLevel(N=2), states:4)
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 2.0 + 0.0im

Since, if the code does what I think it should be doing, the last operation should be

$$(a^\dagger_2 a^\dagger_1 + a^\dagger_1 a^\dagger_2 )|0,0\rangle = ( a^\dagger_2 a^\dagger_1 - a^\dagger_2 a^\dagger_1 )|0,0\rangle = 0$$

but instead it is (wrongly) doing

$$(a^\dagger_2 a^\dagger_1 + a^\dagger_1 a^\dagger_2 )|0,0\rangle = 2* a^\dagger_2 a^\dagger_1 |0,0\rangle = 2*|1,1\rangle$$

as if they were bosons. Interestingly:

julia> (destroy(b_mb,1)*create(b_mb,1)+create(b_mb,1)*destroy(b_mb,1))*basisstate(b_mb, [0,0])
Ket(dim=4)
  basis: ManyBody(onebodybasis=NLevel(N=2), states:4)
 1.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im
 0.0 + 0.0im

which is correct but looks like a fluke judging the above issue.

[amilsted]

I think both results are correct for bosons. I think the only thing fermionstates does is implement "pauli exclusion", so disallowing >1 fermion in each mode. This is understandably misleading. Perhaps @david-pl can correct me?

[david-pl]

I admit to being a bit rusty myself when it comes to ManyBodyBasis. But yes, @amilsted is correct, fermionstates simply limits the number of particles per site (per level of NLevelBasis(2) in the above example) to 1.

With create(b_mb, 1) * create(b_mb, 2) you're creating one fermion that is in state 1 and another one that is in state 2. Since the particle number is not conserved, this is fine and the order of the operations doesn't matter. So the anti-commutation relation isn't applicable here.

[acroscarrillo]

Hey thank you both for chipping in!

With create(b_mb, 1) * create(b_mb, 2) you're creating one fermion that is in state 1 and another one that is in state 2.

The above statement is true, however

Since the particle number is not conserved, this is fine and the order of the operations doesn't matter. So the anti-commutation relation isn't applicable here.

this is not true: fermionic operator anticommute at different sites regardless of any symmetry of the state. That is $c_1 c_2 = - c_2 c_1$ is true at the level of the underlying algebra. An immediate consequence of this relation is that $c_1 c_1 = 0$ and the same goes for it's complex conjugate: this is what limits the occupation number to be between 0 and 1 (no need to hard code it). These fermionic operators can be implemented via a Jordan Wigner transformation of the Pauli matrices.

Im not an expert on the package, but I believe this behaviour needs to be implemented to correctly characterise fermions (that is if it isn't already implemented and I am missing something!)

[david-pl]

@acroscarrillo no this asymmetry relation for fermions is not implemented.

[acroscarrillo]

I see, I think that so long as a heads up is given in the documentation it should be okay.
For what is worth, I wrote some code to sketch how one could implement fermionic operators via the jordan wigner transformation. Unfortunately, Im not familiar enough with the library to write anything more in line with it but I hope someone with more experience could easily translate it if so wished. Thank you all for your time and dedication :)

using SparseArrays

const sigma_x = [[0.0 1.0]; [1.0 0.0]]
const sigma_y = [[0.0 -1.0 * im]; [1.0 * im 0.0]]
const sigma_z = [[1.0 0.0]; [0.0 -1.0]]
const sigma_p = 0.5 * (sigma_x + 1 * im * sigma_y)
const sigma_m = sigma_p'
const sigma_0 = [[1.0 0.0]; [0.0 1.0]]

⨷(a, b) = kron(a, b)

function pauli_string(L, j)
    if j > L
        throw(ArgumentError("Requested pauli string at site j with j>L."))
    end

    if j > 1
        string = sparse(-sigma_z)
        for _ in 2:(j - 1)
            string = string ⨷ (-sigma_z)
        end
        for _ in j:L
            string = string ⨷ sigma_0
        end
        return string
    else
        return sparse(I(2^L))
    end
end

function create(L, j)
    I_left = sparse(I(2^(j - 1)))
    I_right = sparse(I(2^(L - j)))
    return pauli_string(L, j) * (I_left ⨷ sparse(sigma_p) ⨷ I_right)
end

function destroy(L, j)
    return create(L, j)'
end

[aryavorskiy]

Okay, I may have some ideas on how to implement the correct anticommutation relations, but this might be a breaking change, so RFC:

Let us define the filling number states as follows:

$$| N_1^{i_1} N_2^{i_2} \ldots \rangle = C \cdot (\hat{a}^\dagger_{i_1})^{N_1} (\hat{a}^\dagger_{i_2})^{N_2} \ldots |0\rangle, \\ \\ \\\ i_1 < i_2 < \ldots$$

The ordering of the creation operators here is unimportant for bosons, but extremely important for fermions. Also for fermions we must have N_j = 1, but this is already taken care of by the "Pauli exclusion" that is already implemented.

Thus, if we want to find the result of the creation operator, we can do the following:

$$\hat{a}^\dagger_j | 1^{i_1} 1^{i_2} \ldots \rangle = \hat{a}^\dagger_j \hat{a}^\dagger_{i_1} \hat{a}^\dagger_{i_2} \ldots |0\rangle = (-1)^k \hat{a}^\dagger_{i_1} \ldots \hat{a}^\dagger_{i_k} \hat{a}^\dagger_j \hat{a}^\dagger_{i_{k+1}} \ldots |0\rangle, \\ \\ \\\ i_k < j < i_{k+1}$$

The result of the annihilation operator can be carried out in a similar way.

The simplest way of somehow setting the particle statistics is adding a isbosestatistics::Bool field to the ManyBodyBasis struct. However, we can implement Abelian anyons by remembering the exp(im * theta) phase instead:

$$\hat{a}_i \hat{a}_j - e^{i \theta} \hat{a}_j \hat{a}_i = 0, \\ \\ \\\ \hat{a}_i \hat{a}^\dagger_j - e^{i \theta} \hat{a}^\dagger_j \hat{a}_i = \delta_{i j}$$

But do we really need that complication?

Will file a PR ASAP if no objections.

[Krastanov]

Thank you @acroscarrillo for reporting this and thank you @aryavorskiy for the fix. It is now merged (in #160 ) and a release will be made shortly.