User Guide

This section provides a brief introduction to the theory behind gem and worked examples showing how to set up and run calculations.

Theoretical Background

A proper presentation of the method has been given in many works for zero [1] and finite [2] temperature. Here, we give a brief recap.

The Ghost Gutzwiller Approximation

The Gutzwiller approximation (GA) is a variational method for lattice models of strongly correlated electrons with local interaction [3]. The trial state is obtained by applying a linear map \(\hat{P} = \prod_i \hat{P}_i\) (often called a projector, despite not being one strictly speaking) to an uncorrelated Slater determinant \(|\Psi_0\rangle\). Each on-site operator \(\hat{P}_i\) re-weights the local many-body configurations, suppressing or enhancing occupancies relative to the non-interacting reference and provided the first explanation for the Mott transition via a renormalization of the hopping parameters.

The standard GA is limited because its embedding Hamiltonian is an interacting single-impurity model with at most one bath orbital per correlated orbital. This restricts the representable spectral structures, missing high-energy features like Hubbard bands, and prevents a systematic improvement of the approximation toward a dynamical description of strongly correlated systems as in Dynamical Mean-Field Theory (DMFT) [4].

Ghost-GA [1] extends the embedding Hamiltonian by adding auxiliary (ghost) orbitals to the uncorrelated Slater determinant \(|\Psi_0\rangle\), considering a number of electronic levels that is a multiple \(B\) of the physical ones.

• At \(B = 1\) the method reduces to the standard GA.

• As \(B \to \infty\) it has been proved that the method converges to DMFT [2].

• For finite \(B > 1\) it provides a controlled, systematically improvable approximation.

The variational parameters are the bath hybridisation matrix and the on-site embedding energies. One of the key advantages of this method is that its self-consistency is rooted in static expectation values rather than dynamical correlators (as in DMFT), enabling much faster calculations of correlated matter.

In Ref. [2] it has been shown that Ghost-GA admits a free-energy functional formulation that connects it with DMFT and allows for a natural finite-temperature extension. This framework embeds the correlated site in an effective Anderson impurity model at temperature \(T\), and the self-consistency equations remain rooted in static — now thermal — expectation values. The method interpolates smoothly between zero-temperature ghost-GA (\(T \to 0\)) and conventional DMFT (\(B \to \infty\)), retaining the low computational cost of ghost-GA for moderate \(B\).

The ghost embedding functional

In Ref. [2], a free-energy functional formulation of the method has been derived, which we briefly recap here. Consider a general multi-orbital Hubbard Hamiltonian written as

(1)\[\hat H = \hat H_0 + \hat H_{\mathrm{int}}\]

where the interaction part is assumed to be local,

(2)\[\hat H_{\mathrm{int}} = \sum_{i=1}^{\mathcal{N}} \hat H^{i}_{\mathrm{int}}\]

and the one-body part is

(3)\[\hat H_0 = \sum_{i,j=1}^{\mathcal{N}} \sum_{\alpha=1}^{\nu_i} \sum_{\beta=1}^{\nu_j} [h_0]_{i\alpha,j\beta}\, c^\dagger_{i\alpha} c_{j\beta} \, .\]

Here \(i=1,\ldots, \mathcal{N}\) labels system fragments, each hosting \(\nu_i\) physical fermionic modes with annihilation operators \(c_{i\alpha}\) (\(\alpha=1,\ldots,\nu_i\)), including both spin and orbital degrees of freedom.

We represent the one-body matrix \(h_0\) in block form as

(4)\[\begin{split}h_0 = \begin{pmatrix} \epsilon_1 & t_{12} & \dots & t_{1\mathcal{N} } \\ t_{21} & \epsilon_2 & \dots & \vdots \\ \vdots & \vdots & \ddots & \vdots \\ t_{\mathcal{N}1} & \dots & \dots & \epsilon_{\mathcal{N}} \end{pmatrix}\end{split}\]

where \(\epsilon_i \in \mathbb{C}^{\nu_i \times \nu_i}\) are Hermitian local one-body blocks and \(t_{ij} \in \mathbb{C}^{\nu_i \times \nu_j}\) (\(i \neq j\)) are the hopping blocks satisfying \(t_{ji}=t_{ij}^{\dagger}\) as a consequence of the Hermiticity of \(h_0\).

We also define the hopping matrix \(t\) as the block matrix collecting the off-diagonal blocks of \(h_0\):

(5)\[\begin{split}t = \begin{pmatrix} \mathbf{0} & t_{12} & \dots & t_{1\mathcal{N} } \\ t_{21} & \mathbf{0} & \dots & \vdots \\ \vdots & \vdots & \ddots & \vdots \\ t_{\mathcal{N}1} & \dots & \dots & \mathbf{0} \end{pmatrix} \, ,\end{split}\]

