SyTen
syten-sql-mps-fermi-hubbard-kns.cpp File Reference

Creates a $${\mathbb{Z}_W} \times \mathrm{U}(1) \times \mathrm{SU}(2)$$ (momentum × charge × spin)-symmetric or a $${\mathbb{Z}_W} \times \mathrm{U}(1) \times \mathrm{U}(1)$$ (momentum × charge × z-spin)-symmetric 2D real/momentum space Fermi-Hubbard model. More...

#include "inc/util/bpo_helper.h"
#include "inc/util/output.h"
#include "inc/util/exit_codes.h"
#include "inc/util/constants.h"
#include "inc/dense/dense_operators.h"
#include "lat/mps/sql-helpers.h"
#include "lat/mps/kspace-helpers.h"
#include "lat/mps/reordering-helpers.h"
#include "lat/mps/su2u1z.h"
#include "lat/mps/u1u1z.h"
Include dependency graph for syten-sql-mps-fermi-hubbard-kns.cpp:

## Namespaces

syten
Syten namespace.

## Functions

MPS::Operator syten::build_hopping (MPS::Lattice const &lat, std::string const &sym, Index const first, Index const second)
Returns the MPO rep of the (non-hermitian/‘half’) hopping term c^\dagger_first \cdot c_second. More...

void syten::create_dg_kinetic_term (MPS::Lattice &lat, std::string const &sym, Index const length, Index const width, bool const bpo_quiet, std::ostream &out, std::function< Index(int, int)> coord, std::function< Index(int)> kfun, std::vector< Index > idx)
Creates a kinetic term in the diagonal lattice lat. More...

void syten::create_interaction_term (MPS::Lattice &lat, std::string const &sym, Index const length, Index const width, bool const bpo_quiet, std::ostream &out, std::function< Index(int, int)> coord, std::function< Index(int)> kfun, std::vector< Index > idx)
Creates the interaction term in the lattice lat. More...

void syten::create_ring_repulsion (MPS::Lattice &lat, std::string const &sym, Index const length, Index const width, bool const bpo_quiet, std::ostream &out, std::function< Index(int, int)> coord, std::function< Index(int)> kfun, std::vector< Index > idx)
Creates a ring-wise quadratic interaction term which favours equidistribution between different rings. More...

void syten::create_rs_occ_operators (MPS::Lattice &lat, std::string const &sym, Index const length, Index const width, bool const bpo_quiet, std::ostream &out, std::function< Index(int, int)> coord, std::function< Index(int)> kfun, std::vector< Index > idx)
For every point of the lattice, creates a real-space occupation number operator to facilitate the calculation of real-space expectation values. More...

void syten::create_rs_spin_operators (MPS::Lattice &lat, std::string const &sym, Index const length, Index const width, bool const bpo_quiet, std::ostream &out, std::function< Index(int, int)> coord, std::function< Index(int)> kfun, std::vector< Index > idx)
For every point of the lattice, creates a real-space spin operator to facilitate the calculation of real-space expectation values. More...

void syten::create_sq_inter_ring_kinetic_term (MPS::Lattice &lat, std::string const &sym, Index const length, Index const width, bool const bpo_quiet, std::ostream &out, std::function< Index(int, int)> coord, std::function< Index(int)> kfun, std::vector< Index > idx)
Creates a length-wise kinetic term in the square lattice lat. More...

void syten::create_sq_intra_ring_kinetic_term (MPS::Lattice &lat, std::string const &sym, Index const length, Index const width, bool const bpo_quiet, std::ostream &out, std::function< Index(int, int)> coord, std::function< Index(int)> kfun, std::vector< Index > idx)
Creates a width-wise kinetic term in the square lattice lat. More...

## Detailed Description

Creates a $${\mathbb{Z}_W} \times \mathrm{U}(1) \times \mathrm{SU}(2)$$ (momentum × charge × spin)-symmetric or a $${\mathbb{Z}_W} \times \mathrm{U}(1) \times \mathrm{U}(1)$$ (momentum × charge × z-spin)-symmetric 2D real/momentum space Fermi-Hubbard model.

## Square Lattice (PBC/OBC parallel to lattice vectors)

The model lives on a cylinder of length $$L$$ and width/circumference $$W$$. Usually, $$W < L$$.

We then can have three terms in the Hamiltonian $$\hat H = t_W \hat H_{t_W} + t_L \hat H_{t_L} + U \hat H_U$$. $$t_W$$ is the hopping coefficient within a single ring of the cylinder; $$t_L$$ is the inter-ring hopping coefficient and $$U$$ is the on-site interaction.

