Construct outer product of mixed-type delta's as an MPS

Hi,

Given the following (mixed) sites (2 spin + 2N bosons):

using ITensors
sites = append!(siteinds("S=1/2",2), siteinds("Boson", N * 2 ; dim=d)

I’d like to construct the following MPS:

|\Psi\rangle = \left(|00\rangle + |11\rangle\right)\otimes \left(|00\rangle + |11\rangle +\cdots |d-1,d-1\rangle\right)\otimes\cdots\otimes \left(|00\rangle + |11\rangle +\cdots |d-1,d-1\rangle\right).

I am wondering what is the best way to do this in iTensor.

A little bit background: the vector |\Psi\rangle is the left vector used in third quantization, the inner product of which is equivalent to taking trace in the operator space.

I found a way to manually do this. Admittedly not very elegant, but it seems to be doing the trick. Please let me know what you think and if there is a more elegant / efficient way to do this…

N = 10
dim = 4
sites = append!(siteinds("S=1/2",2), siteinds("Boson", N * 2 ; dim=dim))

vac_left = MPS(sites,"0")

spin_id = productMPS(sites[1:2],["0","0"]) + productMPS(sites[1:2],["1","1"])

vac_left[1] = spin_id[1]
s2 = Array(spin_id[2], siteind(spin_id,2), linkind(spin_id,1))
s2 = reshape(s2, 1,2,2)
vac_left[2] = ITensor(s2, linkind(psi1,2), siteind(spin_id,2), linkind(spin_id,1))

for i in 1:N
    boson_id = productMPS(sites[2*i+1:2*i+2],["0","0"]) 
    for d in 2:dim
        boson_id = boson_id + productMPS(sites[2*i+1:2*i+2],[string(d-1),string(d-1)]) 
    end
    t1 = reshape(diagm(ones(dim)),1,dim,dim)
    vac_left[2*i+1] = ITensor(t1, linkind(psi1,2*i), siteind(boson_id,1), linkind(boson_id,1))
    if i<N
        vac_left[2*i+2] = ITensor(t1, linkind(psi1,2*i+2), siteind(boson_id,2), linkind(boson_id,1))
    else
        vac_left[2*i+2] = boson_id[2]
    end
end