Product of many local operators

Question

Hello!

I'd like to ask a question about constructing the MPO operator being the product of local operators. In my problem I consider two types of bosons "a" and "b" with total number @@N = N_a + N_b@@ conserved, but @@N_a - N_b@@ is not conserved. Total number of sites is @@M = N@@.

I model my system by creating an artificial lattice with @@M = 2N@@ sites, and each two lattice sites form a unit cell where the first site corresponds to boson @@a@@ and second to boson @@b@@.

On such a lattice I can define spin operators @@S^x_i, S^y_i, S^z_i@@ as follows:

auto tmp_Sx = AutoMPO(sites);
auto tmp_Sy = AutoMPO(sites);
auto tmp_Sz = AutoMPO(sites);

for (int j = 1; j < M_sites; j+=2){
  tmp_Sx +=  0.5, "Adag",j,"A", j+1; 
  tmp_Sx +=  0.5, "Adag",j+1,"A", j; 

  tmp_Sy +=  0.5/Cplx_i, "Adag",j,"A", j+1;
  tmp_Sy += -0.5/Cplx_i, "Adag",j+1,"A", j;

  tmp_Sz +=  0.5, "N", j;
  tmp_Sz += -0.5, "N", j+1;
}
auto Sx = toMPO(tmp_Sx);
auto Sy = toMPO(tmp_Sy);
auto Sz = toMPO(tmp_Sz);

I'm interested in calculating expectation value of the @@D = \prod_j d_j@@
operator, where @@d_j = S^y_{j} + i S^z_{j}@@. I've tried the following construction (let's assume N = 2 particles):

vector<MPO> d_vec(N);  
int j = 0;
for(int i = 1; i<M; i+=2){ 
  auto tmp_d_i = AutoMPO(sites);
  tmp_d_i +=  0.5/Cplx_i, "Adag",i,  "A", i+1;
  tmp_d_i += -0.5/Cplx_i, "Adag",i+1,"A", i;
  tmp_d_i +=  0.5*Cplx_i, "N", i;
  tmp_d_i += -0.5*Cplx_i, "N", i+1;
  d_vec[j] = toMPO(tmp_d_i);    
  j = j+1;
}

auto D = nmultMPO(d_vec[1],prime(d_vec[0]), "Method=",nmultMPO_method,"MaxDim",nmultMPO_MaxDim,"Cutoff",nmultMPO_cutoff});

However, it doesn't work (I have exact-diagonalization results to compare).

Can someone point me to what I'm missing? Thanks in advance!

Best,
Marcin

Answer 1

Hi Marcin,
Thanks for the question. It's helpful to see the structure of your operator, which as I understand is a product of single-site operators. So then actually you should not have to use an MPO at all, or any machinery to sum operators or multiply MPOs together.

(1) Instead you can just apply each operator onto each tensor of your MPS to compute the product D|psi>. If you call the resulting (unnormalized) state |phi> = D|psi> you can then compute the overlap <phi|psi> to obtain the expectation value <psi|D|psi>.

(2) Alternatively you could think of D as an MPO, but which has a bond dimension of 1, so it is just a product of operators. In that approach, AutoMPO is not really the best tool to use to make the MPO. Instead it would be better to make the MPO "by hand" by setting each MPO tensor one by one.

I would recommend approach #1 above. To perform the action of D onto |psi> you could write code similar to in this "code formula" page to apply each operator d_j onto your MPS:
http://itensor.org/docs.cgi?vers=cppv3&page=formulas/mps_onesite_op
By looping over all j from 1 to N you will end up with the state |phi> = D|psi>.
Then you can just use the function inner to compute inner(phi,psi) to get your expectation value.

Hope that helps - we can discuss more if you have questions.

Miles

Comments

Hi Miles,

thanks for your reply! In fact, the local operator @@d\_j@@ is a two-site operator acting on neighboring sites, so I applied an approach from http://itensor.org/docs.cgi?vers=cppv3&page=formulas/gate, i.e.

    //|phi> = D|psi> = \prod_{j} d_j |psi>, d_j acts on {j,j+1} sites
    auto phi = psi;
    for (int j = 1; j < M_sites;  j+=2) {        
    //application of two-site gate to MPS: http://itensor.org/docs.cgi?vers=cppv3&page=formulas/gate
    phi.position(j);  
    auto d_j = 0.5/Cplx_i*( sites.op("Adag",j)*sites.op("A",j+1) - sites.op("Adag",j+1)*sites.op("A",j));
    d_j = d_j + 0.5*Cplx_i*( sites.op("N",j)*sites.op("Id",j+1) - sites.op("N",j+1)*sites.op("Id",j) );
    auto phi_j_j_plus_1 = phi(j)*phi(j+1);              
    phi_j_j_plus_1 *= d_j;
    phi_j_j_plus_1.noPrime();
    auto [U,S,V] = svd(phi_j_j_plus_1,inds(phi(j)),{"Cutoff=",1E-8});
    phi.set(j,U);
    phi.set(j+1,S*V);      
    }   

Where "Id" operator is defined by me in boson.h file as

    if(opname == "Id" || opname == "id"){
      for(auto n : range1(1+maxOcc)){
        Op.set(s=n,sP=n,1);
       }
    }

Is there something missing from my side?

Thanks for help
Marcin

---

Hi Marcin, that's helpful info. You could see why I thought d_j might act only on site j. Apart from the code itself, I would say to definitely be sure that you have taken into account that the operators on overlapping sites may not commute, so to make sure that you are applying them in an ordering that gives the result you want on paper. I'm sure you've thought of that but just mentioning for completeness.

I think you mentioned a bug in an earlier version of your comment? Hope you have found it now since I don't see a mention of that anymore.

But overall, yes this looks like a great approach for doing what you want.

Best,
Miles