Greek indices are used for length-wise coordinates, $$j$$ is used for width-wise real-space coordinates and $$k, l, m, n$$ are width-wise momentum-space coordinates.

Specifically, we have (with the identification that $$j = W \equiv j = 0$$):

\begin{align*} H_{t_W} & = -\sum_{\alpha=1}^L \sum_{j=1}^W c^\dagger_{\alpha,j} \cdot c_{\alpha,(j+1)} + c^\dagger_{\alpha,(j+1)} \cdot c_{\alpha, j} \\ H_{t_L} & = -\sum_{\alpha=1}^{L-1} \sum_{j=1}^W c^\dagger_{\alpha,j} \cdot c_{\alpha+1,j} + c^\dagger_{\alpha+1,j} \cdot c_{\alpha, j} \\ H_{U} & = \frac{1}{2} \sum_{\alpha=1}^L \sum_{j=1}^W \left( c^\dagger_{\alpha,j} \cdot c_{\alpha,j} \times c^\dagger_{\alpha,j} \cdot c_{\alpha,j} - c^\dagger_{\alpha,j} \cdot c_{\alpha,j} \right) \end{align*}

In the above and all the following $$A^\dagger \cdot B$$ denotes the dot product between $$A$$ and $$B$$ (i.e. the code uses dot(B, A)). $$A \times B$$ is a product between trivially-transforming operators and has lower precendence than the dot product.

Note that the last term is essentially $$n_i^2 - n_i$$ with the full number operator $$n_i \in \{ 0, 1, 2 \}$$. The usual term $$n_{i,\uparrow} n_{i,\downarrow}$$ is not $$\mathrm{SU}(2)$$-invariant.

We can then do a Fourier transform along the rings of the cylinder as

\begin{align*} c^\dagger_{\alpha,j} & = \frac{1}{\sqrt{W}} \sum_{k=1}^W e^{-i k j} c^\dagger_{\alpha,k} \\ c_{\alpha,j} & = \frac{1}{\sqrt{W}} \sum_{k=1}^W e^{i k j} c_{\alpha,k} \end{align*}

Note the different signs in the exponent. The numeric prefactor of $$k$$ etc. can be fixed later. For the first term $$H_{t_W}$$, we get (with $$\frac{1}{N} \sum_j e^{ijk} = \delta_{k,0}$$):

\begin{align*} & -\sum_{\alpha=1}^L \sum_{j=1}^W \left( \frac{1}{W} \sum_{kl=1}^W e^{-ikj} e^{il(j+1)} c^\dagger_{\alpha,k} \cdot c_{\alpha,l} + \sum_{kl=1}^W e^{-ik(j+1)} e^{ilj} c^\dagger_{\alpha,k} \cdot c_{\alpha,l} \right) \\ = & -\sum_{\alpha=1}^L \frac{1}{W} \sum_{kl=1}^W \sum_{j=1}^W \left( e^{ij(l-k)} e^{il} + e^{ij(l-k)} e^{-ik} \right) c^\dagger_{\alpha,k} \cdot c_{\alpha,l} \\ = & -\sum_{\alpha=1}^L \sum_{k=1}^W \sum_{l=1}^W \left( \delta_{l,k} e^{il} + \delta_{l,k} e^{-ik} \right) c^\dagger_{\alpha,k} \cdot c_{\alpha,l} \\ = & -\sum_{\alpha=1}^L \sum_{k=1}^W \left( e^{ik} + e^{-ik} \right) c^\dagger_{\alpha,k} \cdot c_{\alpha,k} \\ = & -\sum_{\alpha=1}^L \sum_{k=1}^W 2 \cos\left(2 \pi \frac{k}{W}\right) c^\dagger_{\alpha,k} \cdot c_{\alpha,k} \end{align*}

In the last line, we snuck in the proper prefactor for $$k$$. The second term gives us

\begin{align*} & -\sum_{\alpha=1}^{L-1} \sum_{j=1}^W c^\dagger_{\alpha,j} \cdot c_{\alpha+1,j} + c^\dagger_{\alpha+1,j} \cdot c_{\alpha, j} \\ = & -\sum_{\alpha=1}^{L-1} \sum_{j=1}^W \frac{1}{N} \left( \sum_k \sum_l \left[ e^{-ikj} e^{ilj} c^\dagger_{\alpha,k} \cdot c_{\alpha+1,l} \right] + \sum_k \sum_l \left[ e^{-ikj} e^{ilj} c^\dagger_{\alpha+1,k} \cdot c_{\alpha,l} \right] \right) \\ = & -\sum_{\alpha=1}^{L-1} \frac{1}{N} \left( \sum_k \sum_l \left[ N \delta_{k,l} c^\dagger_{\alpha,k} \cdot c_{\alpha+1,l} \right] + \sum_k \sum_l \left[ N \delta_{k,l} c^\dagger_{\alpha+1,k} \cdot c_{\alpha,l} \right] \right) \\ = & -\sum_{\alpha=1}^{L-1} \sum_{k=1}^W \left( c^\dagger_{\alpha,k} \cdot c_{\alpha+1,k} + c^\dagger_{\alpha,k} \cdot c_{\alpha+1,k} \right) \end{align*}

