How to apply a single creation/annihilation operator to a finite MPS ?

[Jason-zo]

Dear Maartenvd,

Thanks for you and other contributors' great work!
I have a problem that how to apply a single creation/annihilation operator to a finite MPS and get a new MPS?

Usually, we need not to consider this issue because the creation operator and the annihilation operator's auxiliary legs A are connected (TensorMap(..., V ← V $\otimes$ A) and TensorMap(..., A $\otimes$ V ← V) ) that constructs a single site or two sites operator like $\hat{O_i}=c^\dagger_i c_i$ or $\hat{O_{ij}} = c^\dagger_i c_j$. In this case, it is easy to calculate expectation or correlation value $\langle 0|\hat{O}|0\rangle$.

However, when it comes to the time evolution, such as $\langle 0|S^+(t) S^-(0) |0\rangle$ or $\langle 0|c(t) c^\dagger(0) |0\rangle$, the above question makes sense that we need to evolve the MPS $S^- |0\rangle$ or $c^\dagger |0\rangle$ into $e^{-iHt} S^- |0\rangle$ or $e^{-iHt}c^\dagger|0\rangle$.

To obtain the new MPS ($S^- |0\rangle$ or $c^\dagger |0\rangle$), one need to deal with the single auxiliary vector space A of the creation/annihilation operator. Please forgive me as a beginner, my naive idea is to add identity operators on the other sites, where the left and right virtual spaces of the identity operators are the auxiliary vector space A. I construct

$c^\dagger_1 \mathbb{1}_2$
$c_1 \mathbb{1}_2$

as

    ID1 = TensorMap(ones, ComplexF64, As ⊗ Ps ← Ps)
    ID2 = TensorMap(ones, ComplexF64, Ps ← As ⊗ Ps)
    if ifdagger
        @planar twosite[-1 -2; -3 -4] := c⁺[-1; -3 1] * ID1[1 -2; -4]
    else
        @planar twosite[-1 -2; -3 -4] := c[-1 1; -3] * ID2[-2; 1 -4]
    end
    return FiniteMPO(twosite)

and

$c^\dagger_1 \mathbb{1}_2 \mathbb{1}_3 \mathbb{1}_4$
$c_1 \mathbb{1}_2 \mathbb{1}_3 \mathbb{1}_4$

as

    ID1 = TensorMap(ones, ComplexF64, As ⊗ Ps ← Ps)
    ID14 = TensorMap(ones, ComplexF64, As ⊗ Ps ← Ps ⊗ As)
    ID2 = TensorMap(ones, ComplexF64, Ps ← As ⊗ Ps)
    ID24 = TensorMap(ones, ComplexF64, Ps ⊗ As ← As ⊗ Ps)
    if ifdagger
        @planar foursite[-1 -2 -3 -4; -5 -6 -7 -8] := c⁺[-1; -5 1] * ID14[1 -2; -6 2] * ID14[2 -3; -7 3] * ID1[3 -4; -8]
    else
        @planar foursite[-1 -2 -3 -4; -5 -6 -7 -8] := c[-1 1; -5] * ID24[-2 2; 1 -6] * ID24[-3 3; 2 -7] * ID2[-4; 3 -8]
    end
    return FiniteMPO(foursite)

I check it by taking the inner product：
I wish dot(ψ₀2, $c^\dagger_1 \mathbb{1}_2$ * $c_1 \mathbb{1}_2$, ψ₀2)=1 and dot(ψ₀4, $c^\dagger_1 \mathbb{1}_2 \mathbb{1}_3 \mathbb{1}_4$ * $c_1 \mathbb{1}_2 \mathbb{1}_3 \mathbb{1}_4$, ψ₀4)=1 (both ψ₀2 and ψ₀4 have one particle per site), however, I got dot(ψ₀2, $c^\dagger_1 \mathbb{1}_2$ * $c_1 \mathbb{1}_2$, ψ₀2)=1 and dot(ψ₀4, $c^\dagger_1 \mathbb{1}_2 \mathbb{1}_3 \mathbb{1}_4$ * $c_1 \mathbb{1}_2 \mathbb{1}_3 \mathbb{1}_4$, ψ₀4) $\neq$ 1.

