Bug in multiplying MPOs and tracing MPOs

ITensor Julia Questions
,
1

Hi admin,

I would like to implement a code that, given an MPS | \psi_{MPS} \rangle, and assuming that I have correctly obtained density matrix \rho = | \psi \rangle \langle \psi | and the reduced density operator \rho_A for some bipartition A/B (I use 60 sites below and use exactly half bipartition), calculate the quantity

Tr( \rho_{A, Q}^2 ) = \int_{- 2 \pi}^{2 \pi} d \lambda ( - \frac{| \lambda |}{(2 \pi)^2} + \frac{1}{2 \pi} ) Tr( \rho_A e^{i \lambda Q} \rho_A e^{- i \lambda Q} )

where Q = \sum_{i=1}^{L_A -1} \sigma_i^z \sigma_{i+1}^z should ideally be an MPO so that we can apply its exponential to the reduced density matrix \rho_A effectively.

I have 2 questions:
(1) I have error when I’m trying to apply MPOs e^{i \lambda \sigma_i^z \sigma^z_{i+1}} to \rho_A. I think I cast both of them as MPOs but I suspect my code earlier might have error.

(2) For full trace, is there a built in function in iTensor Julia that takes in a final MPO and spits out the trace?

(3) Is my way, i.e., orthogonalize and then contract (for both reduced density operator and full trace) the most efficient way?

The code:

# (1) Construct Half Chain Reduced Density Matrix via contraction
    
# Orthogonalize so that we are in left-right canonical form
orthogonalize!(psi, 30)
    
# Density matrix as outer product of mixed canonical MPS
rho = outer(psi', psi)
sites = siteinds(rho)
    
L = ITensor(1.0)
for i=31:N
        # Need to check if delta is needed here
        L *= rho[i] * delta(sites[i][1], sites[i][2])
        
        # This printing check no of indices as sanity check
        # @show inds(L)
        # pause()
end
    
# Absorb the partial trace result to site 30
rho[30] *= L
            
# Cast the partial traced rho into MPO form
 rho_reduced = MPO(30)
 for i=1:30
     rho_reduced[i] = rho[i]
 end
    
# (2) Calculate Half Chain Symmetry Resolved Renyi Entropy
# via numerical integration
    
 Renyi_Sym = 0.0
 lambda = -2*pi
 prev_tr = 0.0
 curr_tr = 0.0
 new_rho_reduced = copy(rho_reduced)
    
# k = No of integration steps/splits
for k=1:101
    # Apply exponential of Q operator to rho_reduced for lambda
    # via rho_A e^(-i lambda Q) rho_A e^(i lambda Q)
            
    temp1_rho_reduced = deepcopy(rho_reduced)
    temp2_rho_reduced = deepcopy(rho_reduced)
        
    # Apply MPO of e^{i sigma^z_i sigma^z_{i+1}} for sites in A
    for alpha=1:29
        s1 = s[alpha]
        s2 = s[alpha+1]
        os = 4 * op("Sz", s1) * op("Sz", s2)
        szszp = exp(1.0im * lambda * os)
        szszm = exp(-1.0im * lambda * os)
        SzSzp = MPO([szszp])
        SzSzm = MPO([szszm])
            
        # Error here
        temp1_rho_reduced = apply(SzSzp, temp1_rho_reduced)
        temp2_rho_reduced = apply(SzSzm, temp2_rho_reduced)
    end
                    
    new_rho_reduced = apply(temp1_rho_reduced, temp2_rho_reduced)
    orthogonalize!(new_rho_reduced, 30)
    temp_sites = siteinds(new_rho_reduced)
            
    # Full trace for lambda
    for site=1:30
        new_rho_reduced[site] = new_rho_reduced[site] * delta(temp_sites[site][1], temp_sites[site][2])
    end
        
    curr_tr = (- (abs(lambda)/(2*pi)^2) + (1/(2*pi)))*(new_rho_reduced)

    # Perform trapezoid integration
    if k>1
        Renyi_Sym += 0.5*(prev_tr + curr_tr)*((2*pi)/50)
    end
            
    # Update
    lambda += (2*pi)/100
    prev_tr = curr_tr
end
    
Renyi_Sym = - log2(Renyi_Sym)
2

Sorry for the very slow reply from us. (One tip also about asking questions on here is to just ask them one at a time.)

About question (2): yes, you can use the function tr to trace an MPO. It was not documented so I have just submitted a PR on Github to add documentation for it.

When you ask if (3) is the most efficient way, what objective are you referring to? (I.e. most efficient way to do which thing?)

3

Hi Miles, thanks for the reply. I managed to resolve the bug, after being careful with priming the outer products, contracts, etc. Now looking back I realised question (3) makes no sense.
I ended up using tr function like you said, though I will also appreciate a documentation of the tr function.