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.
We start with the general AIM Hamiltonian of the form:
\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$.
\begin{equation} F_{\alpha,i,j,k}(\tau'-\tau)= \sum_\tau G(\tau'-\tau)[\tau,(\alpha,i)]Njc[\tau,(j,k)] \end{equation}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}