Sites and operators in a purification/twin-space approach

Hello,

is there a recommended way of how to implement sites and operators in purification/twin-space/reshuffled thermofield dynamics, that is the extension of Hilbert space by a copy of the local Hilbert space

\mathcal{H} = \bigotimes_j \mathcal{H}^j \rightarrow \mathcal{H} = \bigotimes_j \mathcal{H}^j \otimes \mathcal{H}^j?

What I am currently concerned with is, that in such an approach one would have to define custom MPOs with a fermionic algebra enforced only on certain sites. E.g. for annihilation operator c_i, the standard implementation would be

c_i = \left(\bigotimes_{j=1}^{i-1}\sigma^z \right)\otimes c \otimes \left(\bigotimes_{j=i+1}^{N} I \right)

but in twin-space, next to each physical site there is an ancilla site, for which the CAR does not hold. For the operator acting on a physical site one would have

\hat{c}_i = \left(\bigotimes_{j=1}^{i-1}\sigma^z \otimes I \right)\otimes c \otimes I \otimes \left(\bigotimes_{j=i+1}^{N} I \otimes I \right)

and, conversely, an operator acting on the ancilla sites would read

\tilde{c}_i = \left(\bigotimes_{j=1}^{i-1} I \otimes\sigma^z \right)\otimes I \otimes c \otimes \left(\bigotimes_{j=i+1}^{N} I \otimes I\right).

Thanks in advance!

Rudi

So what I’ve tried so far is to define a 1D lattice of L sites, and assign to it the twin space with 2L indices (for simplicity, fermionic ones). Then, I tried to avoid using AutoMPO() or OpSum() to have full control over the parity terms (which is stupid, because OpSum() will be invoked in MPO eventually). For this, I used Prod{Op}() to construct “Twin-space operators” in the way one would construct fermionic operators naively


ITensors.op(::OpName"Id",::SiteType"Fermion") = [1 0; 0 1.]
L = 3
i = 2
j = 3
sites = ITensors.siteinds("Fermion",2*L)

function ITensor_Operator_twinSpace(op_p::String,i::Int64,s::Int64,L::Int64)::Prod{Op}
    A = Prod{Op}()
    k = 1
    for j in 1:i-1
        if s == 1 # Physical space
            A *= "F",k
            A *= "Id",k+1
            k += 2
        else # Ancilla space
            A *= Op("Id",k)
            A *= "F",k+1
            k += 2 
        end
    end
    if s == 1
        A *= "$(op_p)",k
        A *= "Id",k+1
        k += 2
    else
        A *= "Id",k
        A *= "$(op_p)",k+1
        k += 2 
    end
    for j in i+1:L
        A *= "Id",k
        A *= "Id",k+1
        k += 2
    end
    return A  
end

With this function one can create “single body operators”. Next, one needs a routine to properly take the product of multiple “single body operators” and convert the product back to an MPO

function ITensor_Prod(ops,k)
    if length(ops) == 1
        return ITensors.terms(ops[1])[k]
    else
        return ITensors.terms(ops[1])[k] * ITensor_Prod(ops[2:end],k)
    end
end

function ITensor_Prod_twinSpace(ops::Vector{Prod{Op}},sites::Vector{Index{Int64}})::MPO
    m = length(ops[1])   
    M = Vector{ITensor}()
    for j in 1:m
        Mp = MPO(ITensor_Prod(ops,j),sites)[j]
        push!(M,Mp)
    end
    return MPO(M)        
end

I hoped that this will fix the initial issue, which is that single fermionic operators (or simply parity-odd terms) are not supported yet. The same error message is obtained when executing

cdag_ts_p = ITensor_Operator_twinSpace("Cdag",i,1,L)
c_ts_p = ITensor_Operator_twinSpace("C",j,1,L)
ops = [cdag_ts_p,c_ts_p]

print(ITensor_Prod_twinSpace(ops,sites))

The problem is obviously term at index 3 as it is a term which is c \cdot \sigma^z with odd parity.

Is there a way to circumvent this limitation of only having parity-even terms? I checked the

Hi, could you explain why you are using OpSum to make single fermion operators? (If I’m understanding correctly?) If you are just wanting to act with, say, a single c^\dagger or c operator a better approach would be to just call the op function to obtain these as an ITensor and just act with that. You’d then also need to act with the Jordan-Wigner string operator “F” on all sites to the left of the c or c^\dagger.

(It gets a bit more complicated if spin is involved, though not too much more.)

Hi Miles,

for giving a little bit of context on the initial question: I would like to simulate the time-evolution (using TDVP) of a vectorized density matrix in twin-space. The „super-Hamiltonian“ in such a description may have terms which act on both physical and ancilla sites, which effectively describe dissipative processes. However, this super-Hamiltonian may contain terms which have odd parity, which seems not to be allowed by the current MPO constructors. And one has to take into account that fermionic operators in twin-space do not fulfill the usual algebra fermionic operators in Hilbert space do. As described in my initial post, for an operator acting on a physical (ancilla) index, F operators have to be placed only on the physical (ancilla) sites, and identities on the ancilla (physical) site.

Hey Miles,

thanks for your suggestion. Using

ITensors.MPO(sites::Vector{Index{Int64}}, ops::Vector{String})

works well for my purpose and leads the correct Hamiltonian (checked against the actual Kronecker product). However, when using ITensorTDVP, the maximum bond dimension of the time-evolved state is always one, which I find quite strange.