How to calculate two-operator correlation functions with conserved quantum numbers

Question

Hi! I am doing some simulations with conserved quantum numbers. I am interested in computing the two-operator correlation functions. I was basing on this tutorial:

http://itensor.org/docs.cgi?vers=cppv3&page=tutorials/correlations

My code for the Sxi Sxj correlation looks as follows

auto first = [...] // This is just index of the first site
auto second = [...] // This is just index of the other site

//
// Calculate the ⟨Sx_i Sx_j⟩ (nearest-neighbor spin-spin interaction).
//

// Re-gauge psi to get ready to measure at 'first' position.
psi.position(first);
// Get the site index without QNs
auto sj_first = removeQNs(sites(first));
auto sj_second = removeQNs(sites(second));
// Create a SpinOneSite from the Index
// without QNs to use in the op function
auto op_first_x = op(SpinOneSite(sj_first),"Sx");
auto op_second_x = op(SpinOneSite(sj_second),"Sx");

// Create the bra/dual version of the MPS psi.
auto psidag = dag(psi);

// Prime the link indices to make them distinct from the original ket links.
psidag.prime("Link");

// Index linking 'first'-1 to 'first':
auto li_1 = leftLinkIndex(psi,first);

auto C = prime(psi(first),li_1)*op_first_x;
C *= prime(psidag(first),"Site");
for (int k = first+1; k < second; ++k) {
    C *= psi(k);
    C *= psidag(k);
}
// Index linking 'second' to 'second'+1:
auto lj = rightLinkIndex(psi,second);

C *= prime(psi(second),lj)*op_second_x;
C *= prime(psidag(second),"Site");

auto spin_spin_inter_x = elt(C); // <-- final value that I wanted to compute

Is my approach correct? Or maybe I should also remove QNs on the sites between first and second?

Answer 1

Hi, so instead of removing QNs at all, I would recommend using the fact that Sx = (S+ + S-)/2. Then if you measure four correlation functions:

<S+ S+>, <S+ S->, <S- S+>, <S- S->

, you can recombine those numbers to get

<Sx Sx>

. It doesn’t require any extra code if you wrap your measurement code in a function and pass the operator names into this function in the four combinations above.

Does that sound good to you?

Miles

Comments

That's a great note! I have two follow up questions:

1) Does the relation Sx = (S+ + S-)/2 still hold, when we calculate first the expectation values, and then just put them into the formula? So, is <Sx> = (<S+> + <S->)/2 is a correct formula?

2) How would a similar one look in the case of two-site correlators? Is it something like this <Sx_i Sx_j> = (<S+_i S+_j> + <S-_i S-_j>)/2, or is it more complicated?

---

1) For single site expectation values, yes you can just calculate them separately as you show.

2) For multi-site correlation functions, you would need to expand it out, i.e. <Sx_i Sx_j> = (<S+_i S+_j> + <S+_i S-j> + <S-_i S+_j> + <S-_i S-_j>)/4.

-Matt

---

That's fantastic! Thank you!

---

Ok glad it’s making sense! Sorry my answer got formatted very strangely, but I just fixed it.

Yes, as Matt’s answer explains, if you just replace Sx with (S+ + S-)/2 and then expand all the terms using the linearity of correlation functions, then you get the expression he wrote (it is just a linear function because the wavefunction is a vector and the operators are linear operators on the vector space).

Best regards,
Miles’

---

One more question, how can I calculate the different-time correlator with QN breaking operators like green function? Here the two operators are seperated by the time-evolution.

---

Hi, thanks for asking, though it is a very general question and I'm not sure which part you have a question about. Generally speaking you compute such a correlator <0|A(t) B|0> by first acting with B on the ground state |0> (in the case of zero temperature), then time evolving the state B|0>, and act with A. I'm glossing over some important details such as using the fact that |0> is an eigenstate in various places, and using some nice tricks such as actually time evolving by an amount t/2 on the states B|0> and A|0> separately to save on computational costs. There are many papers explaining different schemes for doing this, and all of the steps can be done with ITensor library tools (e.g. see our sample codes on doing time evolution or acting an operator on an MPS).

Does that help?

Miles

---

Thank you for your answer, Miles. ^_^
I have problems on Green function calculation.
When I define the B as "Cdagup", and A as "Cup" in t-J model, it says

I = (2|id=414|l=32,Link) <Out>
  1: 2 QN()
Q = QN({"Nf",-1,-1},{"Sz",-1})
From line 682, file index.cc

Index does not contain given QN block.

Index does not contain given QN block.


and here's my code
    auto ampoA = AutoMPO(sites);
    ampoA += "Cdagup", N/2;
    auto A = toMPO(ampoA); //QNbug

---

Hi, good question. So unfortunately here you have run into a known limitation of the AutoMPO system, which is that when QN conservation is turned on it has a problem making sums of just a single non-diagonal operator. While this is something we should fix, the good news is that there is a much more efficient way to apply this operator to your state that you should definitely use instead.

This better way, which we could discuss more, is to obtain the operator first as an ITensor like this:
auto cdup = sites.op("Cdagup",N/2);

then you can apply it to the N/2 site of your MPS by writing similar code as in the following example here:
http://itensor.org/docs.cgi?vers=julia&page=formulas/mps_onesite_op

Does that allow you to do what you want? It is a better (more efficient) way to apply a single-site operator to an MPS than using an MPO would be.

Best,
Miles