Tower Of States attractive Hubbard model

Author Hannes Karlsson

We perform a tower of states (TOS) analysis [1],[2] on the attractive Hubbard model on a square lattice. The (particle-hole symmetric) model can be written as $$ \mathcal{H} = -t\sum_{\langle i j \rangle} c_i^{\dagger} c_j - U\sum_{i} (n_{i\uparrow}-\frac{1}{2})(n_{i\downarrow}-\frac{1}{2}). $$

The attractive Hubbard model is known to display a superconducting ground state, which in mean-field theory is well described by a fermionic U\((1)\) coherent state $$ \ket{\psi} = \prod_{k} (u_k + v_k c_{k\uparrow}^\dagger c_{-k\downarrow}^\dagger)\ket{0}, $$ i.e, the ground state is degenerate in the number of particle pairs. We thus expect that in the superfluid state, in the neighbourhood of the most favourable particle number (determined by the chemical potential), there will be a tower of almost degenerate energy states, containing different number of superconducting pairs, but otherwise the same.

To perform the TOS analysis, we gather the eigenvalues using either full exact diagonalization (ED) or the Lanczos algorithm in each symmetry sector, defined by the particle number. Here we use the Lanczos algorithm with \(3\) eigenvalues well converged.

As we can see in the spectra above (computed for a \(4\times 4\) lattice), for large enough values of \(U\), the ground state energy is almost degenerate in the different particle number sectors, providing strong evidence for U\((1)\) spontaneous symmetry breaking.

C++Julia

#include <xdiag/all.hpp>

using namespace xdiag;

int main() try {
  say_hello();
  int Lx = 4;
  int Ly = 4;
  int Nmin = 4;
  int Nmax = 2 * Lx * Ly - 4;
  std::vector<double> Us = {0.0, -2.0, -10.0};

  for (double U : Us) {

    // Construct Hamiltonian
    auto ops = OpSum();

    ops += "U" * Op("HubbardU"); // apply Hubbard interaction to all sites

    for (int x = 0; x < Lx; x++) {
      for (int y = 0; y < Ly; y++) {
        int s = x * Ly + y;
        int sx = ((x + 1) % Lx) * Ly + y;
        int sy = x * Ly + (y + 1) % Ly;

        ops += "t" * Op("Hop", {s, sx}); // hopping in x-direction
        ops += "t" * Op("Hop", {s, sy}); // hopping in y-direction
        ops += "mu" * Op("Ntot", s);     // chemical potential
      }
    }
    ops["t"] = 1;
    ops["mu"] = -U / 2;
    ops["U"] = U;

    for (int N = Nmin; N <= Nmax; N += 2) {
      int nup = N / 2;
      int ndn = N / 2;
      auto block = Electron(Lx * Ly, nup, ndn); // create Hilbert space

      auto res = eigs_lanczos(ops, block, 3);

      std::string outfile = fmt::format(
          "data/tos_ahm/U({})_N({})_Lx({})_Ly({}).h5", U, N, Lx, Ly);
      auto f = FileH5(outfile, "w!");
      f["spectrum"] = res.eigenvalues;
    }
  }
  return 0;
} catch (Error e) {
  error_trace(e);
}

using XDiag
using HDF5

"""Struct for lattice sites"""
struct LatticeBond
    s1::Int64
    x1::Float64
    y1::Float64
    s2::Int64
    x2::Float64
    y2::Float64
end

"""Returns list of all bonds on the square lattice"""
function square_lattice(Lx::Int,Ly::Int;xperiodic=true,yperiodic=true)

    bonds = Vector{LatticeBond}()

    for x1 in 1:Lx, y1 in 1:Ly
        s1 = (x1-1)*Ly+y1
        for (x2,y2) in [( mod(x1,Lx)+1, y1 ),( x1, mod(y1,Ly)+1 )]
            s2 = (x2-1)*Ly+y2
            set = reduce(hcat,[[s1,x1,y1],[s2,x2,y2]])
            set = set[:,sortperm(set[1,:])]     # sort according to s1 and s2
            b = LatticeBond(set[:,1]...,set[:,2]...)
            push!(bonds,b)
        end
    end
    return bonds
end

"""Returns the Hamiltonian for the Hubbard model"""
function ham(U::Float64,Lx::Int64,Ly::Int64)

    lattice = square_lattice(Lx,Ly)
    ops = OpSum()

    for b in lattice # hopping term
        ops += "t" * Op("Hop", [b.s1,b.s2])     # both for spin up and down
    end

    ops += "U" * Op("HubbardU")     # applies Hubbard interaction over the entire lattice
    for i in 1:Lx*Ly
        ops += "mu" * Op("Ntot",i)     # chemical potential
    end

    ops["t"] = 1.0
    ops["U"] = U
    ops["mu"] = -U/2

    return ops

end

function main(U::Float64,Lx::Int64,Ly::Int64,Ns::StepRange{Int64, Int64})

    for N in Ns

        ops = ham(U,Lx,Ly) # create the Hamiltonian
        nup = N%2 == 0 ? N ÷ 2 : (N+1) ÷ 2
        ndn = N%2 == 0 ? N ÷ 2 : (N-1) ÷ 2

        block = Electron(Lx*Ly, nup, ndn)   # create Hilbert space
        res = eigs_lanczos(ops, block, neigvals = 3)
        eigs = res.eigenvalues

        filename = "data/tos_ahm/U($U)_N($N)_Lx($Lx)_Ly($Ly).h5"

        h5open(filename,"w") do f
            write(f,"spectrum",eigs)
        end

    end

end

Us = [0.0,-2.0,-10.0]
Lx, Ly = (4,4)
Nmin = 4
Nmax = 2*Lx*Ly - 4
Ns = Nmin:2:Nmax

say_hello()
for U in Us
    println("Now doing U=$U")
    main(U,Lx,Ly,Ns)
end

references

[1] P. W. Anderson, An Approximate Quantum Theory of the Antiferromagnetic Ground State, Phys. Rev. 86, 694 (1952)

[3] Alexander Wietek, Michael Schuler and Andreas M. Läuchli, Studying Continuous Symmetry Breaking using Energy Level Spectroscopy, arXiv:1704.08622 (2017).