These are somewhat expected: electrons with momentum $$k$$ along the ring gain energy according to the usual dispersion relation. Electrons hopping from one ring to the other keep their momentum in the orthogonal ring direction.

The interaction term is rather annoying:

\begin{align*} & \frac{1}{2} \sum_{\alpha=1}^L \sum_{j=1}^W \left( \left[ \sum_{klmn} \frac{1}{W^2} e^{-ikj} e^{ilj} e^{-imj} e^{inj} c^\dagger_{\alpha,k} \cdot c_{\alpha,l} \times c^\dagger_{\alpha,m} \cdot c_{\alpha,n} \right] - \left[ \sum_{pq} \frac{1}{W} e^{-ipj} e^{iqj} c^\dagger_{\alpha,p} \cdot c_{\alpha,q} \right] \right) \\ = & \frac{1}{2} \sum_{\alpha=1}^L \left( \left[ \sum_{klmn} \delta_{k-l+m,n} c^\dagger_{\alpha,k} \cdot c_{\alpha,l} \times c^\dagger_{\alpha,m} \cdot c_{\alpha,n} \right]- \left[ \sum_{pq} \delta_{p,q} c^\dagger_{\alpha,p} \cdot c_{\alpha,q} \right] \right) \\ = & \frac{1}{2} \sum_{\alpha=1}^L \left( \left[ \sum_{klm} \frac{1}{W} c^\dagger_{\alpha,k} \cdot c_{\alpha,l} \times c^\dagger_{\alpha,m} \cdot c_{\alpha,k-l+m} \right]- \left[ \sum_p c^\dagger_{\alpha,p} \cdot c_{\alpha,p}\right] \right) \\ = & \frac{1}{2} \sum_{\alpha=1}^L \sum_{k=1}^W \left( \left[ \sum_{lm=1}^W \frac{1}{W} c^\dagger_{\alpha,k} \cdot c_{\alpha,l} \times c^\dagger_{\alpha,m} \cdot c_{\alpha,k-l+m} \right] - c^\dagger_{\alpha,k} \cdot c_{\alpha,k} \right) \end{align*}

Note that here, if we expand the scalar products fully again, we get the following expression:

\begin{align*} & c^\dagger_{\alpha,k} \cdot c_{\alpha,l} \times c^\dagger_{\alpha,m} c_{\alpha,k+m-l} \\ = & c^\dagger_{\alpha,k,\uparrow} c_{\alpha,l,\uparrow} c^\dagger_{\alpha,m,\uparrow} c_{\alpha,k+m-l,\uparrow} + c^\dagger_{\alpha,k,\uparrow} c_{\alpha,l,\uparrow} c^\dagger_{\alpha,m,\downarrow} c_{\alpha,k+m-l,\downarrow} + c^\dagger_{\alpha,k,\downarrow} c_{\alpha,l,\downarrow} c^\dagger_{\alpha,m,\uparrow} c_{\alpha,k+m-l,\uparrow} + c^\dagger_{\alpha,k,\downarrow} c_{\alpha,l,\downarrow} c^\dagger_{\alpha,m,\downarrow} c_{\alpha,k+m-l,\downarrow} \end{align*}

The mixed terms are exactly the ones we want. We can immediately commute the spin-up terms to the front of the second term and then by the following relabelling:

\begin{align*} & c^\dagger_{\alpha,k,\downarrow} c_{\alpha,l,\downarrow} c^\dagger_{\alpha,m,\uparrow} c_{\alpha,k+m-l,\uparrow} & \\ = \quad & c^\dagger_{\alpha,m,\uparrow} c_{\alpha,k+m-l,\uparrow} c^\dagger_{\alpha,k,\downarrow} c_{\alpha,l,\downarrow} & m \to a, k \to b, l \to c \\ = \quad & c^\dagger_{\alpha,a,\uparrow} c_{\alpha,a+b-c,\uparrow} c^\dagger_{\alpha,b,\downarrow} c_{\alpha,c,\downarrow} & b+a-c = d \Rightarrow c = a+b-d \\ = \quad & c^\dagger_{\alpha,a,\uparrow} c_{\alpha,d,\uparrow} c^\dagger_{\alpha,b,\downarrow} c_{\alpha,a+b-d,\downarrow} \\ = \quad & c^\dagger_{\alpha,k,\uparrow} c_{\alpha,l,\uparrow} c^\dagger_{\alpha,m,\downarrow} c_{\alpha,k+m-l,\downarrow} \end{align*}

