Getting RDM of a RDM Eigenvector

Question

Hi, I have a similar question.

My goal is to calculate the entanglement entropy for a bipartition of an eigenstate of the RDM of a bipartition of the system.

In other words, I have a system of size L and I want to find the RDM of the subsystem [0,La]. I would then like to find an eigenstate of that system (corresponding to the largest eigenvalue) and find the entanglement entropy of the sub-sub-system [0,Lb], (Lb < La).

If I understand correctly I should do something like this?:

psi.position(La); //gauge mps
auto wf1 = psi.A(La)*psi.A(La+1); 
auto spectrum1 = svd(U,S,V); //U gives me the RDM of the [0,La]
auto ir = commonIndex(psi.A(La),psi.A(La+1))
auto selt1 = setElt(ir(1));
auto rdmEvec = U*selt1; //this gives me the eigenvector corresponding to largest eigenvalue
rdmEvec.position(Lb);
auto wf2 = rdmEvec.A(Lb)*rdmEvec.A(Lb+1);
auto spectrum = svd(wf2,U,S,V); //gives me spectrum I want

Thank you in advance for the help.

Answer 1

Good question. Yes, the sample code you wrote is very much on track with what you should do. The main issue you might run into with that code, though, is that rdmEvec will not be an MPS, but just a single tensor, so you can't call .position on it (though I know what you were meaning by writing that in your code).

Instead, what you should do is prepare a new MPS object of size La, and set the tensors of it to be those of your original MPS together with the modified tensor at the edge where you "cut" and project onto the RDM eigenvector.

Then you can call .position on this new, shorter, MPS and do the rest of the calculation the way you wrote it.

I'd definitely suggest testing on a system where you already know the correct answer to ensure everything is working! Or to test in some other way before just assuming the code is correct, even if it runs ok.

Comments

Hi, thanks for the reply. I have tried what you suggested  but my eigenvalues are larger than 1. This is the code I am using:
    //Create new MPS of size La
         SpinHalf sites(La);
         MPS psi2 = MPS(sites);
 
         //Create Tensor to pick out RDM eigenvector
         psi.position(La);
         Index ir = commonIndex(psi.A(La),psi.A(La+1));
         ITensor selt1 = setElt(ir(1));
 
         ITensor wf = psi.A(La)*psi.A(La+1);
         ITensor U = psi.A(La),S,V;
         svd(wf,U,S,V);
 
         ITensor rdmEVec = U*selt1;
         for(int i=1; i<La; i++){
                 psi2.setA(i,psi.A(i));
         }
         psi2.setA(La,rdmEVec);
         psi2.position(Lb);
         ITensor wf2 = psi2.A(Lb)*psi2.A(Lb+1);
 
         ITensor U1=psi2.A(Lb),S1,V1; //set last tensor to the rdm eigenvector
 
         auto spectrum = svd(wf2,U1,S1,V1);
 
         printf("T=%d", time);
         for(auto p: spectrum.eigs()){
                 printf("\t%.12e", p);
         }

---

The code looks good, so I’m not sure why you are getting an eigenvalue greater than 1. My starting point, though, would be to assume that all of the ITensor components are working correctly (since we have tested them so much), so it’s more likely a logical error with the procedure you are trying to do.

So I’d recommend putting small tests throughout your code and printing out various things. Like are all of the tensors of your new psi2 MPS left orthogonal except the last one? What is the norm of the last tensor of psi2?

---

I believe the issue was the common index was wrong. I need it between U and S not between psi.A(La) and ps.iA(La+1). This is at least giving me proper eigenvalues. I still need to verify the results with ED. Thank you for the help.