Obtain 2-nd Renyi negativity without calculating reduce density matrix explicitly.

ITensor Julia Questions
1

Hello everyone!

I am trying to compute the second-order Renyi negativity E_2=-\log(\text{tr}(ρ_{AB}^{T_B})^2/\text{tr}(ρ²)) of a subsystem from a given MPS wavefunction, based on the framework of ITensorEntropyTools.jl. The standard approach uses density_matrix_sites to build a high-rank reduced density matrix (RDM) tensor, then performs a partial transpose and computes E_2. However, I am running into two issues:

  1. A tensor-by-tensor contraction method (without building the full RDM) gives a different result from the RDM-based calculation.
  2. An MPO representation of the RDM behaves differently under partial transpose compared to the high-rank tensor representation.

Code and Issues

1. RENYI-2 negativity function (based on ITensorEntropyTools.jl )

function renyi_2_nega(psi::MPS, st::Int64, en::Int64, region::Vector)
    sites_p = siteinds(psi)
    rho_st_en = density_matrix_sites(psi, collect(st:1:en))
    tr_rho_2 = tr(inner(rho_st_en, rho_st_en))
    for ind in eachindex(region)
        j_swap = region[ind]
        swapinds!(rho_st_en, sites_p[j_swap], prime(dag(sites_p[j_swap])))
    end
    return -log(tr(apply(rho_st_en, rho_st_en)) / tr_rho_2)
end

Question 1: I also attempted to compute \text{tr}(ρ_{AB}^{T_B})^2 by contracting tensors one by one (without building the full RDM), which I expected to be equivalent. However, the result from this contraction disagrees with the result obtained by building the RDM via density_matrix_sites and where the 2-nd Renyi negativity is obtained bytr(apply(rho_2, rho_2)) after applying partial transpose swapinds!(rho_2,siteind(psi_g,3),dag(prime(siteind(psi_g,3)))). What could cause this discrepancy? Is there a subtlety in the index contraction order or in how inner vs apply handles the indices? The code to calculate 2-nd Renyi negativity is posted below and also the schematic diagram.

image
image2820×1692 305 KB

function trace_rho_sq_pt(psi::MPS, subsystem::Vector{Int}, transposed_sites::Vector{Int})
    N = length(psi)
    sites = siteinds(psi)
    first = subsystem[1]
    last = subsystem[end]

    site_tensor=1.0
    for n_ind in eachindex(subsystem)

        n=subsystem[n_ind]

        k1 = psi_g[n];  b1 = dag(prime(psi_g[n], 1))
        k2 = prime(psi_g[n], 2);  b2 = dag(prime(psi_g[n], 3))
        si = siteind(psi_g, n)

        if n==1
            if first > 1
                l = linkind(psi, first - 1)
                site_tensor *= delta(l, prime(l, 2)) * delta(prime(l, 1), prime(l, 3))
            end
        end
        
        if n in transposed_sites
            d1 = delta(dag(si), dag(prime(si, 2)))
            d2 = delta(prime(si, 1), prime(si, 3))
            site_tensor *= (k1 * d1 * k2) * (b1 * d2 * b2)
        else 
            d1 = delta(dag(si), prime(si, 3))
            d2 = delta(prime(si, 1), dag(prime(si, 2)))
            site_tensor *= (k1 * d1 * b2) * (b1 * d2 * k2)
        end

    end

    
    if last < N
        r = linkind(psi, last)
        site_tensor *= delta(dag(r), prime(r, 1)) 
        site_tensor *= delta(dag(prime(r, 2)), prime(r, 3))
    end

    return real(scalar(site_tensor))
end

2. MPO representation of RDM

Because density_matrix_sites(psi, collect(i:j)) constructs a high-rank tensor that can be very memory-intensive when i and j are far apart, I constructed an MPO-form RDM:

function reduced_density_MPO(psi::MPS, subsystem::Vector{Int})
    sorted_sub = sort(subsystem)
    first_sub = first(sorted_sub)
    last_sub = last(sorted_sub)
    n_sub = length(sorted_sub)

    orthogonalize!(psi, first_sub)

    subsys_sites = [siteind(psi, i) for i in sorted_sub]
    rho = MPO(subsys_sites)

    for (j, i) in enumerate(sorted_sub)
        A = psi[i]
        rho[j] = A * dag(prime(A, 1))

        if j == 1 && first_sub > 1
            ll = linkind(psi, first_sub - 1)
            rho[1] *= delta(ll, dag(prime(ll, 1)))
        end

        if j == n_sub
            rl = linkind(psi, last_sub)
            rho[n_sub] *= delta(dag(rl), prime(rl, 1))
        end
    end

    return rho
end

This works correctly for \text{tr}(\rho_{AB}^2), where tr(inner(rho_1mrho_1))=tr(inner(rho_2,rho_2)):

rho_1 = reduced_density_MPO(psi_g, [1,2,3])
rho_2 = density_matrix_sites(psi_g, [1,2,3])

3. Partial transpose on MPO vs tensor form

However, after partial transpose on site 3:

MPO form

swapinds!(rho_1[3], siteind(psi_g, 3), dag(prime(siteind(psi_g, 3))))
tr(apply(rho_1, rho_1))

Tensor form

swapinds!(rho_2, siteind(psi_g, 3)dag(prime(siteind(psi_g, 3))))
tr(apply(rho_2, rho_2))

Question 2:For the tensor form, after partial transpose, tr(apply(rho_2, rho_2)) ≠ tr(inner(rho_2, rho_2)). For the MPO form, tr(apply(rho_1, rho_1)) =tr(inner(rho_1, rho_1)) . After partial transpose, the RDM is no longer Hermitian, so these two quantities should differ. What is the difference between apply ( ) acting on higher-order tensors and MPOs resulting this inconsistence？