I am not sure if this is the right way to deal with this problem and I would be very grateful for your help.

Thanks!

[lkdvos]

Hi Jason,

This is a good question, as it has quite a lot of subtle details. Let me start by mentioning that while your approach is a good attempt, but will not produce the correct results. The problem is that this auxiliary space really does change the structure of the tensors, and it is not possible to create an identity operator that has an odd auxiliary charge. The reason is that the total charge of a tensor needs to be trivial, such that this "identity" you created actually will exchange even with odd, and is actually some kind of flip operator.

I don't think we currently have a good system to deal with these kind of setups. In principle, the state you try to create is charged, i.e. the total charge is not trivial. Thus, to create these MPSs, the MPS should have an auxiliary index that carries this charge, i.e. one of the mps tensors should have an additional charged index. In MPSKit, we can absorb this charged index in the last index of the final mps tensor, or the first index of the first mps tensor, as these already are auxiliary spaces, but there is currently no automated method to create these states. I attached a small drawing to clarify what I mean.

This is definitely something that is on my to-do list, but I currently have some other things I need done first, so it might take me a while to get back to this (say a couple weeks). If you manage to implement it yourself, I would be more than happy to check, or answer additional questions to help out. I'll definitely keep this issue open, such that I am reminded to get this done.

[Jason-zo]

Hi Lukas,

Thank you for the prompt and kind response, and the ideas you conceived! I will try to resolve this issue and share with you my difficulties or progress.

[Jason-zo]

Hi Lukas!

I apologize for not being able to consider this issue for a long time due to some matters I had at hand. I create charged MPS and absorb the charge index at the finial bond of the MPS according to your idea:

function chargedMPS(state::FiniteMPS, opt::TensorMap, site::Integer)
    pspaces = [codomain(state.AL[i])[2] for i in 1:length(state)]
    @assert (length(domain(opt))==2)&&(in(domain(opt)[1],pspaces)) "Physical and virtual space should be set at the left and right side of domain, respectively"
    temp = [state.AL[i] for i in 1:length(state)]
    @tensor T[-1 -2; -3 -4] := opt[-2; 1 -3] * state.AL[site][-1 1; -4]
    temp[site] = TensorMap(T.data, codomain(T), fuse(domain(T)))
    ide = TensorMap(ones, ComplexF64, domain(opt)[2] ⊗ pspaces[1] ← pspaces[1] ⊗ domain(opt)[2])
    for i in (site+1):length(state)
        @tensor T[-1 -2 -3; -4 -5] := ide[-2 -3; 1 -4] * state.AL[i][-1 1; -5]
        temp[i] = TensorMap(T.data, fuse(codomain(T)[1],codomain(T)[2]) ⊗ codomain(T)[3], fuse(domain(T)))
    end
    return FiniteMPS(temp)
end

where opt is the creation/annihilation operator and ide is used to transfer the nontrivial quantum number to the final bond of the MPS.
For examples, we have a ground state gs:

4-site FiniteMPS (ComplexF64, GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, U1Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, U1Irrep}}, Int64}}):
┌── AC[4]: TensorMap((Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1)) ← Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1))
├── AL[3]: TensorMap((Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>4, (0, 2, 0)=>1, (0, -2, 0)=>1, (1, 1, 1)=>2, (1, -1, 1)=>2, (1, 1, -1)=>2, (1, -1, -1)=>2, (0, 0, 2)=>1, (0, 0, -2)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1)) ← Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1))
├── AL[2]: TensorMap((Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1)) ← Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>4, (0, 2, 0)=>1, (0, -2, 0)=>1, (1, 1, 1)=>2, (1, -1, 1)=>2, (1, 1, -1)=>2, (1, -1, -1)=>2, (0, 0, 2)=>1, (0, 0, -2)=>1))
└── AL[1]: TensorMap((Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1)) ← Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1))

