Overlaps with odd number of fermion operators

Question

I am using the Julia version of ITensor (v"0.2.4") to study fermionic systems where a degeneracy of the ground state is expected. As a simple example let us consider a doubly degenerate ground state with two wavefunctions say @@ |\Psi \rangle @@ and @@ |\Phi \rangle @@. Due to the nature of the problem, the wavefunctions differ by parity. In order to characterize the real space structure of the ground state I would like to calculate, among other things, overlaps of the following form:

$$ \langle \Psi | c^\dagger_i | \Phi \rangle $$

In order to calculate this overlap, I tried to define, using AutoMPO(), a MPO containing only a single local fermionic operator, which resulted in the following error:

"Parity-odd fermionic terms not yet supported by AutoMPO"

Is there a way around this? Is this possible in ITensor at all (either Julia or C++)?

Comments

Hi, this is a limitation of AutoMPO we would like to lift in the future. However, for an operator of this type, using a single ITensor to represent c^\dagger would be much more efficient than using a whole MPO. Would that solve your problem? We could discuss how to do this.

Miles

---

Miles, thanks for the replay! I would be glad if you could provide me with some pointers.

---

The simplest thing you could do is something like

    s = siteinds(phi)
    cdag_i = op("Cdag", s, i)
    cphi = apply(cdag_i, phi)
    inner(psi, cphi)

Basically all this does is create the operator `cdag_i` and then contract it onto the @@i@@th site of the state @@phi@@, then computes the overlap of the result with @@psi@@.

---

Thanks, Matt. The method Matt recommends is a good one. The other method would be more "manual" and would involve contracting all of the tensors of the two states together, like one is computing their dot product, except at some site "i" one puts the operator cdag_i (made as above) in between. Basically that is what the function "inner" does above, so it would be nearly a reimplementation of that but with an operator put in a certain place.

I would recommend Matt's method above!

---

Dear Matt and Miles!
Thanks, this is exactly what I needed! I was halfway down the rabbit hole already with looking into how I could extract the appropriate ITensors and figuring out how to manually build up the Cdag-s :) I guess, given my lack of detailed knowledge of ITensor, this would have been a painfully long albeit educative experience :)
Thanks again for this illuminating example!  I don't know what the usual practice here is but I think we can consider this question answered!

---

Glad to hear it’s now working for you!

---

By the way, we used exactly this technique to visualize the real space properties of Majorana edge states of topological 1D systems in this paper, just in case you’re interested to see (Fig. 3): https://arxiv.org/abs/1104.5493

---

Actually, the system I am looking at is something similar :). At some point, I also wanted to ask about that paper. I guess this is not the best place for it, but I was wondering about the actual dmrg parameters used like bond dimensions and precision goals for the svd. This, just as my problem, involves superconducting regions with spin-orbit coupling, so only parity remains from the possible fermionic symmetries. How far did you had to go with convergence parameters in order to get converged results?

---

Hm that’s a pretty detailed topic and the answer is “it depends”. Usually the idea is that you raise the maxdim as large as needed for each point in the phase diagram until the physics properties stop changing and the reported truncation error is less than at least 1e-5. Ideally it will be even smaller than that by a few orders of magnitude.

---

For 1D problems including this one, usually only bond dimensions of up to many hundreds are needed. But again it depends on things like the size of the gap and correlation length.

---

I think I roughly understand what to look for in order to ensure the convergence of a problem. I was wondering about the paper you mentioned above. I might have overlooked it but I did not find these convergency parameters in the paper.

---

Correct - I just checked and that information isn’t in the paper unfortunately. But the good thing is that DMRG tells you how accurate it is! So if you follow certain best practices (and always just do extra sweeps if your runs aren’t taking too long overall) and check that results look physical and smooth you should be ok. If you have more detailed questions about how to converge DMRG feel free to post a new question, though also you may find some discussion about this in previous questions.

---

I came across another issue regarding my original question. There might be still something amiss. I performed exact diagonalization on a small system and also calculated some low energy states with ITensor. It seems that the spectrum is the same for both calculations to 6-7 digits. However, the overlaps differ quite a bit. I realize that without divulging further details about my Hamiltonian it might be a bit hard for you to help me with debugging it. I will compile a short gist in the coming days and will share it here.

---

Do both codes conserve parity or just one of them?

---

In the exact diagonalization code, I have a sparse matrix representation of the Hamiltonian and just use scipy's eigsh, so in general, the solver does not know anything about parity. With ITensor I explicitly perform DMRG in the two parity sectors. I could however separate the exact diagonalization Hamiltonian into the two parity sectors to bring it closer to what happens in DMRG. I will try that and see if I get something closer to ITensor. Although I have to say the ED results look more symmetric, and look more physical to me.. this is just a hunch.

---

I think if the ED mixes parity sectors that could explain why the overlap is low while the energy is very close. For these systems it’s very important for the edge physics that the states have a well defined parity (eg note that c^dagger maps from one parity sector to the other, so to get a sensible result for your matrix element that you are computing it could be quite important to use parity conserving states for your two ground states).

---

The states seem to be a fixed parity in ED as well, that is @@\langle\Psi |c^\dagger_i|\Psi \rangle = \langle\Phi |c^\dagger_i|\Phi \rangle = 0 @@. My problem is that the spatial structure of  @@ \langle\Phi |c^\dagger_i|\Psi \rangle @@ is not the same in ED and DMRG.
Moreover, the ED structure is nicely symmetric on both edges while the DMRG is somewhat asymmetric. I will come back to you with a gist so you can see what I mean.

---

Just checked with explicit fermion parity preservation, and ED results stay the same, still different from DMRG.

---

I compiled a gist:
https://gist.github.com/oroszl/d113711a009512e398844d5aa8d06bd4
it has two notebooks, one Julia with the ITensor results and one Python with the ED results.

---

Hm it's a good question why the spatial structure is different, but unfortunately it's really hard to say without sitting down with the code and running a lot of different tests.

I would recommend thinking of any other things you can check: for example you can measure some more 'tame' properties too, just of each ground state individually, such as the density on each site. For this you can use the `expect` function in  ITensor to easily calculate this and any other local properties (expected values of local operators on each site). If these match the ED, then you have a mystery but it narrows down the possibilities. If these don't match ED, then very likely they just aren't the same wavefunctions despite having the same energies, which also constrains the logical possibilities. (E.g. although having a zero overlap with c^\dagger is consistent with having a well-defined parity, I'm not sure it proves that the parity is actually well-defined.)

But from the way you're (correctly) using our QN system, the two ground states should have well-defined parity. You can do `flux` on them again after DMRG runs to verify this.

---

Hi Miles,
thanks for your suggestions, I will go through them.

---

Hi guys!
I made some further benchmarks. Now I considered a simple 1D, 4 site Hubbard model with a little magnetic field.
- I performed ED calculations and DMRG calculations with fixed fermion numbers.
- The system according to ED has a non-degenerate ground state with 4 electrons, while the first excited state is doubly degenerate, one state has 3 while the other has 5 electrons.
- I calculated "tame" quantities as suggested by Miles, like charge distribution and magnetization. These quantities match DMRG results perfectly.
- I again calculated my odd fermion overlaps, and I got a puzzling result, namely, the quantities at the first site match perfectly between ED and DMRG, however, the following sites are wrong.
- I compiled again a gist:
  https://gist.github.com/oroszl/604aa7ee74e49bb408dbc031c03f3509

Based on this I have a hunch that this may be a "Jordan-Wigner string" issue, but I can not put my finger on it precisely.

---

Oh! Yes this was a big oversight of mine when answering your question. You are right: you have to explicitly include Jordan-Wigner string when using operators such as C^\dagger at the level of individual tensors in ITensor currently. (Other systems, like OpSum/AutoMPO include the string for you. This creates a somewhat confusing situation, but we are in the process of installing a much more 'automatic' fermion system - stay tuned.)

So then your code for measuring C^\dagger, if you are following Matt's recipe, should include extra lines like this:

    s = siteinds(phi)
    for j=1:i-1
      Fj = op("F",s,j)
     phi = apply(Fj,phi)
    end

    cdag_i = op("Cdag", s, i)
    cphi = apply(cdag_i, phi)
    inner(psi, cphi)

---

Right, I totally forgot about the Jordan-Wigner string. Thanks for the correction Miles.

---

Ok, this is Great!  Are there any extra precautions for spin-full systems?
Also, Miles in your code above does the line
apply(Fj,phi)  
automatically update phi in-place?  or it should rather be something like:
phi=apply(Fj,phi) # I am not sure what this does in Julia but it is legit in python

---

(1) good catch, yes it should be phi = apply(Fj,phi). I edited my comment above to reflect this.

(2) for spinful systems there can a subtlety depending on how the Jordan-Wigner mapping is defined. You can still just apply "F" on all of the sites to the left of site i. But then on site i the operator you need to apply can require on extra piece of string in certain cases i.e. when the Jordan-Wigner transformation maps completely to a bosonic operator obeying bosonic commutation relations.

*However* the operators "Cdagup", "Cdagdn", "Cup", and "Cdn" in the "Electron" site type in ITensor already include these extra, local minus signs that correspond to the "piece of string" I mentioned. So really the code above should 'just work' for spinful electrons too, as long as you just write "Cdagup" or "Cdagdn" for the operator whose matrix element you want and leave the rest of the code the same.

---

Great! Now Everything works like a charm!
Thank you guys for the quick help! ITensor was the best thing that happened to me in 2021!

---

Nice to hear that! Glad it's all working now.

Answer 1

(Please see discussion above.)