Physical background

Matsubara functions

A Matsubara function \(G(\mathrm i\nu_k)\) is a matrix-valued function of size \(N_{\text{orb}} \times N_{\text{orb}}\), where \(N_{\text{orb}}\) is the number of orbitals, and \(\{\mathrm i\nu_k, k\in\mathbb Z\}\) are the Matsubara frequencies:

\[\begin{split}\mathrm i\nu_k = \left\{ \begin{aligned} \mathrm i \frac{(2n+1)\pi}{\beta}, & \text{ for fermions},\\ \mathrm i \frac{2n\pi}{\beta}, & \text{ for bosons}. \end{aligned} \right.\end{split}\]

Here \(\beta\) is the inverse temperature.

The Matsubara function admits the following Lehmann representation:

\[\begin{equation} G(\mathrm i\nu_k) = \int_{-\infty}^{\infty} \frac{1}{\mathrm i\nu_k - w} \rho(w)\mathrm dw, \tag{1} \label{Lehmann} \end{equation}\]

Here \(\rho(w)\) is a matrix-valued function, and could either be sum of Dirac delta functions, or a continuous function.

The Matsubara function in the time domain, denoted as \(G(\tau)\), is related to \(G(\mathrm i\nu_k)\) by the Fourier transform:

\[G(\mathrm i\nu_k) = \int_{0}^{\beta} G(\tau)\mathrm{e}^{\mathrm i\nu_k\tau} \mathrm d\tau.\]

In other words, \(G(\mathrm i\nu_k)\) is the Fourier coefficient of \(G(\tau)\), and \(G(\tau)\) satisfies the Kubo-Martin-Schwinger (KMS) (anti-)periodic condition: \(G(\tau+\beta) = \pm G(\tau)\).

Pole fitting

It has been well known that any Matsubara functions \(G(\mathrm i\nu_k)\) (i.e. functions defined by \(\eqref{Lehmann}\)) could be well approximated by a finite sum of poles,

\[\begin{equation} G(\mathrm i\nu_k) \approx \sum_{p=1}^{N_p} \frac{M_p}{\mathrm i\nu_k - E_p}, \tag{2} \label{pole} \end{equation}\]

where \(N_p\) is the number of poles, \(M_p\) is a matrix-valued function, and \(E_p\) are the poles.

Our goal is to obtain the poles \(E_p\) and the weights \(M_p\) (also known as the residues) numerically. We hope to obtain such an approximation with \(N_p\) as small as possible.

Hybridization fitting

The hybridization fitting is a crucial subroutine in dynamical mean-field theory calculations. In DMFT, the hybridization function \(\Delta(\mathrm i\nu_k)\), which is obtained in each iteration of the DMFT self-consistent loop, is a matrix-valued Matsubara function of size \(N_{\text{orb}}\times N_{\text{orb}}\), where \(N_{\text{orb}}\) is the number of impurity orbitals. After obtaining \(\Delta(\mathrm i\nu_k)\), one conduct the following hybridization fitting:

\[\begin{equation} \Delta(\mathrm i\nu_k) = \sum_{j=1}^{N_b} \frac{V_jV_j^{\dagger}}{\mathrm i\nu_k - E_j}. \tag{3} \label{hybfit} \end{equation}\]

Here \(E_j\in\mathbb R\) are bath energies, and \(V_j\) are impurity-bath coupling coefficients, and are vectors of size \((N_{\text{orb}},1)\). Using \(E_p\) and the vectors \(V_p\), one can construct the new impurity-bath Hamiltonian \(\hat H\):

\[\hat H =\hat H_{\text{loc}} + \sum_{j=1}^{N_b} E_j \hat c_j^{\dagger} \hat c_j + \sum_{j=1}^{N_b} \sum_{\alpha=1}^{N_{\mathrm{orb}}} (V_{j\alpha}\hat c_j^{\dagger}\hat c_{\alpha} + \text{h.c.}).\]

Here \(\hat H_{\text{loc}}\) is the local impurity Hamiltonian.

Note that \(N_b\) is the number of bath orbitals, and is also the number of terms that are used in the hybridization fitting \(\eqref{hybfit}\).

Analytical continuation

The Matsubara function \(G(\mathrm i\nu_k)\) (as defined in Eq. \(\eqref{Lehmann}\)) could be extended to a function that is defined on the entire complex plane:

\[G(z) = \int_{-\infty}^{\infty} \frac{1}{z - w} \rho(w)\mathrm dw.\]

\(G(z)\) is analytic in the upper half-plane, and in the lower half plane, but has a branch cut in the real axis. In the special case that \(\rho(w)\) is sum of Dirac delta functions, \(G(z)\) is a rational function and is analytic in the entire complex plane except for the poles.

More importantly, \(\rho(w)\) and \(G(z)\) are related by the following formula:

\[\rho(w) = -\frac{1}{\pi} \lim_{\eta\to 0^+} \operatorname{Im} G(w+\mathrm i\eta).\]

Matsubara Green’s function and spectral function

If \(G(\mathrm i\nu_k)\) is a single-particle Green’s function, (or similarly, any correlation functions), the Lehmann’s representation states that it can be written as a sum of poles:

\[G_{ij}(\mathrm i\nu_k) =\frac{1}{Z} \sum_{r, s} \frac{\left\langle\Psi_s\left|\hat{c}_i\right| \Psi_r\right\rangle\left\langle\Psi_r\left|\hat{c}_j^{\dagger}\right| \Psi_s\right\rangle}{z+E_s-E_r}\left(\mathrm{e}^{-\beta E_s} \mp \mathrm{e}^{-\beta E_r}\right),\]

where \(\hat{c}_i, \hat{c}_i^{\dagger}\) are the annihilation and creation operator for the \(i\)-th orbital, \(\left|\Psi_s\right\rangle\) is the \(s\)-th eigenvector of the Hamiltonian \(\hat H\) with energy \(E_s\).

By grouping the \((r,s)\) indices into a single index \(p\), we have a more compact form:

\[G(\mathrm i\nu_k) = \sum_{p=1}^{N_p} \frac{V_pV_p^{\dagger}}{\mathrm i\nu_k - E_p},\]

and in the thermodynamic limit, the sum becomes an integral:

\[G(\mathrm i\nu_k) = \int_{-\infty}^{\infty} \frac{1}{\mathrm i\nu_k - w} \rho(w),\]

Here \(\rho(w)\) is positive semi-definite for all \(w\). In the scalar case (\(N_{\text{orb}}=1\)), \(\rho(w)\) is the spectral function. In the matrix case, the spectral function is:

\[{\mathrm{Spec}}(w) = \sum_{i=1}^{N_{\text{orb}}}\rho_{ii}(w).\]

Using the analytic continuation properties, the spectral function could be evaluated by calculating the Green’s function on the real axis with infinitesimal broadening:

\[\mathrm{Spec}(w) = -\frac{1}{\pi} \lim_{\eta\to 0^+} \sum_{i=1}^{N_{\text{orb}}} \operatorname{Im} (G_{ii}(w+\mathrm i\eta)).\]