a spin-up creation operator c⁺u:

TensorMap(Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1) ← (Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 1)=>1))):
* Data for sector ((FermionParity(1) ⊠ Irrep[U₁](1) ⊠ Irrep[U₁](0)),) ← ((FermionParity(0) ⊠ Irrep[U₁](0) ⊠ Irrep[U₁](-1)), (FermionParity(1) ⊠ Irrep[U₁](1) ⊠ Irrep[U₁](1))):
[:, :, 1] =
 1.0 + 0.0im
* Data for sector ((FermionParity(0) ⊠ Irrep[U₁](0) ⊠ Irrep[U₁](1)),) ← ((FermionParity(1) ⊠ Irrep[U₁](-1) ⊠ Irrep[U₁](0)), (FermionParity(1) ⊠ Irrep[U₁](1) ⊠ Irrep[U₁](1))):
[:, :, 1] =
 1.0 + 0.0im

and a spin-down creation operator c⁺d:

TensorMap(Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1) ← (Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, -1, 1)=>1))):
* Data for sector ((FermionParity(1) ⊠ Irrep[U₁](-1) ⊠ Irrep[U₁](0)),) ← ((FermionParity(0) ⊠ Irrep[U₁](0) ⊠ Irrep[U₁](-1)), (FermionParity(1) ⊠ Irrep[U₁](-1) ⊠ Irrep[U₁](1))):
[:, :, 1] =
 1.0 + 0.0im
* Data for sector ((FermionParity(0) ⊠ Irrep[U₁](0) ⊠ Irrep[U₁](1)),) ← ((FermionParity(1) ⊠ Irrep[U₁](1) ⊠ Irrep[U₁](0)), (FermionParity(1) ⊠ Irrep[U₁](-1) ⊠ Irrep[U₁](1))):
[:, :, 1] =
 1.0 + 0.0im

I apply the creation operators sequentially on the ground state at the 1st, 3rd, and 2nd lattice sites

cgs₁ = chargedMPS(gs, c⁺u, 1)
cgs₂ = chargedMPS(cgs₁, c⁺d, 3)
cgs₃ = chargedMPS(cgs₂, c⁺d, 2)

and obtain

4-site FiniteMPS (ComplexF64, GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, U1Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, U1Irrep}}, Int64}}):
┌ CL[5]: TensorMap(Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 1)=>1) ← Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 1)=>1))
├── AL[4]: TensorMap((Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (0, 0, 1)=>1, (0, 2, 1)=>1, (1, 1, 2)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1)) ← Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 1)=>1))
├── AL[3]: TensorMap((Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>2, (0, 2, 0)=>1, (1, 1, 1)=>2, (1, -1, 1)=>1, (1, 1, -1)=>1, (0, 0, 2)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1)) ← Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (0, 0, 1)=>1, (0, 2, 1)=>1, (1, 1, 2)=>1))
├── AL[2]: TensorMap((Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (0, 0, 1)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1)) ← Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>2, (0, 2, 0)=>1, (1, 1, 1)=>2, (1, -1, 1)=>1, (1, 1, -1)=>1, (0, 0, 2)=>1))
└── AL[1]: TensorMap((Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1)) ← Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (0, 0, 1)=>1))

for cgs₁,

4-site FiniteMPS (ComplexF64, GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, U1Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, U1Irrep}}, Int64}}):
┌ CL[5]: TensorMap(Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 2)=>1) ← Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 2)=>1))
├── AL[4]: TensorMap((Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 1)=>1, (1, 1, 2)=>1, (1, -1, 2)=>1, (0, 0, 3)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1)) ← Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 2)=>1))
├── AL[3]: TensorMap((Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>2, (0, 2, 0)=>1, (1, 1, 1)=>2, (1, -1, 1)=>1, (1, 1, -1)=>1, (0, 0, 2)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1)) ← Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 1)=>1, (1, 1, 2)=>1, (1, -1, 2)=>1, (0, 0, 3)=>1))
├── AL[2]: TensorMap((Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (0, 0, 1)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1)) ← Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>2, (0, 2, 0)=>1, (1, 1, 1)=>2, (1, -1, 1)=>1, (1, 1, -1)=>1, (0, 0, 2)=>1))
└── AL[1]: TensorMap((Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1)) ← Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (0, 0, 1)=>1))

