docs: def gf_lehmann(eig_e, eig_states, d_dag, beta, omega, d=None):
      """Outputs the lehmann representation of the greens function
         omega has to be given, as matsubara or real frequencies"""
  
      ew = np.exp(-beta*eig_e)
      zet = ew.sum()
      G = np.zeros_like(omega)
      basis_create = np.dot(eig_states.T, d_dag.dot(eig_states))
      if d is None:
          tmat = np.square(basis_create)
      else:
          tmat = np.dot(eig_states.T, d.T.dot(eig_states))*basis_create
  
      tmat *= np.add.outer(ew, ew)
      gap = np.add.outer(-eig_e, eig_e)
  
      N = eig_e.size
      for i, j in product(range(N), range(N)):
          G += tmat[i, j] / (omega + gap[i, j])
      return G / zet