Filling for Bosonic Green's Functions

[andrewkhardy]

I am interested in constructing a bosonic Green's function in TRIQS (a holon to be precise) and finding the chemical potential to fix the filling at some number $n_B < 1$ for a finite temperature. I can compute the density for a corresponding Green's function using

$$n_b = -G_B(\tau = 0^+) - 1$$.

The definition for the bosonic-holon Green's function I am using is

$$G(i \nu_m)=\sum \frac{d^2 k}{(2 \pi)^2}[i \nu_m + \mu - \varepsilon(k)]^{-1} \equiv \sum \rho(\epsilon) [i \nu_m + \mu - \epsilon]^{-1}$$.

A defining feature of the Bose-Einstein distribution is that I need $$\varepsilon(k) - \mu > 0, \forall k$$. If this holon has a dispersion with half-bandwidth $$D= 1$$, this means that I need $$\mu < -1.0$$.

I can compute the BE-distribution exactly for a square lattice dispersion $$\varepsilon(k) = -2t [\cos(k_x) + \cos(k_y)]$$) and find a sensible $$n_b $$ versus $$\mu$$ curve. ($$t = 0.25$$ for $$D = 1$$ (Code supplied).

I can also compute $$n_b $$ versus $$\mu$$ curve for a "local" holon with fixed $$\epsilon = -1$$. However, if I try to construct a holon with such a $$\varepsilon(k)$$, I can never get a positive $$n_b$$.

I am unsure if this is a conceptual problem with working with bosonic Green's functions, or an error in my implementation in TRIQS, but I would be incredibly grateful if anyone could provide guidance.

The minimal code I can pair down to explain the issue is here

from triqs.gf import *
from triqs.gf.tools import *
from triqs.gf.block_gf import *
from triqs.plot.mpl_interface import *
import matplotlib.pyplot as plt
import numpy as np
from triqs_tprf.tight_binding import TBLattice
from triqs.lattice.tight_binding import dos 
from triqs_tprf.lattice import *
from triqs_tprf.lattice_utils import *
class Hamiltonian:
    def __init__(self, H, params):
        self.H = H
        self.kmesh = H.get_kmesh(n_k=(params["n_k"],params["n_k"], 1))
        self.e_k = H.fourier(self.kmesh)
        self.bz = H.bz
        dos_array = dos(self.H.tb, params["n_k"]*20, params["n_e"], "DOS")
        self.dos = dos_array[0]
def H_square(params):
    gf_struct = []
    for o in range(0, -params["n_orbs"]):
        gf_struct.append("up-" + str(o))
        gf_struct.append("dn-" + str(o))
    H = TBLattice(
    units = [(1, 0, 0), (0, 1, 0)],
    hopping = {
        ( 0,+1): -params["hopping"] * np.eye(2*params["n_orbs"]),
        ( 0,-1): -params["hopping"] * np.eye(2*params["n_orbs"]),
        (+1, 0): -params["hopping"] * np.eye(2*params["n_orbs"]),
        (-1, 0): -params["hopping"] * np.eye(2*params["n_orbs"]),
        },
    orbital_positions = [(0,0,0)]*2*params["n_orbs"],
    orbital_names = gf_struct,
    )
    Ham = Hamiltonian(H,params)
    return Ham
params = {
    "n_k": 64,
    "n_e": 1000,
    "beta": 10.0,
    "hopping": 0.25,
    "n_orbs": 1
}
mu_arr = np.linspace(-1.5,-1.0, 100)
nb = np.zeros(100)
dens = np.zeros(100)
dens_2 = np.zeros(100)
h_0 = H_square(params)
ek = np.array([np.linalg.eigvalsh(h_0.e_k(h_0.kmesh[i].value)) for i in range(params["n_k"]**2)])[:,0]
plt.figure(figsize=(8,8))
plt.plot(h_0.dos.eps, h_0.dos.rho)
plt.xlabel(r"$\epsilon$")
plt.ylabel(r"$\rho(\epsilon)$")
plt.show()
for (i,x) in enumerate(mu_arr):
    nb[i] = np.sum(1/(np.e**(params["beta"]*(ek-x))-1))/(2*np.pi)**2
iw_mesh = MeshDLRImFreq(beta=params["beta"], statistic='Boson', w_max= 60, eps = 1e-14, symmetrize = True)
G_iw = Gf(mesh = iw_mesh, target_shape=(1,1))
params["mu"] = 0
G_dummy = G_iw.copy()
G_k = G_dummy.copy()
mu_mat = np.zeros((1,1), dtype = complex)                                     
for i in range(1):                                                  
    mu_mat[i,i] = params["mu"]                                             
for i,rho in enumerate(h_0.dos.rho):
    G_k << rho*1.0*inverse(iOmega_n + mu_mat-h_0.dos.eps[i]) 
    G_dummy += G_k*2/params["n_e"]
G_iw.data[:,0,0] = G_dummy[0,0].data[:]
G_2 = G_iw.copy()
for (i, mu) in enumerate(mu_arr):
    G_iw.data[:,0,0] = G_dummy[0,0].data[:]
    G_2 << inverse(iOmega_n + 1.0001 + mu )
    G_iw << inverse(inverse(G_iw)+mu)
    G_dlr = make_gf_dlr(G_iw)
    G_2_dlr = make_gf_dlr(G_2)
    dens[i] = -G_dlr(0.0)[0,0].real-1
    dens_2[i] = -G_2_dlr(0.0)[0,0].real - 1

plt.figure(figsize=(8,8))
plt.plot(mu_arr,nb, label = "Analytical")
plt.plot(mu_arr,dens, label = "Local Green's Function")
plt.plot(mu_arr,dens_2, label = "Dispersive Green's Function")
plt.ylim(-0.2,2)
plt.xlabel(r"$\mu$")
plt.ylabel(r"$n_B(\mu)$")
plt.legend()
plt.show()