for cgs₂, and

4-site FiniteMPS (ComplexF64, GradedSpace{ProductSector{Tuple{FermionParity, U1Irrep, U1Irrep}}, TensorKit.SortedVectorDict{ProductSector{Tuple{FermionParity, U1Irrep, U1Irrep}}, Int64}}):
┌ CL[5]: TensorMap(Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, -1, 3)=>1) ← Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, -1, 3)=>1))
├── AL[4]: TensorMap((Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, -1, 2)=>1, (0, 0, 3)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1)) ← Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, -1, 3)=>1))
├── AL[3]: TensorMap((Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1, (1, 1, 1)=>1, (1, -1, 1)=>1, (0, 0, 2)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1)) ← Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, -1, 2)=>1, (0, 0, 3)=>1))
├── AL[2]: TensorMap((Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (0, 0, 1)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1)) ← Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1, (1, 1, 1)=>1, (1, -1, 1)=>1, (0, 0, 2)=>1))
└── AL[1]: TensorMap((Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((0, 0, 0)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1)) ← Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (0, 0, 1)=>1))

for cgs₃.
I check the finial bond CL[5] are correct for these MPS, at least at the level of quantum numbers of spin and charge. I mean the accumulated quantum numbers (0,0,0)+(1,1,1) is equal to (1,1,1) for cgs₁,(0,0,0)+(1,1,1)+(1,-1,1)is equal to (0,0,2) for cgs₂ and (0,0,0)+(1,1,1)+(1,-1,1)+(1,-1,1) is equal to (1,-1,3) for cgs₃.

I thought its pretty good until I checked their inner product and compared it to the exact diagonalization (ED) results.
Their results and ED results are shown at the left and right sides of //, respectively:

dot(cgs₁,cgs₁):
0.5000000525834422  //  0.4999999999999999
dot(cgs₂,cgs₂):
0.3941804497742264  //  0.19709024478006926
dot(cgs₃,cgs₃):
0.8092385994548575  //  0.15949303621597344

As it stands, there still seem to be some errors, but it feels like it is very close to the correct result. The inner product result of cgs₂ shows that it is exactly twice the result of ED, but from the cgs₃ result, it doesn’t appear to be a simple factor of 2.

I suspect the issue might lie in two areas. First, with @tensor since replacing it with @planar or @plansor doesn’t work either; to be honest, I haven’t yet figured out the relationship between these three. Second, the problem might be with the returned FiniteMPS(temp), as charged MPS might require properties beyond those of FiniteMPS.

I kindly request your help or advice, please.

[Jason-zo]

I attached a figure to clarify how the functionchargedMPSworks.

[maartenvd]

In general I would caution against using @tensor with fermionic symmetries. If you keep it planar, you will have to be more verbose in writing down the tensor contraction but you will end up with the result that you would have expected. If you use @tensor, then you may have to compensate with parity matrices in some place.

To keep things planar, let's look at the first tensor contraction. The codomain of B1 (pre-fusing) should be: first virtual space x first physical space. The domain should be the domain of the original A1 tensor times the virtual leg of your c+ operator. But if you print the domain of the T tensor you find: (Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 1)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1)) - you have flipped those two spaces! That is why @planar complains.

The way you fuse spaces is also unsafe. There may be refined machinery in tensorkit that does it efficiently but the easiest way to get it correct is to explicitly multiply by isomorhpisms.

