Susceptibility of the Hubbard-Zeeman model

In this tutorial we will compute the dynamic, momentum-resolved generalized susceptibility \(\chi_{abcd}(\omega,\mathbf{q})\) of the Hubbard-Zeeman model on the square lattice, along the lines of Commun. Phys. 6, 289 (2023). Unlike previous tutorials, in this tutorial we will compute and look at the full tensor structure of the susceptibility. Both the BSE and DBSE will be used.

The Hubbard-Zeeman Hamiltonian has the form

\[H = -t \sum_{\langle i,j \rangle, \sigma} c^\dagger_{i\sigma} c_{j\sigma} + U \sum_i \hat{n}_{i\uparrow} \hat{n}_{i\downarrow} + B \sum_i(\hat{n}_{i\uparrow}-\hat{n}_{i\downarrow})\]

with nearest neighbor hopping parameter \(t\), a Hubbard interaction \(U\), and a Zeeman magnetic field \(B\) which lifts the degeneracy between spin up (\(\sigma=\uparrow\)) and spin down (\(\sigma=\downarrow\)). The non-degenerate spin components induces, in turn, a richer spin structure of the susceptibility including non-trivial (spin) off-diagonal elements, e.g. \(\chi_{\downarrow\uparrow\uparrow\downarrow} \,\, (\omega, \mathbf{q}) = \chi_{S^- S^+}(\omega, \mathbf{q})\ne 0\).

Here we solve the model within DMFT using TRIQS/cthyb as an impurity solver for the DMFT self-consistency loop, and then compute the two-particle vertices using W2Dynamics. The overall structure of the scripts is similar to previous tutorials:

common.py is a module with general functionality,
calc_sc_dmft.py runs the DMFT self-consistency loop, and
calc_g2.py, calc_tri.py, calc_chi.py measure the three required correlation functions.
Finally, calc_bse.py evaluates the Bethe-Salpeter equation with different fermionic cut-offs nwf.

If you want to try the dual Bethe-Salpeter equation without running the impurity solver yourself, you can download a dataset from the Zenodo repository.

Convergence of Susceptibility

The following plot script shows the convergence with the fermionic frequency box size \(N_\nu\) (nwf) of the normal and dual Bethe-Salpeter equation, in color and black respectively. We focus on one particular component of the orbital tensor and one non-zero bosonic frequency \(\Omega=3\). The linear scaling with \(1/N_\nu\) for the normal Bethe-Salpeter equation is clearly visible, which makes extrapolation to \(N_\nu \rightarrow \infty\) necessary to obtain a quantitative result. On the other hand, for the dual Bethe-Salpeter equation, the calculations with finite frequency box size are very close to the extrapolated DBSE result.

In this plot, the extrapolated values of the BSE and DBSE are not identical, due to non-negligible Monte Carlo and self-consistency errors. However, the BSE and DBSE does converge to the same result once

the Monte Carlo error bars on all impurity Green’s functions are negligible
the DMFT self-consistency loop is converged perfectly, and
the fermionic box sizes are sufficiently large that extrapolation can be done.

In practice, given a finite amount of computational resources for the Monte Carlo, the two results are expected to be close.

plot_components.py

Larmor precession

In the presence of a Zeeman field \(B\), the total spin of the system will undergo Larmor precession. This is visible in the in-plane (x,y) components of the spin susceptibility at \(\mathbf{q}=\Gamma\). The frequency of the oscillation is fixed by \(B\) and is thus known analytically. This is a good check on the numerical results.

plot_chiG.py

From left to right, this plot shows the \(S^x S^x\), \(S^x S^y\), \(S^z S^z\) and \(N N\) components of \(\chi(q=\Gamma,\omega_m)\) as a function of \(\omega_m\). The orange lines are the BSE result, with darker lines indicating a larger fermionic frequency box. The blue lines are the DBSE results for the same frequency box sizes, which show essentially perfect convergence. In the first two plots, the black lines are the analytically known exact result. For more information, see Commun. Phys. 6, 289 (2023).