Efficient calculation of magnetization susceptibility

Question

Hi Miles,

I am trying to calculate the magnetization susceptibility for a quantum ferromagnet and in the process, I need the expectation value of $M_N^2$, where $$M_N \equiv \frac{1}{N} \sum_i S_i.$$

This implies that one requires the expectation value of $S_i S_j$ for all $i, j$. Although this scales terribly ($\sim L^4$) for a 2D system, in principle, I understand how to compute this using the same methods documented for computing two-point correlation functions.

However, I was wondering if there is a more efficient way of computing this. In particular, suppose I define an operator $$O \equiv \frac{1}{N^2} \sum_{i,j} S_i S_j,$$ which is nothing but $M_N^2$. Then, is it possible to simply construct this full operator as an MPO and contract it with the MPS ground-state yielded by DMRG? If so, I am not quite sure of how to gauge the MPS correctly.

Thanks a lot for your time.

Best,
RS

Answer 1

Hi RS,
Good question. Yes: making an MPO is definitely the way to go. I've used this technique before and it works very well.

In fact, there is some code here you can use or copy to make the MPO:
https://github.com/emstoudenmire/finiteTMPS/blob/master/S2.h

(It is ITensor v2 code, which shouldn't be hard to upgrade to v3 if you read our guide, or we can chat about it).

For measuring expectation values of MPOs, there's no need to think about the gauge of the MPS. The reason is that the algorithm to compute the expectation value contracts the entire MPS with the MPO with another (entire) copy of the MPS. So the gauge simply plays no role.

In ITensor (v3) you can use the function inner(psi,W,psi) where W is your MPO to compute <psi|W|psi>. Use innerC if your MPS is complex-valued.

Best regards,
Miles

Comments

Hi Miles,

Thank you very much for your quick reply and for sharing the link to your code—it is really helpful to have that as a reference!

At the moment, since I am not dealing with conserved quantum numbers, I was thinking of something much simpler:

    auto sites = SpinHalf(N); //make a chain of N spin 1/2's
    // sites is the same site-set that I used in constructing my Hamiltonian
    auto ampo = AutoMPO(sites);
    for(int i = 1; i < N; ++i)
        {
        for (int j = 1; j < N; ++j)
            {
            ampo += (1/(N*N)),"Sz",i,"Sz",j;
            }
        }
    auto M2 = MPO(ampo);

Is there a reason why AutoMPO might not work in this context? Sorry if I’m missing something obvious here. Thanks a lot!

Best,
RS

---

Good question. So even with QNs AutoMPO ought to work just fine for this. Definitely it will work without QNs for an MPO of this type.

The only issue is a possibly mild one: for MPOs of this type, the AutoMPO back end will scale as N^2. So you won't be able to treat extremely large systems. Give it a try though and it may be just fine on the speed. In the future I have some ideas of how to make AutoMPO scale only as N for MPOs like this, but it's a work in progress.