There is also this "identity tensor" that you insert, to pull the virtual leg all the way to the end. That tensor too may be wrong - why should filling the fused blocks with one result in an actual identity tensor? We have a special symbol in @planar for the correct braiding tensor:

function chargedMPS(state::FiniteMPS, opt::TensorMap, site::Integer)
  pspaces = [codomain(state.AL[i])[2] for i in 1:length(state)]
  @assert (length(domain(opt))==2)&&(in(domain(opt)[1],pspaces)) "Physical and virtual space should be set at the left and right side of domain, respectively"
  temp = [state.AL[i] for i in 1:length(state)]
  @planar T[-1 -2; -3 -4] := opt[-2; 1 -4] * state.AL[site][-1 1; -3] # had to switch the two legs to keep it planar
  temp[site] = T * isomorphism(domain(T),fuse(domain(T)))
  for i in (site+1):length(state)
    iso_1 = isomorphism(fuse(codomain(state.AL[i])[1],space(ide,1)),codomain(state.AL[i])[1]*space(ide,1))
    @planar T[-1 -2; -4 -5] := iso_1[-1;2 3]*τ[3 -2; 1 -5] * state.AL[i][2 1; -4] # had to switch 4/5 here as well
    temp[i] = T * isomorphism(domain(T),fuse(domain(T)))
  end
  return FiniteMPS(temp)
end

Ultimately though, I fear that your operators are simply not the correct operators.

[Jason-zo]

Hi Maartenvd,
Thank you very much for your prompt, detailed, and kind response and your response cleared up my confusion about @planar and @tensor, and reminded me of the correct use of the isomorphism and braiding tools. I also reconsidered my operator definition might be incorrect as you mentioned, and actually you are right, I find a minus sign is missing in my c⁺d operator. I fix it, and now,c⁺d and its Array form are as follow :

c⁺d
TensorMap(Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1) ← (Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, -1, 1)=>1))):
* Data for sector ((FermionParity(1) ⊠ Irrep[U₁](-1) ⊠ Irrep[U₁](0)),) ← ((FermionParity(0) ⊠ Irrep[U₁](0) ⊠ Irrep[U₁](-1)), (FermionParity(1) ⊠ Irrep[U₁](-1) ⊠ Irrep[U₁](1))):
[:, :, 1] =
 1.0 + 0.0im
* Data for sector ((FermionParity(0) ⊠ Irrep[U₁](0) ⊠ Irrep[U₁](1)),) ← ((FermionParity(1) ⊠ Irrep[U₁](1) ⊠ Irrep[U₁](0)), (FermionParity(1) ⊠ Irrep[U₁](-1) ⊠ Irrep[U₁](1))):
[:, :, 1] =
 -1.0 + 0.0im

convert(Array, c⁺d)
┌ Warning: FermionParity Arrays do not preserve categorical properties.
└ @ TensorKit ~/.julia/packages/TensorKit/Fz2UW/src/sectors/fermions.jl:53
4×4×1 Array{ComplexF64, 3}:
[:, :, 1] =
  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
  0.0+0.0im  0.0+0.0im  0.0+0.0im  1.0+0.0im
 -1.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im
  0.0+0.0im  0.0+0.0im  0.0+0.0im  0.0+0.0im

Then, I used the modified function

function chargedMPS(state::FiniteMPS, opt::TensorMap, site::Integer)
  pspaces = [codomain(state.AL[i])[2] for i in 1:length(state)]
  @assert (length(domain(opt))==2)&&(in(domain(opt)[1],pspaces)) "Physical and virtual space should be set at the left and right side of domain, respectively"
  temp = [state.AL[i] for i in 1:length(state)]
  @planar T[-1 -2; -3 -4] := opt[-2; 1 -4] * state.AL[site][-1 1; -3] # had to switch the two legs to keep it planar
  temp[site] = T * isomorphism(domain(T),fuse(domain(T)))
  for i in (site+1):length(state)
    iso_1 = isomorphism(fuse(codomain(state.AL[i])[1],domain(opt)[2]),codomain(state.AL[i])[1]*domain(opt)[2])
    @planar T[-1 -2; -4 -5] := iso_1[-1;2 3] * τ[3 -2; 1 -5] * state.AL[i][2 1; -4] # had to switch 4/5 here as well
    temp[i] = T * isomorphism(domain(T),fuse(domain(T)))
  end
  return FiniteMPS(temp)
