Kronecker product and thermal states

Question

Hi All,
I am trying to implement a thermal state in the form @@exp(-\beta H/2)|I\rangle@@ where
$$
I = \sum_{k_1,k_2...} | k_1>|k_1> \otimes | k_2>|k_2> ...
$$
Here the main problem is to compute the single MPS @@ \sum_{k_i} | k_i>| k_i>@@ and then to take their Kronecker product. What would be the best way to compute such a state in itensor?

Thanks
Raffaele

Comments

I am just adding that all my sites are boson sites and the summation runs over a quite large number of basis states. I took a look at the Finite T code using itensor but I can't really see how to reuse it to use boson sites.
Thanks.

Answer 1

Hi Raffaele,
Good question: I think the tutorial/finiteT/ancilla.cc code you referenced would be a good starting point for doing what you are asking, but for other site types.

That code has a portion like this:

auto psi = MPS(sites);
for(int n = 1; n <= 2*N; n += 2)
    {
    auto s1 = sites(n);
    auto s2 = sites(n+1);
    auto wf = ITensor(s1,s2);
    wf.set(s1(1),s2(2), ISqrt2);
    wf.set(s1(2),s2(1), -ISqrt2);
    ITensor D;
    psi.ref(n) = ITensor(s1);
    psi.ref(n+1) = ITensor(s2);
    svd(wf,psi.ref(n),D,psi.ref(n+1));
    psi.ref(n) *= D;
    }

The main parts that would need to be generalized are:

1. instead of defining a singlet state |up,dn> - |dn,up> as in the above code, you'll probably want to define some other kind of maximally entangled state. As you probably know, the ancilla/purification method doesn't depend on the precise initial state used, as long as it's maximally entangled (unitarily equivalent to a Bell-pair state, meaning up to a unitary acting on the ancilla sites). A more general state that's good to define is one where the quantum numbers of the ancilla sites are different from the ones on the physical site in such a way that the overall entangled state has well-defined total quantum numbers. And you can just make all the coefficients the same sign. For the spin case above, the resulting state would be the Sz=0 triplet state rather than the singlet state.

2. you'll need to add two additional loops inside the code, where instead of writing each wf.set line by hand, you write something like:

   for(int n : range1(d))
   {
   wf.set(s1(n),s2(d+1-n), 1.0);
   }
   wf /= norm(wf); //normalize to 1.0

where above the d+1-n argument to s2 is because of the idea I mentioned in step #1 above about pairing physical-ancilla states together to make the total quantum number always the same (assuming here that your sites have the property that state n and state d+1-n together have a fixed number of particles).

Also note above d is the site dimension.

Even if d is very large, as you said is the case for you, the above code should run rather quickly. If it doesn't then you probably have a bug or we'll need to have a closer look.

Best,
Miles

Comments

Hi Miles,
really thanks for your detailed comment. I am trying the approach you suggested. You mentioned "two additional loop". What would the second loop be?
I mean if I understand correctly the code should be
    
    for(int n = 1; n <= 2*N; n += 2)
        {
        auto s1 = sites(n);
        auto s2 = sites(n+1);
        auto wf = ITensor(s1,s2);
        for(int k : range1(d))
        {
            wf.set(s1(k),s2(d+1-k), 1.0);
        }
        wf /= norm(wf); //normalize to 1.0
        ITensor D;
        psi.Aref(n) = ITensor(s1);
        psi.Aref(n+1) = ITensor(s2);
        svd(wf,psi.Aref(n),D,psi.Aref(n+1));
        psi.Aref(n) *= D;
        }

Is this correct or am I missing something?
Raffaele

---

Hi Raffaele, yes thanks for checking on this. I misspoke when I said "two more loops", because I was still thinking through the correct method and forgot to go back and edit my earlier text.


Your code above looks potentially correct. I say potentially because if you are conserving quantum numbers, you'll still need to check that the QN associated with s1(k) plus the QN associated with s2(d+1-k) always add up to the same total QN. That way the resulting wavefunction will have a well-defined total QN and you can use it in a QN conserving time evolution code. If you aren't conserving QNs, then your code above ought to work just fine.

Please try it out and check! It would be good if you separately made a code that produces the exact results for a small system to compare to. (E.g. for four physical sites you may be able to fully diagonalize H and construct exp(-beta*H) and trace it with operators etc.)

Best,
Miles