## CTQMC and the S mode¶

When S mode is used in CTQMC solver, the self-energy is computed from the two particle vertex function. The derivation of this approach is available in seminal paper by Baym Kadanoff (PRB 124, 287 (1961) ) Equation 30.

This trick was used in Bulla's work on single impurity model within NRG (arXiv 9804224), hence it was named Bulla's trick.

In this ctqmc, it is available for a general impurity problem, and in the new version works for both the Ising as well as for Full (rotationally invariant) Coulomb repuslion.

\begin{eqnarray} H= \varepsilon_{\alpha\beta} \psi^\dagger_\alpha\psi_\beta + \frac{1}{2}U_{\alpha\beta\gamma\delta}\psi^\dagger_\alpha\psi^\dagger_\beta\psi_\gamma\psi_\delta + \varepsilon_k c_k^\dagger c_k + V_{k\alpha} c^\dagger_k \psi_\alpha + V_{\alpha k}\psi^\dagger_\alpha c_k \end{eqnarray}

and the Green's function

\begin{eqnarray} G_{\alpha\beta}(\tau'-\tau)=-\langle T_\tau \psi_\alpha(\tau') \psi_\beta^\dagger(\tau)\rangle \end{eqnarray}

We use Equation of motion to generate the two-time two particle vertex function, i.e.,