We get twice the term $$c^\dagger_{\alpha,k,\uparrow} c_{\alpha,l,\uparrow} c^\dagger_{\alpha,m,\downarrow} c_{\alpha,k+m-l,\downarrow}$$ as required from a Fourier transformation of the original $$n_\uparrow n_\downarrow$$.

The terms with only spin-up or only spin-down components are removed as follows: By commuting the creators to the front, we get:

$c^\dagger_{\alpha,k,\uparrow} \left( \delta_{m,l} - c^\dagger_{\alpha,m,\uparrow} c_{\alpha,l,\uparrow} \right) c_{\alpha,k+m-l,\uparrow} + c^\dagger_{\alpha,k,\downarrow} \left( \delta_{m,l} - c^\dagger_{\alpha,m,\downarrow} c_{\alpha,l,\downarrow} \right) c_{\alpha,k+m-l,\downarrow}$

Note that the $$\delta_{m,l}$$ terms cancel exactly with the subtracted terms resulting from $$c^\dagger_{\alpha,k} \cdot c_{\alpha,k}$$. Then, we are left with the terms $$c^\dagger_{\alpha,k,\sigma} c^\dagger_{\alpha,m,\sigma}$$ and sums over both $$k$$ and $$m$$. For $$k = m$$, the term is zero (two creators). For $$k \neq m$$, the term cancels with the corresponding $$m \neq k$$ term (symmetric sum over antisymmetric term; thanks Leo!).

Hence, in total we are left exactly with

$\frac{1}{W} \sum_{klm=0}^W c^\dagger_{\alpha,k,\uparrow} c_{\alpha,l,\uparrow} c^\dagger_{\alpha,m,\downarrow} c_{\alpha,k+m-l,\downarrow}$

## Diagonal Lattice (PBC/OBC at 45° to lattice vectors)

For rungs α and legs j of the ladder, the kinetic Hamiltonian in real-space is given as

\begin{align*} H_t = & - \sum_{\alpha} \sum_{j} \left( c^\dagger_{\alpha, j} \cdot c_{\alpha+1, j} + c^\dagger_{\alpha+1, j} \cdot c_{\alpha, j} \right) \\ & - \sum_{\textrm{even }\alpha} \sum_{j} \left( c^\dagger_{\alpha, j} \cdot c_{\alpha+1, j+1} + c^\dagger_{\alpha+1, j+1} \cdot c_{\alpha, j} \right) \\ & - \sum_{\textrm{odd }\alpha} \sum_{j} \left( c^\dagger_{\alpha, j} \cdot c_{\alpha+1, j-1} + c^\dagger_{\alpha+1, j-1} \cdot c_{\alpha, j} \right) \end{align*}

The interaction term is as before

$H_U = \frac{1}{2} \sum_{\alpha} \sum_{j} c^\dagger_{\alpha,j} \cdot c_{\alpha, j} \times c^\dagger_{\alpha,j} \cdot c_{\alpha, j} - c^\dagger_{\alpha,j} \cdot c_{\alpha, j}$

with the total Hamiltonian $$H = t H_t + U H_U$$.

The transformations are the same as before along the rungs of the ladder:

\begin{align*} c^\dagger_{\alpha,j} & = \frac{1}{\sqrt{W}} \sum_k e^{-ikj} c^\dagger_{\alpha,k} \\ c_{\alpha,j} & = \frac{1}{\sqrt{W}} \sum_k e^{ikj} c_{\alpha,k} \\ \delta_{k,0} & = \frac{1}{W} \sum_j e^{ikj} \end{align*}

The kinetic part of the Hamiltonian is then transformed as follows:

