Algorithms

Definition of the susceptibility

The momentum-dependent susceptibility computed by ChiQ is defined by

(1)\[\chi_{ij,kl}(\boldsymbol{q}, i\Omega_m) = \frac{1}{N} \sum_{\boldsymbol{r}} \int_{0}^{\beta} d\tau \langle O_{ij}(\boldsymbol{r}, \tau) O_{lk}(\boldsymbol{0}, 0) \rangle e^{i\Omega_m \tau} e^{-i\boldsymbol{q} \cdot \boldsymbol{r}},\]

where \(\boldsymbol{q}\) represents the momentum, \(\Omega_m\) represents the bosonic Matsubara frequency, \(\beta\) is the inverse temperature, and \(\boldsymbol{r}\) denotes the lattice positions (the origin of each unit cell) with \(N\) being their total number. The superscript \(i, j, k, l\) indicates the combined spin and orbital indices. Here, the term “orbital” includes the site within a unit cell, as defined in the tight-binding model in the Wannier90 format. The angle brackets denote the thermal average. The operator \(O_{ij}(\boldsymbol{r}, \tau)\) is defined by

\[O_{ij}(\boldsymbol{r}) = c_{i}^{\dagger}(\boldsymbol{r}) c_{j}(\boldsymbol{r}),\]

and the argument \(\tau\) indicates the time evolution in the Heisenberg picture.

Calculation of the susceptibility

Bethe-Salpeter equation (BSE)

This part is based on Section II.A in Ref. [Otsuki et al. 2019].

In DMFT, we compute the two-particle Green’s function \(X_{ij,kl}(i\omega_n, i\omega_{n'}; i\Omega_m)\), which has two additional fermionic Matsubara frequencies \(i\omega_n\) and \(i\omega_{n'}\). The susceptibility \(\chi_{ij,kl}(\boldsymbol{q}, i\Omega_m)\) in Eq. (1) is obtained by summing over the fermionic frequencies as

(2)\[\chi_{ij,kl}(\boldsymbol{q}, i\Omega_m) = T \sum_{nn'} X_{ij,kl}(i\omega_n, i\omega_{n'}; \boldsymbol{q}, i\Omega_m).\]

We introduce a matrix notation by combining spin-orbital indices and the frequency index as

\[[\hat{X}(\boldsymbol{q}, i\Omega_{m})]_{(ijn),(kln')} \equiv X_{ij,kl}(i\omega_{n}, i\omega_{n'}; \boldsymbol{q}, i\Omega_{m}).\]

Then, the Bethe-Salpeter (BS) equation that \(\hat{X}(\boldsymbol{q}, i\Omega_{m})\) follows is given by

\[\hat{X}^{-1}(\boldsymbol{q}, i\Omega_{m}) = \hat{X}_0^{-1}(\boldsymbol{q}, i\Omega_{m}) - \hat{\Gamma}_{\text{loc}}(i\Omega_{m}),\]

where \(\hat{X}_0(\boldsymbol{q}, i\Omega_{m})\) is the bare two-particle Green’s function