and the corresponding block-diagonal local one-body matrix

(6)\[\begin{split}\epsilon = \begin{pmatrix} \epsilon_1 & \mathbf{0} & \dots & \mathbf{0} \\ \mathbf{0} & \epsilon_2 & \dots & \vdots \\ \vdots & \vdots & \ddots & \vdots \\ \mathbf{0} & \dots & \dots & \epsilon_{\mathcal{N}} \end{pmatrix}, \qquad h_0 = \epsilon + t \, .\end{split}\]

Grouping the local one-body terms and the local interactions into a single local operator

(7)\[\hat H^{i}_{\mathrm{loc}} \big[c^\dagger_{i\alpha}, c_{i\alpha}\big] = \sum_{\alpha,\beta=1}^{\nu_i} [\epsilon_i]_{\alpha\beta}\, c^\dagger_{i\alpha} c_{i\beta} + \hat H^{i}_{\mathrm{int}} \, ,\]

the full Hamiltonian at Eq. (1) can be rewritten as:

(8)\[\begin{split}\hat H = \sum_{\substack{i,j=1 \\ i \neq j}}^{\mathcal{N}} \sum_{\alpha=1}^{\nu_i} \sum_{\beta=1}^{\nu_j} [t_{ij}]_{\alpha\beta}\, c^\dagger_{i\alpha} c_{j\beta} + \sum_{i=1}^{\mathcal{N}} \hat H^{i}_{\mathrm{loc}} \big[c^\dagger_{i\alpha}, c_{i\alpha}\big] \, ,\end{split}\]

The method is built on the following dynamical functional:

(9)\[\begin{split}\begin{align} \mathcal{L}_\beta [ {\Sigma_i(i\omega_n)} , {\Delta_i(i\omega_n) }] = -& T\sum_{n}e^{i\omega_n0^+} \text{Tr} \ln \big(i \omega_n \mathbf{1}-h_0-\Sigma(i \omega_n ) \big) \nonumber \\ +& T\sum_{i=1}^{\mathcal{N}} \sum_{n}e^{i\omega_n0^+}\text{Tr} \ln \big(i \omega_n \mathbf{1} - \epsilon_i - \Delta_i(i \omega_n ) - \Sigma_i(i \omega_n )\big) \nonumber \\ +& \sum_{i=1}^{\mathcal{N}} \Omega_{\rm imp}^i[\Delta_i(i \omega_n )] \end{align}\end{split}\]

Where \(\Sigma(i\omega_n)= \text{diag}(\Sigma_1(i\omega_n), \ldots, \Sigma_{\mathcal{N}}(i\omega_n))\) The saddle point equations of this functional are the DMFT self-consistency conditions. Introducing the following parametrization of the self-energy and hybridization functions:

(10)\[\begin{split}\begin{align} \Sigma_i(z) &= z\mathbf{1}_{\nu_i} -\epsilon_i -\Big[ \mathcal{R}_i^\dagger \big( z\mathbf{1}_{B\nu_i}-\Lambda_i\big)^{-1}\mathcal{R}_i \Big]^{-1} \\ \Delta_i(z) &= \mathcal{D}_i^{T}\big(z\mathbf{1}+\Lambda_i^c\big)^{-1}\mathcal{D}_i^{*} \end{align}\end{split}\]

where \(\mathcal{R}_i, \mathcal{D}_i\) are complex matrices of shape \(B \nu_i \times \nu_i\) and \(\Lambda_i, \Lambda_i^c\) are complex hermitian matrices, the functional can be rewritten as:

(11)\[\begin{split}\begin{align} \mathcal{L}_\beta[\mathcal{R},\Lambda,\mathcal{D},\Lambda^c] =& \Omega_{\rm qp}[\mathcal{R},\Lambda] \nonumber \\ -&\sum_{i=1}^{\mathcal{N}}\Omega_{0,{\rm emb}}^i[\mathcal{D}_i,\Lambda_i^c;\mathcal{R}_i,\Lambda_i] \nonumber \\ +&\sum_{i=1}^{\mathcal{N}}\Omega_{\rm emb}^i[\mathcal{D}_i,\Lambda_i^c] \end{align}\end{split}\]

where \(\Omega_{\alpha}\) are the free-energies associated to the following Hamiltonian:

(12)\[\begin{split}\begin{align} \hat{H}_{\text{qp}} &= \sum_{ij} \sum_{a=1}^{\mathcal{B} \nu_i} \sum_{b=1}^{\mathcal{B} \nu_j} f^\dagger_{i a} \, [ \delta_{ij} (\Lambda^i)_{ ab} + \sum_{\alpha \beta} R_{i, a \alpha} \, t^{ij}_{\alpha\beta} \, R^\dagger_{j,\beta b} ] \, f_{ j b}\\ \hat{H}^{i}_{\text{emb}} &= \hat{H}^{i}_{\text{int}}[ c^\dagger_i ,c_i ] + \sum_{ a=1}^{\mathcal{B} \nu_i} \sum_{ \alpha=1 }^{\nu_i} [ (D^i)_{ a \alpha } \, b^\dagger_{i, a} c_{i, \alpha} + \text{h.c.} ] + \sum_{\alpha\beta=1}^{\mathcal{B} \nu_i} (\Lambda_c^i)_{\alpha\beta} \, f^\dagger_\alpha b_\beta \\ \hat{H}^{i}_{0,\text{emb}} &= \sum_{a,b=1}^{B\nu_i} \big[\Lambda^i\big]_{ab}\, f^\dagger_{ia} f_{ib} + \sum_{a,b=1}^{B \nu_i} \Big(\big[D^i R_i^{T}\big]_{ab}\, f^\dagger_{ib} b_{ia}+\text{H.c.}\Big)+ \sum_{a,b=1}^{B\nu_i}\big[\Lambda^i_c\big]_{ab}\,b^\dagger_{ib} b_{ia} \end{align}\end{split}\]

The saddle point equations of the functional at \(T=0\) retrieve the ghost-GA self-consistency equations and generalize them to finite temperature, while the limit \(B \to \infty\) retrieves the DMFT functional and self-consistency equations.

Self-Consistency Loop

The self-consistency loop of the method is based on the saddle point equations of the functional at Eq. (11), which are given in the following:

(13)\[\frac{\partial \mathcal{L}_\beta}{\partial [\mathcal{R}_i]_{a \alpha} } = 0 \Rightarrow \sum_{j=1}^{\mathcal{N}} \sum_{\beta=1}^{\nu_j} \sum_{b=1}^{B \nu_j} [t^{ij}]_{\alpha \beta} [\mathcal{R}^\dagger_j ]_{\beta b} \langle f^\dagger_{ia} f_{j b} \rangle_{\text{qp}} = \sum_{b=1}^{B \nu_i} [\mathcal{D}_i]_{b \alpha} \langle f^\dagger_{i a} b_{i b} \rangle_{\text{0,emb,i}}\]

(14)\[\frac{\partial \mathcal{L}_\beta}{\partial [\Lambda_i]_{a b}} = 0 \Rightarrow \langle f^\dagger_\alpha f_\beta \rangle_{\text{qp}} = \langle f^\dagger_\alpha f_\beta \rangle_{\text{0,emb,i}}\]

(15)\[\frac{\partial \mathcal{L}_\beta}{\partial [\mathcal{D}_i]_{a \alpha} } = 0 \Rightarrow \langle c^\dagger_{i \alpha} b_{ib} \rangle_{\text{emb,i}} = \langle \Big( \sum_{a=1}^{B\nu_i}[ \mathcal{R}_i ]_{a \alpha} f^\dagger_{ia} \Big)b_{ib} \rangle_{\text{0,emb,i}}\]

(16)\[\frac{\partial \mathcal{L}_\beta}{\partial [\Lambda_i^c]_{a b} } = 0 \Rightarrow \langle b_\alpha b^\dagger_\beta \rangle_{\text{emb,i}} = \langle b_\alpha b^\dagger_\beta \rangle_{\text{0,emb,i}}\]

A typical ghost-GA calculation for a single correlated fragment proceeds as follows:

1. Build the dispersion — construct the array eks of shape (n_k, nimp, nimp) and the corresponding k-point weights wks.

2. Create the lattice — Lattice(eks, wk_list=wks).

3. Choose a solver and create a fragment — Fragment(nimp, nbath, eloc, Utensor, solver).

4. Iterate until convergence:

   1. lattice.solve_qp([fragment], T=T) — solve the quasiparticle problem and update the quasiparticle target expectation values.

   2. fragment.update_hybridization(T=T) — Update the parameters \(D\) and \(\Lambda_c\) of the hybridization function via the saddle point equations (15) and (16).

   3. fragment.solve_impurity(mu, T=T) — solve the embedding Hamiltonian and update the impurity target expectation values.

   4. fragment.update_self_energy(T=T) — Update the parameters \(R\) and \(\Lambda\) of the self-rnergy via the saddle point equations (13) and (14). parameters \(R\) and \(\Lambda\) from the density matrix.

   5. Mix old and new parameters and check convergence.

5. Extract observables — fragment.compute_Z() for the quasiparticle weight, fragment.denMat for the impurity density matrix.

Examples