\begin{align*} H_t = & - \frac{1}{W} \sum_{\alpha} \sum_{j} \sum_k \sum_l \left( c^\dagger_{\alpha, k} \cdot c_{\alpha+1, l} e^{i(l-k)j} + c^\dagger_{\alpha+1, k} \cdot c_{\alpha, l} e^{i(l-k)j} \right) \\ & - \frac{1}{W} \sum_{\textrm{even }\alpha} \sum_{j} \sum_k \sum_l \left( c^\dagger_{\alpha, k} \cdot c_{\alpha+1, l} e^{i(l-k)j} e^{il} + c^\dagger_{\alpha+1, k} \cdot c_{\alpha, l} e^{i(l-k)j} e^{-ik} \right) \\ & - \frac{1}{W} \sum_{\textrm{odd }\alpha} \sum_{j} \sum_k \sum_l \left( c^\dagger_{\alpha, k} \cdot c_{\alpha+1, l} e^{i(l-k)j} e^{-il} + c^\dagger_{\alpha+1, k} \cdot c_{\alpha, l} e^{i(l-k)j} e^{ik} \right) \\ = & - \sum_{\alpha} \sum_k \left( c^\dagger_{\alpha, k} \cdot c_{\alpha+1, k} + c^\dagger_{\alpha+1, k} \cdot c_{\alpha, k} \right) \\ & - \sum_{\textrm{even }\alpha} \sum_k \left( c^\dagger_{\alpha, k} \cdot c_{\alpha+1, k} e^{ik} + c^\dagger_{\alpha+1, k} \cdot c_{\alpha, k} e^{-ik} \right) \\ & - \sum_{\textrm{odd }\alpha} \sum_k \left( c^\dagger_{\alpha, k} \cdot c_{\alpha+1, k} e^{-ik} + c^\dagger_{\alpha+1, k} \cdot c_{\alpha, k} e^{ik} \right) \\ = & - \sum_{\alpha} \sum_k \left( c^\dagger_{\alpha, k} \cdot c_{\alpha+1,k} \left(1 + e^{(-1)^\alpha ik 2 \pi / W}\right) + c^\dagger_{\alpha+1,k} \cdot c_{\alpha,k} \left(1 + e^{(-1)^{\alpha+1} ik 2 \pi / W} \right) \right) \end{align*}

The shape of the interaction term does not change but stays the same as before:

$H_U = \frac{1}{2} \sum_\alpha \sum_k \left( \sum_l \sum_m \left[ \frac{1}{W} c^\dagger_{\alpha,k} \cdot c_{\alpha,l} \times c^\dagger_{\alpha,m} \cdot c_{\alpha, k-l+m} \right] - c^\dagger_{\alpha,k} \cdot c_{\alpha,k} \right)$

## Real-Space Spin and Particle Number Operators (both lattices)

Given a real-space site $$j$$, we can write the spin operator as $$\hat S^a_j = c^\dagger_j \sigma^a c_j \sqrt{\frac{3}{4}}$$ where $$\sigma^a$$ is a vector of the Pauli matrices. In code, this is just the dot product of c with itself in the S=1, n=0 irrep.

With the usual Fourier transform, we arrive at

$S^a_{\alpha,j} = \sqrt{\frac{3}{4}} \frac{1}{W} \sum_k \sum_l e^{\frac{2\pi i}{W} j (l-k)} c^\dagger_{\alpha,k} \sigma^a c_{\alpha,l}$

Note that this operator does not transform uniquely in k-space, as the mixed terms can contribute arbitrary k-values. However, the dot product of this operator with another one $$S^a_{\beta, i}$$ will give the correct result.

Similarly, we can write the particle number operator as the dot product of c with itself in the S=0, n=0 irrep. Again, this operator will not transform uniquely wrt k-values.

## (Re)ordering

There are a number of transformations done during the creation of this lattice.

We start with two-dimensional hybrid space coordinates $$(x, k)$$. The momentum ordering function (--momentum) or the reordering specification (--reordering) decide how $$k$$ is mapped to a logical position $$y = \mathrm{kfun}(k)$$ on a 2D lattice.

Then, given the 2D lattice with coordinates $$(x, y = \mathrm{kfun}(k))$$, the 2D→1D map (--map) specifies how the one-dimensional MPS ‘snake’ is placed through this lattice, resulting in a 1D coordinate $$i = \mathrm{coord}(x, y = \mathrm{kfun}(k))$$.

Finally, the full reordering (--full-reordering) takes the 1D coordinate $$i$$ and maps it to another 1D coordinate $$j = \mathrm{fr}(\mathrm{coord}(x, y = \mathrm{kfun}(k)))$$ which is then used to build up the MPS chain.

This is useful since we can independently specify the ordering of momentum sites within each rung in the first step, then specifiy how we want to map the 2D lattice into a 1D chain and finally apply arbitrary transformations as they might arise from Dynamically-Reordering DMRG.