end

to recalculate these inner product again, but unfortunately got the same incorrect results:

dot(cgs₁,cgs₁):
0.5000000088791062  //  0.4999999999999999
dot(cgs₂,cgs₂):
0.394180485365641   //  0.19709024478006926
dot(cgs₃,cgs₃):
0.8092386084268487  //  0.15949303621597344

This is a bit frustrating. Do you think this issue could be related to the FiniteMPS constructor itself? Or perhaps, at this point, the inner product is not quite like the inner product of FiniteMPS, but should be handled separately, similar to how correlator are processed?

[maartenvd]

Creating the correct operators is by far the most cumbersome part of using tensorkit (and by extension mpskit), I can definitely sympathize! I have limited information about your particular usecase, but I assume that you have spinful fermions but don't want to impose su(2) x u(1) x fermion (which would be easier, as I have already struggled through that case). With the minus sign, your operators do make sense:

cpu = TensorMap(ones,ComplexF64,Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1) ← (Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 1)=>1)))
cpd = TensorMap(ones,ComplexF64,Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1) ← (Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1) ⊗ Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, -1, 1)=>1)))
blocks(cpd)[(FermionParity(0) ⊠ Irrep[U₁](0) ⊠ Irrep[U₁](1))].*=-1;

Can you share what the groundstate should be analytically?

For a random test groundstate with one-site density matrix rho:

TensorMap(Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1) ← Vect[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])]((1, 1, 0)=>1, (1, -1, 0)=>1, (0, 0, 1)=>1, (0, 0, -1)=>1)):
* Data for sector ((FermionParity(1) ⊠ Irrep[U₁](1) ⊠ Irrep[U₁](0)),) ← ((FermionParity(1) ⊠ Irrep[U₁](1) ⊠ Irrep[U₁](0)),):
 0.05157321226771339 + 0.0im
* Data for sector ((FermionParity(1) ⊠ Irrep[U₁](-1) ⊠ Irrep[U₁](0)),) ← ((FermionParity(1) ⊠ Irrep[U₁](-1) ⊠ Irrep[U₁](0)),):
 0.7011118056268212 + 0.0im
* Data for sector ((FermionParity(0) ⊠ Irrep[U₁](0) ⊠ Irrep[U₁](1)),) ← ((FermionParity(0) ⊠ Irrep[U₁](0) ⊠ Irrep[U₁](1)),):
 0.036147179235762134 + 0.0im
* Data for sector ((FermionParity(0) ⊠ Irrep[U₁](0) ⊠ Irrep[U₁](-1)),) ← ((FermionParity(0) ⊠ Irrep[U₁](0) ⊠ Irrep[U₁](-1)),):
 0.21116780286970435 + 0.0im

You should find that

cgs₁ = chargedMPS(gs, cpu, 1)
dot(cgs₁,cgs₁)

prints

<vacuum| c ⁻d c ⁻u (c ⁻u c⁺u) c⁺u c⁺d|vacuum> *  0.036147179235762134
+<vacuum| c ⁻d (c ⁻u c⁺u) c⁺d|vacuum> *   0.7011118056268212
+<vacuum| c ⁻u (c ⁻u c⁺u) c⁺u|vacuum> *  0.05157321226771339
+<vacuum|  (c ⁻u c⁺u) |vacuum> *  0.21116780286970435
 =  0.7011118056268212 + 0.21116780286970435 = 0.9122796084965259

and indeed

dot(cgs₁,cgs₁) = 0.9122796084965259

similarly

cgs₂ = chargedMPS(gs, cpd, 1)
dot(cgs₂,cgs₂)

