Best way to generate the MPO of an entangled state

ITensor Julia Questions
1

Hey there,
I’m trying to generate an initial state composed of the tensor product of L Bell states |\psi\rangle = (\frac{|00>+|11>}{\sqrt{2}})^{\otimes L}, and then obtain the MPO of its density matrix \rho = |\psi\rangle\langle\psi|. Initially, I tried the following code:

using ITensors
Bell = [1; 0; 0; 1]/√(2)
state = 1
for i = 1:L
    global state = kron(state, Bell)
end
Idx_list = siteinds("Qubit", 2*L)
psi = MPS(state, Idx_list; cutoff = 1E-8, maxdim = 1000)
rho = outer(psi,psi')

It works fine when L<7, but displays a warning when I set L=8:

Contraction resulted in ITensor with 16 indices, which is greater than or equal to the ITensor order warning threshold 14. You can modify the threshold with macros like @set_warn_order N, @reset_warn_order, and @disable_warn_order or functions like ITensors.set_warn_order(N::Int), ITensors.reset_warn_order(), and ITensors.disable_warn_order().

I tried ITensors.disable_warn_order() but the system automatically kills the run due to large memory.

I’m not sure if I can create an empty MPO first and assign each element to the corresponding ITensor of the Bell states (in the density matrix form), as each contains 4 indices and is entangled. I could also start with a product state, and then apply two-site gates to evolve it to the state that I want, but that’ll be too much trouble. So I’m wondering what the best way is.

Thank you so much for your time!

2

I tried the last method that I mentioned and succeeded, below is my code:

using ITensors
L=2
Idx_list = siteinds("Qubit", 2*L)
P_u = [1 0; 0 0]
rho_0 = MPO(2L)
# Initialize a product state |000...>
for j=1:2L
    T = op(P_u, Idx_list[j])
    rho_0[j] = T
end
# Map from |00> to Bell state (|00>+|11>)/√2
U = [1 0 0 0; 0 0 0 0; 0 0 0 0; 1 0 0 0]./√2
for i=1:L
    x_1 = 2i-1
    x_2 = 2i
    U_B = op(U,Idx_list[x_1],Idx_list[x_2])
    rho_0 = apply(U_B,rho_0,apply_dag=true; cutoff=Cutoff, maxdim=Maxdim)
end

Still curious to know if there’s a simpler way though.

3

What you did is fine, although you can initialize the MPO to be |0\rangle\langle0| directly with MPO(Idx_list,"Proj0").

If you combine your two codes, you’ll get you’re after I think

using ITensors, ITensorMPS

L=10 # number of pairs
cutoff = 1e-10
maxdim = 10
sites = siteinds("Qubit", L*2)
# Initialize a product state |000...>
psi = MPS(sites,"0")
# Map from |00> to Bell state (|00>+|11>)/√2
U = [1 0 0 0; 0 0 0 0; 0 0 0 0; 1 0 0 0]./√2
# for each pair of sites, apply U to them
for i=1:2:L*2
    x_1 = i
    x_2 = i+1
    U_B = op(U,sites[x_1],sites[x_2])
    global psi = apply(U_B,psi; cutoff, maxdim)
end
@show maxlinkdim(psi) # should be 2
ρ_0 = outer(psi',psi)
@show maxlinkdim(ρ_0) # should be 4

In general if there is a low bond dimension |\psi\rangle that you want to form \rho = |\psi\rangle\langle\psi| then doing outer(psi',psi) should work, or apply the unitary gates.

1 Like
4

Thank you Ryan! The “Proj0” matrix is very useful. It’s interesting that the warning of contracting a large number of indices in the outer() function doesn’t appear this time.

1 Like
5

Note that, at least for when I run your code, the warning was coming from psi = MPS(state, Idx_list; cutoff = 1E-8, maxdim = 1000) not outer()

1 Like
6

This topic was automatically closed 10 days after the last reply. New replies are no longer allowed.