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

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.

// Index linking 'first'-1 to '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:

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?

+1 vote
edited by

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

commented by (300 points)
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?
commented by (14.1k points)
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
commented by (300 points)
That's fantastic! Thank you!
commented by (70.1k points)
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’
commented by (160 points)
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.
commented by (70.1k points)
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
commented by (160 points)
I have problems on Green function calculation.
When I define the B as "Cdagup", and A as "Cup" in t-J model, it says

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
commented by (70.1k points)
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