should - and does - reduce to:

<vacuum| c ⁻d c ⁻u (c ⁻d c⁺d) c⁺u c⁺d|vacuum> *  0.036147179235762134
+<vacuum| c ⁻d (c ⁻d c⁺d) c⁺d|vacuum> *   0.7011118056268212
+<vacuum| c ⁻u (c ⁻d c⁺d) c⁺u|vacuum> *  0.05157321226771339
+<vacuum|  (c ⁻d c⁺d) |vacuum> *  0.21116780286970435
 =  0.21116780286970435 + 0.05157321226771339 = 0.26274101513741766 

[maartenvd]

@Jason-zo actually come to think of it, you may be running into a confusion I also had when implementing the quantum chemistry hamiltonian. The contraction we are doing here is actually the one drawn in the following picture

Now let's calculate the first operator I drew there explicitly!

@planar blob[-1 -2;-3 -4] := cpu[-1;-3 1]*τ[-4 1;2 3]*conj(cpu[3;-2 2])
stringnames = Dict((FermionParity(1) ⊠ Irrep[U₁](-1) ⊠ Irrep[U₁](0)) => ".d",(FermionParity(1) ⊠ Irrep[U₁](1) ⊠ Irrep[U₁](0)) => "u.",(FermionParity(0) ⊠ Irrep[U₁](0) ⊠ Irrep[U₁](-1)) => "..",(FermionParity(0) ⊠ Irrep[U₁](0) ⊠ Irrep[U₁](1)) => "ud")
for (f1,f2) in fusiontrees(blob)
  if norm(blob[f1,f2])>1e-12
    println(stringnames[f1.uncoupled[1]]," x ",stringnames[f1.uncoupled[2]],"←",stringnames[f2.uncoupled[1]]," x ",stringnames[f2.uncoupled[2]])
    println(blob[f1,f2])
  end
end

prints

u. x .d←.. x ud
ComplexF64[1.0 + 0.0im;;;;]
ud x ..←.d x u.
ComplexF64[-1.0 + 0.0im;;;;]
ud x .d←.d x ud
ComplexF64[1.0 + 0.0im;;;;]
u. x ..←.. x u.
ComplexF64[-1.0 + 0.0im;;;;]

The operator you are actually measuring is therefore:

• c⁺u x c⁻u + 2 * c⁺u x c⁺u c⁺d c⁺d

It is possible to find the two tensors a and b such that dot(chargedMPS(...,a),chargedMPS(...,b)) indeed measures the c⁺u x c⁻u correlator (I can give them to you). You can also measure a whole bunch of inner products and work out what all the correlators are. There might be an easier approach that I can help you with. What exactly are the correlators that you want to measure?

[Jason-zo]

@maartenvd Thank you very much for your two replies. They have indeed greatly helped me in understanding symmetric MPS and resolving the issue. To be honest, it took me some time to fully digest and comprehend your insights.

To start with the conclusion, I found the cause of the incorrect results and solved the problem. Actually, I made a basic mistake when writing the chargedMPS function. I forgot that an MPS with charge has a non-trivial final bond, so when initializing temp, I forgot to contract the final bond onto the last tensor. By adding temp[end] = temp[end] * state.CLs[end] after the line temp = [state.AL[i] for i in 1:length(state)], the issue was resolved.

I would like to share with you and other tensor beginners like me how I discovered this problem. Actually, this was all thanks to your guidance in helping me understand MPS in an analytical way:

I start with a two-site system which is half filled and the total spin is equal to zero.

That is to say the final bond of cgs₁ has no effect on cgs₂, and it let me rethink about the initialization of charged MPS and let the mistake be discovered.

After that, I test

dot(cgs₂,cgs₂)
0.473606797749977 - 3.110084090640804e-18im

and now it is equal to the ED result0.4736067977499788.

And the next problem is to deal with things like dot(chargedMPS(...,a),chargedMPS(...,b)) I think with your second advice I will solve it soon.
Finally, thank you again, Maarten.

