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

**URL:** <https://itensor.discourse.group/t/obtain-2-nd-renyi-negativity-without-calculating-reduce-density-matrix-explicitly/2627>\
**Category:** ITensor Julia Questions\
**Created:** [July 2, 2026, 1:24pm UTC](https://itensor.discourse.group/t/obtain-2-nd-renyi-negativity-without-calculating-reduce-density-matrix-explicitly/2627 "2026-07-02T13:24:24Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![Wang](https://avatars.discourse-cdn.com/v4/letter/w/edb3f5/32.png) [@Wang](https://itensor.discourse.group/u/Wang)\
**Post date:** [July 2, 2026, 1:24pm UTC](https://itensor.discourse.group/t/obtain-2-nd-renyi-negativity-without-calculating-reduce-density-matrix-explicitly/2627/1 "2026-07-02T13:24:24Z")

</div>

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](https://github.com/ryanlevy/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 )

```julia
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 by`tr(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](https://global.discourse-cdn.com/free1/uploads/itensor/original/1X/4a2b7d8cef7663892d580a286828e6a2616d919d.png)

```julia
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:

```julia
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))`:

```julia
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

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

```

#### Tensor form

```julia
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？