\begin{eqnarray} \frac{d}{d\tau}G(\tau'-\tau)=\delta(\tau-\tau')\delta_{\alpha\beta} -\langle T_\tau \psi_\alpha(\tau')\frac{d}{d\tau}\psi^\dagger_\beta(\tau)\rangle =\delta(\tau-\tau')\delta_{\alpha\beta} -\langle T_\tau \psi_\alpha(\tau')[H,\psi^\dagger_\beta(\tau)]\rangle \end{eqnarray}

Where we use the fact that

\begin{eqnarray} \frac{d}{d\tau}\psi^\dagger(\tau) = \frac{d}{d\tau}e^{H \tau}\psi^\dagger e^{-H \tau}=[H,\psi^\dagger(\tau)] \end{eqnarray}

The communtator can be calculated:

\begin{eqnarray} [H,\psi^\dagger_\beta] = \psi^\dagger_i \varepsilon_{i\beta} + c^\dagger_k V_{k\beta}+ \frac{1}{2}(U_{ijk\beta}-U_{ij\beta k})\psi^\dagger_i\psi^\dagger_j \psi_k \end{eqnarray}

and therefore

\begin{eqnarray} \frac{d}{d\tau}G(\tau'-\tau)=\delta(\tau-\tau')\delta_{\alpha\beta} -\langle T_\tau \psi_\alpha(\tau') \psi^\dagger_i(\tau)\rangle \varepsilon_{i\beta} -\langle T_\tau \psi_\alpha(\tau') c^\dagger_k(\tau)\rangle V_{k\beta} -\frac{1}{2}(U_{ijk\beta}-U_{ij\beta k})\langle T_\tau \psi_\alpha(\tau')\psi^\dagger_i(\tau)\psi^\dagger_j(\tau) \psi_k(\tau^-) \rangle \end{eqnarray}

which is equivalent to

\begin{eqnarray} \frac{d}{d\tau}G(\tau'-\tau)=\delta(\tau-\tau')\delta_{\alpha\beta} +G_{\alpha i}(\tau'-\tau)\varepsilon_{i\beta} +G_{\alpha i}V_{i k}(\frac{d}{d\tau}-\varepsilon_k)^{-1}V_{k \beta} -\frac{1}{2}(U_{ijk\beta}-U_{ij\beta k})\langle T_\tau \psi_\alpha(\tau')\psi^\dagger_i(\tau)\psi^\dagger_j(\tau^-) \psi_k(\tau^-) \rangle \end{eqnarray}

and in imaginary time becomes

\begin{eqnarray} G_{\alpha i}(\omega_n)(i\omega_n\delta_{i\beta}-\varepsilon_{i\beta}-\Delta_{i\beta})=\delta_{\alpha\beta} -\frac{1}{2}(U_{ijk\beta}-U_{ij\beta k}) \int_0^\beta e^{i\omega(\tau'-\tau)}\langle T_\tau \psi_\alpha(\tau')\psi^\dagger_i(\tau)\psi^\dagger_j(\tau) \psi_k(\tau^-) \rangle \end{eqnarray}

Or in matrix form:

\begin{eqnarray} (G(i\omega_n-\varepsilon-\Delta)-1)_{\alpha\beta}= -\frac{1}{2}(U_{ijk\beta}-U_{ij\beta k}) \int_0^\beta e^{i\omega(\tau'-\tau)}\langle T_\tau \psi_\alpha(\tau')\psi^\dagger_i(\tau)\psi^\dagger_j(\tau) \psi_k(\tau^-) \rangle \end{eqnarray}\begin{eqnarray} (G\Sigma)_{\alpha\beta}= -\frac{1}{2}(U_{ijk\beta}-U_{ij\beta k}) \int_0^\beta e^{i\omega(\tau'-\tau)}\langle T_\tau \psi_\alpha(\tau')\psi^\dagger_i(\tau)\psi^\dagger_j(\tau) \psi_k(\tau^-) \rangle \end{eqnarray}

We can also write :

\begin{eqnarray} (G\Sigma)_{\alpha\beta}= -\int_0^\beta e^{i\omega(\tau'-\tau)}\langle T_\tau \psi_\alpha(\tau')\psi^\dagger_i(\tau)\psi^\dagger_j(\tau) \psi_k(\tau^-) \rangle \frac{1}{2}(U_{ijk\beta}-U_{ji k\beta}) \end{eqnarray}

We could replace $1/2(U_{ijk\beta}-U_{ji k\beta})$ with simply $U_{ijk\beta}$, because the correlation function in the $\langle\rangle$ is odd with respect to interchange of $i,j$. However, we want to keep this symmetry in the equation, so that this sum over $i,j$ will symmetrize the result, even if we do not symmetrize it when computing the correlation function.

We see that $\alpha$ and $\beta$ correspond to the same block of orbitals. One of the $i$ and $j$ must be from the same block. Withouth loss of generality, we can decide that $i$ is in the block of $\alpha$ and $\beta$ (if it was $j$, we would just exchange $i$ and $j$ as they appear at the same time. But then there is a sign change in $U$, which compensates for the exchange of $i$ and $j$).

Indices $j$ and $k$ must also be from the same block of orbitals, but can come from any block, either the same as $\alpha$, $\beta$ or different.

We enumerate a block of orbitals with index $ifl$, and for this orbital block, we enumerate all times of creation operators with successive index $is$. When we have $Ns$ creation operators $\psi^\dagger$ from bath $ifl$, $is$ would go from 1 to $Ns$. Hence, index $ifl$ and $is$ uniquely determine a creation operator in the time evolution.

We use combined index $ii$ for two baths, coresponding to baths $(j,k)$ in the above equation.

If we want to store the number operator $\psi^\dagger_j \psi_k$, we need an array of type $Njc[ifl](is,ii)$, which correspond to:

$Njc[\tau,ii=(j,k)]= \langle ... |\psi^\dagger_j(\tau)\psi_k(\tau^-)|... \rangle$.

Separately we store the current approximation for the Green's function, in which we store index to the creation operator $is$, i.e.,

$G(\tau'-\tau)[\tau=(is,ifl),(\alpha,i)] = \langle ...| \psi_\alpha(\tau') \psi_i^\dagger(\tau) |... \rangle$

Once $G(\tau'-\tau)[\tau,(\alpha,i)]$ is stored and $Njc[\tau,(j,k)]$ is ready, we can just multiply the two, and sum over all times $\tau$, to get $G\Sigma$.

$$F_{\alpha,i,j,k}(\tau'-\tau)= \sum_\tau G(\tau'-\tau)[\tau,(\alpha,i)]Njc[\tau,(j,k)]$$

Finally, at the end we multply with tensor $U$ to get

\begin{eqnarray} (G\Sigma)_{\alpha\beta} = -\sum_{i,j,k}\frac{1}{2}(U_{ijk\beta}-U_{ij\beta k})F_{\alpha,i,j,k} \end{eqnarray}
In [ ]: