General XXZ Heisenberg Hamiltonian

[RichRick1]

Before coding the Heisenberg Hamiltonian it's good to rewrite it in terms of spin-creation/annihilation operators.
We start from the form that is available in notes:

$$\hat{H}{X X Z}=\sum_p\left(\mu_p^Z-J{p p}^{\mathrm{eq}}\right) S_p^Z+\sum_{p q} J_{p q}^{\mathrm{ax}} S_p^Z S_q^Z+\sum_{p q} J_{p q}^{\mathrm{eq}} S_p^{+} S_q^{-}$$

Using the relationships:
$S_p^{+}=a_{p \alpha}^{\dagger} a_{p \beta}^{\dagger}$
$S_p^{-}=a_{p \beta} a_{p \alpha}$
$S_p^Z=\frac{1}{2}\left(a_{p \alpha}^{\dagger} a_{p \alpha}+a_{p \beta}^{\dagger} a_{p \beta}-1\right)$

We can get:
$H_{x x z}=\sum_p\left(\mu_p^z-J_{p p}^{eq}\right) S_p^z+\sum_{p q} J_{p q}^{a x} S_p^z S_q^z+\sum_{p q} J_{p q}^{e q} S_p^{+} S_q^{-} =$
$=\sum_p\left(\mu_p^z-J_{ pp}^{e q}\right)\cdot \frac{1}{2}\left(\hat{n}_p-1\right) +$ $\sum _{pq} J _{p q}^{a x} \cdot \frac{1}{4}\left(\hat{n}_p-1\right) \left(\hat{n}_q-1\right)$ + $\sum _{p q} J _{p q}^{e q} a _{p \alpha}^{+} a _{p \beta}^{+} a _{q \beta} a _{q \alpha}$
Then:
Zero body term
$h_0 = -\frac{1}{2} \sum _p\left(\mu _p^z-J _{p p}^{e q}\right)+\frac{1}{4} \sum _{p q} J _{p q}^{a x}$

One body term
$h_1 = \frac{1}{2} \sum _p\left(\mu _p^z- J _{p p}^{e q}\right) \cdot \hat{n}_p-\frac{1}{4} \sum _{p q} J _{p q}^{a x}\left(\hat{n}_p+n_q\right)$ =
$= \frac{1}{2} \sum _p\left(\mu _p^z-J _{p p}^{e q}\right) \tilde{n}_p-\frac{1}{4} \sum _{p q} J _{p q}^{a x} \hat{n}_p -\frac{1}{4} \sum _{p q} J _{p q}^{a x} \hat{n}_q =$
$=\frac{1}{2} \sum _p\left(\mu_p^z-J _{p p}^{e q}\right) \hat{n}_p-\frac{1}{4} \sum _{p q}\left(J _{p q}^{a x}+J _{q p}^{a x}\right) \hat{n}_p =$
$=\frac{1}{2} \sum _p\left(\mu_p^z-J _{p p}^{e q}\right) \tilde{n}_p-\frac{1}{2} \sum _{pq} J _{p q}^{a x} \hat{n}_p=$
$=\frac{1}{2} \sum _p\left(\mu_p^z-J _{p p}^{e q}-\sum _{q} J _{p q}^{a x}\right) \hat{n}_p$

Two body term
$h_2=\sum _{p q} J _{p q}^{eq} a _{p \alpha}^{+} a _{p \beta}^{+} a _{q \beta} a _{q \alpha}+\frac{1}{4} \sum _{p q} J _{p q} ^{ax} \hat{n}_p \hat{n}_q$
Considering that
$\hat{n} _p \hat{n} _q=\left(\hat{n} _{p\alpha}+\hat{n} _{p \beta}\right) \left(\hat{n} _{q \alpha}+\hat{n} _{q \beta}\right)=$
$= \hat{n} _{p \alpha} \hat{n} _{q \alpha}+ \hat{n} _{p \alpha} \hat{n} _{q \beta}+\hat{n} _{p \beta} \hat{n} _{q \alpha}+\hat{n} _{p \beta} \hat{n} _{q \beta} =$
$=a _{p \alpha}^{+} a _{p \alpha} a _{q \alpha}^{+} a _{q \alpha} + a _{p\alpha}^{+} a _{p\alpha} a _{q \beta}^{+} a _{q\beta}+ a _{p \beta}^{+} a _{p \beta} a _{q \alpha}^{+} a _{q \alpha}+ a _{p \beta}^{+} a _{p \beta} a _{q \beta}^{+} a _{q \beta}=$
$= a _{p \alpha}^{+} a _{q\alpha}^{+} a _{p \alpha} a _{q\alpha} + a _{p \alpha}^{+} a _{q \beta}^{+} a _{p\alpha} a _{q \beta} + a _{p \beta}^{+} a _{q \alpha}^{+} a _{p \beta} a _{q \alpha}+a _{p \beta}^{+} a _{q \beta}^{+} a _{p \beta} a _{p \beta}$

So we can rewrite links between matrices:
$h _{pp} ^{(1)}= \frac{1}{2}\left(\mu_p^z-J _{p p}^{e q}-\sum _{q} J _{p q}^{a x}\right)$,
$h _{pq} ^{(1)} = 0$ if $p \neq q$

$V _{p q p q} = +\frac{1}{4} J _{p q}^{e q}$
$V _{p(q+n) p (q+n)} = +\frac{1}{4} J _{p q}^{ax}$
$V _{p+n) q (p+n) q} = +\frac{1}{4} J _{p q}^{ax}$
$V _{(p+n)(q+n)(p+n)(q+n)} = +\frac{1}{4} J _{p q}^{ax}$
$V _{p(p+n)(q + n) q} = J _{p q}^{e q}$

Also, we can make a mapping to the PPP+ Hamiltonian:

$$\begin{aligned} & h_{p p}=0, p \neq q \\ & h_{p p}=\left(\mu_p^z-J_{p p}^{e q}-\sum_{q} J_{p q}^{a x}\right) \cdot \frac{1}{2} \\ & Q_p=Q_q=1 \\ & \gamma_{p q}=\frac{1}{2} J_{p q}^{a x} \\ & g_{p q}=J_{p q}^{\text {eq }} \end{aligned}$$

[RichRick1]

@gabrielasd can you please check if the math is correct.
@PaulWAyers, I'm not quite sure if the double sum should be so that $p \neq q$

[gabrielasd]

I agree with Rick's derivation above taking p!=q for all double summations.

[PaulWAyers]

The double-sum $p=q$ term can be included. For the $S_p^Z S_p^Z$ term, ($J_{pp}^{ax}$) this "almost" simplifies to a 1-body term, but not quite since there is a term that is a product of the $n_{p\alpha} and $n_{p\beta}$ occupations. That's an important term in the 2DM which corresponds to an "on site" interaction.

The $S_p^+ S_p^-$ term is basically equivalent (up to one-body terms) but sometimes the Hamiltonian is written using the $+-$ term, and sometimes using the $ZZ$ term. It's good to support both, as the user may like to use one rather than the other. In an Ising-type model, only the first one would be relevant.

The identity

$$ S_p^+ S_p^- = S_p^2 - S_p^Z S_p^Z - S_p^Z $$

basically explains why these two are equivalent.

[RichRick1]

@PaulWAyers @gabrielasd I think I made a mistake at one body term
It should be : $h _{pp} ^{(1)}= \frac{1}{2}\left(\mu_p^z-J _{p p}^{e q}-\sum _{q} J _{p q}^{a x}\right)$
I've just updated

[PaulWAyers]

I didn't double-check all the math. I think it should be clear from running a few tests to be sure FCI for Heisenberg reference data is reproduced here. I didn't see any errors though.

[RichRick1]

According to the definitions of Heisenberg model I saw in couple papers (like this one) it's common to define the HHZ Heisenberg model using neighbourhood interactions only. So I propose to give user an option:

1. provide float values for J_eq and J_ax and our code will build the hamiltonian based on connectivity matrix based on only neighbourhood interactions (Huckel-like approach)
2. provide exact matrices for J_eq and J_ax and then no connectivity is needed, because all the necessary information about interactions is already stored at the matrices.

@PaulWAyers @gabrielasd what do you think?

[PaulWAyers]

I agree. Options 1 and 2 cover the main cases. I can imagine a few special cases but they are best implemented as routines for generating connectivity matrices.

[RichRick1]

Want to clarify on thing. According to the math above one body term has only diagonal component, however in Ising model we have interactions of the closest site so it should impact the off-diagonal terms. Am I right and we are missing something in the definition or everything is ok? @PaulWAyers @gabrielasd
I think the main difference is that we are using SU2 formalism while the most common approach is Jordan-Wigner mapping, but should this affect the resulting energies? @PaulWAyers

[PaulWAyers]

The Ising model uses only occupation numbers, so the two-site term is an 2-body term $\hat{n}_i \hat{n}_j$. So I think it is OK.

The Jordan-Wigner mapping (or similar, like Bravyi-Kitaev) is more common, but the mapping from fermions to spin requires very high-order spin terms, of the order of $O(N_{basis})$ or $O(N_{basis}^2)$.