(3)\[[\hat{X}_0(\boldsymbol{q}, i\Omega_{m})]_{(ijn),(kln')} = -\delta_{nn'} \frac{1}{N} \sum_k G_{ki}(\boldsymbol{k}, i\omega_{n}) G_{jl}(\boldsymbol{k} + \boldsymbol{q}, i\omega_n + i\Omega_{m}).\]

Here, \(N\) denotes the number of lattice sites, and \(G_{ki}(\boldsymbol{k}, i\omega_n)\) is the single-particle Green function within DMFT. \(\Gamma_{\text{loc}}(i\Omega_m)\) is the vertex part. The vertex part is local (no dependence on \(\boldsymbol{q}\)) in DMFT. The vertex part can be obtained from the local two-particle Green’s function \(\hat{X}_{\text{loc}}(i\Omega_m)\) defined in the effective impurity problem in the DMFT:

\[\hat{X}_\text{loc}^{-1}(i\Omega_{m}) = \hat{X}_\text{0,loc}^{-1}(i\Omega_{m}) - \hat{\Gamma}_{\text{loc}}(i\Omega_{m}).\]

Eliminating \(\Gamma_{\text{loc}}(i\Omega_{m})\) from the above two BS equations, we obtain the final equation

(4)\[\hat{X}(\boldsymbol{q}, i\Omega_{m}) = \left[ \hat{X}_\text{loc}(i\Omega_{m})^{-1} - \hat{X}_\text{0,loc}(i\Omega_{m})^{-1} + \hat{X}_{0}(\boldsymbol{q}, i\Omega_{m})^{-1} \right]^{-1}.\]

Implementation notes:

• \(\hat{X}_0(\boldsymbol{q}, i\Omega_{m})\) [Eq. (3)] is computed in dcore_chiq.py.

• \(\hat{X}_\text{loc}(i\Omega_{m})\) is computed in the impurity solver, which is called from dcore_chiq.py.

• \(\hat{X}(\boldsymbol{q}, i\Omega_{m})\) [Eq. (4)] and \(\hat{\chi}(\boldsymbol{q}, i\Omega_{m})\) [Eq. (2)] are computed in C++ part, which is called from chiq_main.py.

• The fermionic Matsubara frequencies are truncated. The number of fermionic Matsubara frequency is specified by the parameter num_wf. The maximum frequency is \((2 \mathtt{num\_wf} + 1) \pi T\).

• The size of the matrices is \((2N_{\text{orb}})^2 \times \mathtt{num\_wf}\), where \(N_{\text{orb}}\) is the number of the total orbitals of the correlated shells in the unit cell.

Strong-coupling limit formula (SCL)

The desciption below is based on Section III.C in Ref. [Otsuki et al. 2024].

The strong-coupling limit formula (SCL) is a simplified approach to compute the momentum-dependent static susceptibilities in DMFT. The derivation is based on the strong-coupling feature of the two-particle Green’s function. The validity of the SCL formula has been demonstrated for the single-orbital and two-orbital Hubbard models [Otsuki et al. 2019] and for the multipolar ordering in CeB₆ [Otsuki et al. 2024].

In the original paper [Otsuki et al. 2019], three levels of approximations termed SCL1, SCL2, and SCL3 were proposed. ChiQ package provides SCL3, which is the simplest method among the three and does not require the vertex part.

The momentum-dependent static (\(\Omega_m=0\)) susceptibility in the SCL formula is given by

(5)\[\hat{\chi}_\text{SCL}(\boldsymbol{q}) = \left[ \hat{\chi}_\text{loc}^{-1} - \hat{I}_\text{SCL}(\boldsymbol{q}) \right]^{-1},\]

where hat indicates the matrix form defined as

\[[\hat{\chi}(\boldsymbol{q})]_{(ij),(kl)} \equiv \chi_{ij,kl}(\boldsymbol{q}).\]

\(\hat{\chi}_\text{loc}\) is the local susceptibility computed in the effective impurity model. \(\hat{I}_\text{SCL}(\boldsymbol{q})\) is the effective intersite interaction defined by

(6)\[\hat{I}_\text{SCL}(\boldsymbol{q}) = T \sum_{n} \phi(i\omega_n)^2 \hat{\Lambda}(i\omega_n, \boldsymbol{q}).\]

The function \(\hat{\Lambda}(i\omega_n, \boldsymbol{q})\) is computed from the bare two-particle Green’s function in the effective impurity and lattice models:

\[\hat{\Lambda}(i\omega_n, \boldsymbol{q}) = \hat{X}_\text{0,loc}(i\omega_n)^{-1} - \hat{X}_{0}(i\omega_n, \boldsymbol{q})^{-1}.\]

\(\phi(i\omega_n)\) is given by

(7)\[\phi(i\omega_n) = \frac{1}{i\omega_n + \Delta_{-}} - \frac{1}{i\omega_n - \Delta_{+}},\]

where \(\Delta_{-}\) and \(\Delta_{+}\) are the excitation energy from \(n\)-electron state to \((n-1)\)-electron states and \((n+1)\)-electron states, respectively. In the single-band Hubbard model, \(\Delta_{-}\) and \(\Delta_{+}\) correspond to the position of the lower and upper Hubbard bands, respectively.

Implementation notes:

• \(\hat{X}_0(\boldsymbol{q}, i\Omega_{m})\) [Eq. (3)] is computed in dcore_chiq.py.

• \(\hat{\chi}_\text{SCL}(\boldsymbol{q})\) and \(\hat{I}_\text{SCL}(\boldsymbol{q})\) [Eqs. (5) and (6)] are computed in the script calc_Iq_scl.py (chiq_main.py is not used).

• The value of \(\Delta_{-}\) and \(\Delta_{+}\) are specified by the parameter delta_minus and delta_plus, respectively.

• The size of the matrices is \((2N_{\text{orb}})^2\).

Random phase approximation (RPA)

The random-phase approximation (RPA) is a simple weak-coupling approximation to compute the momentum-dependent susceptibility. The RPA susceptibility is given by

(8)\[\hat{\chi}_\text{RPA}(\boldsymbol{q}, i\Omega_m) = \left[ \hat{\chi}_0(\boldsymbol{q}, i\Omega_{m})^{-1} - \hat{\gamma}_0 \right]^{-1},\]

where \(\hat{\chi}_0(\boldsymbol{q}, i\Omega_{m})\) is the bare susceptibility, which is calculated from \(\hat{X}_0\) in Eq. (3) by

(9)\[[\hat{\chi}_0(\boldsymbol{q}, i\Omega_{m})]_{(ij),(kl)} = T \sum_{n} [\hat{X}_0(\boldsymbol{q}, i\Omega_{m})]_{(ijn),(kln)}.\]

\(\hat{\gamma}\) is the symmetrized vertex part given by

\[[ \hat{\gamma}_0 ]_{(ij),(kl)} = U_{ijkl} - U_{ijlk}.\]

Implementation notes:

• \(\hat{\gamma}_0\) is generated in dcore_chiq.py.

• \(\hat{X}_0(\boldsymbol{q}, i\Omega_{m})\) [Eq. (3)] is computed in dcore_chiq.py.

• \(\hat{\chi}_\text{RPA}(\boldsymbol{q}, i\Omega_{m})\) and \(\hat{\chi}_0(\boldsymbol{q}, i\Omega_{m})\) [Eqs. (8) and (9)] are computed in C++ part, which is called from chiq_main.py.

• The size of the matrices is \((2N_{\text{orb}})^2\).

Renormalized random-phase approximation (RRPA)

This part is based on Section II.B in Ref. [Otsuki et al. 2019].

The RPA above overestimates the fluctuations. Improvement to the RPA is made by including correlations from DMFT. By taking the Matsubara summation independently for each \(X\) in Eq. (4), we obtain the renormalized RPA (RRPA) susceptibility as

(10)\[\hat{\chi}_\text{RRPA}(\boldsymbol{q}, i\Omega_{m}) = \left[ \hat{\chi}_{0}(\boldsymbol{q}, i\Omega_{m})^{-1} - \hat{\chi}_\text{0,loc}(i\Omega_{m})^{-1} + \hat{\chi}_\text{loc}(i\Omega_{m})^{-1} \right]^{-1}.\]

Compared with the RPA formula in Eq. (8), the bare interaction \(\hat{\gamma}\) is replaced by \(\hat{\chi}_\text{0,loc}(i\Omega_{m})^{-1} - \hat{\chi}_\text{loc}(i\Omega_{m})^{-1}\), which takes the local correlations into account.

Implementation notes:

• \(\hat{X}_0(\boldsymbol{q}, i\Omega_{m})\) [Eq. (3)] is computed in dcore_chiq.py.

• \(\hat{\chi}_\text{RRPA}(\boldsymbol{q}, i\Omega_{m})\) and \(\hat{\chi}_0(\boldsymbol{q}, i\Omega_{m})\) [Eqs. (10) and (9)] are computed in C++ part, which is called from chiq_main.py.

• The size of the matrices is \((2N_{\text{orb}})^2\).

Intersite interactions

This part is based on Section III.B in Ref. [Otsuki et al. 2024].

We can estimate the intersite interactions inversly from \(\hat{\chi}(\boldsymbol{q})\) computed by solving the BSE. Assuming that the final results for \(\hat{\chi}(\boldsymbol{q})\) was obtained from a localized model within the mean-field approximation, we can define the momentum-dependent interactions \(\hat{I}(\boldsymbol{q})\) by

(11)\[\hat{I}(\boldsymbol{q}) = \hat{\chi}_\text{loc}^{-1} - \hat{\chi}(\boldsymbol{q})^{-1}.\]

We note that, in the SCL, we directly compute \(\hat{I}(\boldsymbol{q})\) by Eq. (6).

Performing the Fourier transform, we obtain the interaction coefficients in the real space

(12)\[\hat{I}(\boldsymbol{r}) = \frac{1}{N} \sum_{\boldsymbol{q}} \hat{I}(\boldsymbol{q}) e^{i\boldsymbol{q} \cdot \boldsymbol{r}},\]

where the summation of \(\boldsymbol{q}\) is taken over the Brillouin zone. This quantity defines the effective localized Hamiltonian

\[\mathcal{H}_\text{eff} = -\sum_{\boldsymbol{r} \boldsymbol{r}'} \sum_{ijkl} O_{ij}(\boldsymbol{r}) I_{ij,kl}(\boldsymbol{r} - \boldsymbol{r}') O_{lk}(\boldsymbol{r}').\]

Implementation notes:

• \(\hat{I}(\boldsymbol{q})\) [Eq. (11)] is computed in calc_Iq.py.

• Calculation of Eq. (11) is unstable when the charge fluctuation is tiny. The script calc_Iq.py takes a special care of the charge fluctuation to avoid the instability. See Appendix B in Ref. [Otsuki et al. 2024] for details.

• chiq_main.py also output \(\hat{I}(\boldsymbol{q})\), but using calc_Iq.py is recommended because the special treatment of the charge fluctuation is done only in calc_Iq.py.

• \(\hat{I}(\boldsymbol{r})\) [Eq. (12)] is computed in chiq_fft.py.

Analyzing the susceptibility

\(\chi_{ij,kl}(\boldsymbol{q}, i\Omega_m)\) is a complicated object with many components. This section describes how to analyze it to extract the physical information. For simplicity, we consider only the static component, \(\Omega_{m}=0\).

Eigenvalues

A simple way is to diagonalize \(\chi_{ij,kl}(\boldsymbol{q})\) by taking the combined index \((ij)\) as row and \((kl)\) as column.

(13)\[\chi_{\xi}(\boldsymbol{q}) = \sum_{ijkl} U^{(\xi)}_{ij}(\boldsymbol{q})^{\ast} \chi_{ij,kl}(\boldsymbol{q}) U^{(\xi)}_{kl}(\boldsymbol{q}).\]

The eigenvalues \(\chi_{\xi}(\boldsymbol{q})\) represent physical susceptibilities. The eigenvectors represent the physical quantitiy corresponding to the fluctuation.

Impelementation notes:

• chiq_post.py outputs eigenvalues by default (mode = 'eigen').

• The eigenvalues are sorted in descending order by default. This behavior can be changed by parameter order.

• The eigenvectors are output by setting vector = true (default is false).

Linear combinations

We can define the spin, charge, orbital, or multipole operator by taking linear comination of the density operator \(O_{\gamma} (\boldsymbol{r})\) as

\[O_{\gamma}(\boldsymbol{r}) = \sum_{ij} C^{(\gamma)}_{ij} O_{ij}(\boldsymbol{r}).\]

With this coefficient \(C^{(\gamma)}_{ij}\), we can transform \(\chi_{ij,kl}(\boldsymbol{q})\) as follows:

(14)\[\chi_{\gamma\gamma'}(\boldsymbol{q}) = \sum_{ijkl} C^{(\gamma)}_{ij} \chi_{ij,kl}(\boldsymbol{q}) C^{(\gamma') \ast}_{kl}.\]

Implementation notes:

• chiq_post.py outputs the diagonal componenents \(\chi_{\gamma\gamma}(\boldsymbol{q})\) when mode = 'linear_combination'.

• The coefficients \(C^{(\gamma)}_{ij}\) are input by an external file.

• The number of coefficients can be smaller than the total number of operators.

References

• J. Otsuki, K. Yoshimi, H. Shinaoka, and Y. Nomura, “Strong-coupling formula for momentum-dependent susceptibilities in dynamical mean-field theory”, Phys. Rev. B, 99, 165134 (2019).

• J. Otsuki, K. Yoshimi, H. Shinaoka, and H. Jeschke, “Multipolar ordering from dynamical mean field theory with application to CeB₆“, Phys. Rev. B, 110, 035104 (2024).