1-RDM from Julia version

[Manas]

This is about the question of computing a RDM from an MPS here. Given an MPS over n spin sites , I want to compute the 1-RDM of the kth site only (all n-1 sites with be traced except site index k). I want to retrieve the 1-RDM and its spectrum thereafter. I followed the post below

How to calculate flux of a reduced density matrix?

Hi Miles, Sorry for the late update. It seems i have some troubles for obtaining rdm and its factorization. For obtaining rdm from a MPS |\psi\rangle, a naive way is to call outer(psi',psi) (for example, How to get density matrix operator from MPS wave function? - ITensor Support Q&A) and partially traced out half of the system. (for example, the discussion here trace and partial trace of MPO - #3 by jisutich) However, this method seems to be very slow when the bond dimension is large. The…

However was not able to understand the last response completely. I understand the orthogonalize part but not specifically the lines beyond as mentioned here

How to calculate flux of a reduced density matrix?

Hi Zach, Thanks for the follow-up question here. So for the case of a pure state which is an MPS, there is a shortcut to getting the kind of rdm you want (where sites 1:k are traced over). The shortcut is to “gauge” the MPS so that the tensors on sites 1,2,…,k are all “left-orthogonal”. This can be done in ITensor by calling orthogonalize!(psi,k+1) on an MPS psi. Afterward, the “orthogonality center” will be at site k+1, meaning that the tensors 1,2,…,k will all be left-orthogonal (and tensors k…

What did I try ?

The following is the code I tried

#Function for density matrix
#***************************
function density_matrix(sites, psi::MPS)
    psi_diag=dag(prime(psi,sites))
    prime!(psi_diag, "Link")
    return outer(psi,psi_diag)
end

#Computing 1st order reduced density matrix at each site
#*********************************************************
function one_RDM(sites, state, index)
    rho=density_matrix(sites, state)
    #rho_tilde = copy(rho)
    #orthogonalize!(rho, index)
    for (i,s) in zip(sort!(setdiff(collect(1:length(sites)), index)), sites)
        rho[i] *= delta(s, prime(s))
    end

    L = ITensor(1.0)
    L = prod(rho)   # conversion to Itensor object to apply eigen methods
    #for i in 1:sort!(setdiff(collect(1:length(sites)),index))
    #    L *= rho[i]      #same as prod above
    #end
    @show(tr(L))
    return L
end

It works , but is remarkably slow for large bond dimensions. Following the comment in the above linked post, I am not sure how to proceed beyond the following :

using ITensors
n = 6
s = siteinds("S=1",n;conserve_qns=true)
state = ["Dn" for i = 1:n]
psi = randomMPS(s,state)
k = 4  # k is the physical index where i want the 1RDM not traced (the ^1\rho^k_k' are the indices)

orthogonalize!(psi, k)
for i in k+1:length(s)
     psi[i] = psi[i]*prime(dag(psi[i]),siteinds(psi))
     # What needs to be done next following the discussion to get correct 1-RDM and its spectrum ?

Any help would be appreciated on how to efficiently compute the quantity above (even for larger bond dimensions) Also what would change if I wanted to retrieve the 2-RDM at say site l and m. I can find a C++ post about this ITensor but not a Julia post. Certain commands like rightLinkIndex or leftLinkIndex has no analogue it seems

[miles]

Hi Manas, here are some notes about how to compute the 1-RDM and the reasoning for the code. It helps to draw it as diagrams:

The main points are:

• you only need to prime the site index, as you want to contract over the Link indices
• you should just use the * operator (ITensor contraction) and not the outer function

Regarding your other question, yes in Julia there’s no analogue of leftLinkIndex and rightLinkIndex currently. For simplicity there is just linkind(psi,j) which gives you the Link index between tensors j and j+1.