[maartenvd]

No problem! This is honestly a really tricky part of our code stack... But I don't know of a cleaner alternative that doesn't end up producing even more unexpected result in another setting.

The thing to keep in mind with fermions is that the operators you have are indeed creating excitations, but the state may pick up unexpected phases depending on other fermions elsewhere on the chain. So when you are creating overlaps between "c^+" and "c^-" on the same site, there won't be a problem, the phases have to cancel out. When you are creating overlaps where those operators act on different sites, you may get surprising results. It is essentially still the jordan wigner string, that is very much made explicit by the insertion of those braiding tensors.

[Jason-zo]

@maartenvd Hi Maarten, I'm considering dot((chargedMPS(...,a),chargedMPS(...,b)) on different sites that you mentioned above. I try to replace $c^\dagger_i$ with $c^\dagger_i \cdot Z_{i+1} \cdot Z_{i+2}\cdots Z_{N}$ to deal with jordan wigner string where $Z$ is the fermion parity operator:

function chargedMPS(state::FiniteMPS, opt::TensorMap, site::Integer)
  N = length(state)
  pspaces = [codomain(state.AL[i])[2] for i in 1:N]
  Ps = pspaces[1]
  Z = TensorMap(ones, ComplexF64, Ps, Ps)
  blocks(Z)[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])(1,1,0)] .= -1
  blocks(Z)[(FermionParity ⊠ Irrep[U₁] ⊠ Irrep[U₁])(1,-1,0)] .= -1
  @assert (length(domain(opt))==2)&&(in(domain(opt)[1],pspaces)) "Physical and virtual space should be set at the left and right side of domain, respectively"
  As = [[state.AL[i] for i in 1:N-1]..., state.AC[end]]
  @planar T[-1 -2; -3 -4] := opt[-2; 1 -4] * As[site][-1 1; -3] 
  As[site] = T * isomorphism(domain(T),fuse(domain(T)))
  for i in (site+1):N
      iso = isomorphism(fuse(codomain(As[i])[1],domain(opt)[2]),codomain(As[i])[1]*domain(opt)[2])
      @planar T[-1 -2; -3 -4] := iso[-1;2 3] * Z[-2; 4] * τ[3 4; 1 -4] * As[i][2 1; -3]
      As[i] = T * isomorphism(domain(T),fuse(domain(T)))
  end
  return FiniteMPS(As)
end

The improvement with this approach is that the wave function now differs from the exact result by just an overall phase, rather than differing by different phases in each component as it did before, but it still not correct such as it gives 0.22360679774997888 for $\langle c_1 c^\dagger_2\rangle$ of the two-site system which is negative by the ED result -0.22360679774997905.
Actually, I don't think use Z is a good idea, because the specific form of Z depends on symmetries you select. I want chargedMPS be more general. And it seems to appear redundant and conflicting with FermionParity, which might be the reason why it is fail to deal with $\langle c_1 c^\dagger_2\rangle$ ?

[Jason-zo]

In other libraries, like ITensor, the operators themselves are bosonic, so they insert the fermion parity operator for the Jordan-Wigner transformation. However, TensorKit is so nice that allows the operators to be fermionic, which should require a different approach to this problem.

[Jason-zo]

Finally, I solve this issue with FiniteMPO * FiniteMPS

iso1 = isomorphism(storagetype(operator), vspace, vspace)
iso2 = isomorphism(storagetype(operator), pspace, pspace)
@planar Z[-1 -2; -3 -4] := iso1[-1; 1] * iso2[-2; 2] * τ[1 2; 3 4] * iso2[3; -3] * iso1[4; -4]
I = isomorphism(storagetype(operator), oneunit(vspace)*pspace, pspace*oneunit(vspace))
mpo = FiniteMPO([i < site ? I : i == site ? add_single_util_leg(operator) : Z for i in 1:length(state)])
cgs = mpo*state