Derivation and implementation notes

The CTINT algorithm

The following section describes the CTINT algorithm and its derivation. We first present it for the case of static density-density interactions, and present the extension to general dynamical interactions thereafter. For notations and conventions used throughout these notes, the reader is referred to the attached document (which we will refer to as “vertex conventions”).

Action with static density-density interaction

Let us begin by considering the example of an impurity system with static density-density interaction. The action can then be written as

\[\mathcal{S}=\mathcal{S}_0 + S_{\mathrm{int}}\]

with

\[\mathcal{S}_0 = -\bar{c}_{\bar{\alpha}} \mathcal{G}_{\bar{\alpha}\beta}^{-1} c_\beta\]

and

\[\mathcal{S}_{\mathrm{int}} = \sum_{ab} U_{ab} \int_\tau n_a(\tau) n_b(\tau).\]

Interaction expansion

In CT-INT we proceed by expanding the partition function

(1)\[Z\equiv\int\mathcal{D}[\bar{c},c]e^{\bar{c}_{\bar{\alpha}}\mathcal{G}_{\bar{\alpha}\beta}^{-1}c_\beta-\mathcal{S}_{\mathrm{int}}}\]

in powers of \(\mathcal{S}_{\mathrm{int}}\). One obtains

(2)\[Z=Z_0\sum_{k=0}^\infty\frac{(-1)^k}{k!}\sum_{\mathbf{a},\mathbf{b}}\left(\prod_{i=1}^k U_{a_i b_i}\right)\int_{\boldsymbol{\tau}}\langle n_{a_1}(\tau_1)\dots n_{a_k}(\tau_k)n_{b_1}(\tau_1)\dots n_{b_k}(\tau_k)\rangle_0,\]

where here, and in the following, bold symbols denote vectors. Note that in practice one typically orders the integral times, which allows to cancel the \(1/k!\) prefactor. Using the Wick theorem we can write the non-interacting expectation value in the equation above compactly using a determinant

\[Z=Z_0\sum_{k=0}^\infty\frac{(-1)^k}{k!}\sum_{\mathbf{a},\mathbf{b}}\left(\prod_{i=1}^k U_{a_i b_i}\right)\int_{\boldsymbol{\tau}}\det\hat{\mathcal{G}}_k,\]

Here we have introduced the \(2k\times2k\) matrix

\[[\hat{\mathcal{G}}_k]_{ij}=\mathcal{G}_{x_i\bar{y}_j}.\]

with vectors

\[\begin{split}\boldsymbol{x}[\mathbf{a},\mathbf{b},\boldsymbol{\tau}] & \equiv[(c_1,\tau_1),\ldots,(c_{2k},\tau_{2k})]\equiv[(a_1,\tau_1),\ldots,(a_k,\tau_k),(b_1,\tau_1),\ldots,(b_k,\tau_k)],\\ \bar{\boldsymbol{y}}[\mathbf{a},\mathbf{b},\boldsymbol{\tau}] & \equiv[(\bar{c}_1,\tau_1),\ldots,(\bar{c}_{2k},\tau_{2k})]\equiv[(a_1,\tau_1),\ldots,(a_k,\tau_k),(b_1,\tau_1),\ldots,(b_k,\tau_k)].\end{split}\]

In summary we can write

(3)\[Z = Z_0\Xi\]

with

(4)\[\Xi\equiv\sum_{k=0}^\infty\sum_{\mathbf{a},\mathbf{b}}\int_{\boldsymbol{\tau}}A_k[\mathbf{a},\mathbf{b}]\times\det\hat{\mathcal{G}}_k[\mathbf{a},\mathbf{b},\boldsymbol{\tau}]\]

and

\[A_k[\mathbf{a},\mathbf{b}]=\frac{(-1)^k}{k!}\left(\prod_{i=1}^k U_{a_i b_i}\right)\]

Note that \(A_k\) is dependent on the interaction \(U\), while the determinant is dependent on the bare propagator \(\mathcal{G}\).

The \(\alpha\) function and the auxiliary spin

In order to handle the trivial sign problem originating from the minus-signs associated with the expansion of each static interaction vertex we shift the quartic part of the action \(\mathcal{S}_{\mathrm{int}}\rightarrow\mathcal{\tilde{S}}_{\mathrm{int}}\) with

\[\begin{split}\tilde{\mathcal{S}}_{\mathrm{int}} &= \sum_l U_l \int_\tau \bar{c}_{\bar{a}_l}(\tau) c_{b_l}(\tau) \bar{c}_{\bar{c}_l}(\tau) c_{d_l}(\tau) - n_{\bar{a}_l b_l}(\tau) \alpha_{11}^l - n_{\bar{c}_l d_l}(\tau) \alpha_{00}^l + n_{\bar{a}_l d_l}(\tau) \alpha_{10}^l + n_{\bar{c}_l b_l}(\tau) \alpha_{01}^l \\ &= \sum_l U_l \int_\tau \left( n_{\bar{a}_l b_l}(\tau) - \alpha_{00}^l \right) \left( n_{\bar{c}_l d_l}(\tau) - \alpha_{11}^l \right) + n_{\bar{a}_l d_l}(\tau) \alpha_{10}^l + n_{\bar{c}_l b_l}(\tau) \alpha_{01}^l + \mathrm{const} \\ &= \sum_l U_l \int_\tau - \left( n_{\bar{a}_l d_l}(\tau) - \alpha_{01}^l \right) \left( n_{\bar{c}_l b_l}(\tau) - \alpha_{10}^l \right) - n_{\bar{a}_l b_l}(\tau) \alpha_{11}^l - n_{\bar{c}_l d_l}(\tau) \alpha_{00}^l + \mathrm{const}\end{split}\]

This quadratic shift has to be accounted for in both the kinetic and the interaction term of the action such that the overall action remains unchanged, i.e.

\[\mathcal{S}=\tilde{\mathcal{S}}_0+\tilde{\mathcal{S}}_{\mathrm{int}},\]

with

\[\tilde{\mathcal{S}}_0=-\bar{c}_{\bar{\alpha}}\tilde{\mathcal{G}}_{\bar{\alpha}\beta}^{-1}c_\beta.\]

Here we introduced the non-interacting propagator by the quadratic cross terms of the product above

\[\tilde{\mathcal{G}}_{\bar{e}f}^{-1}(i\omega)\equiv\mathcal{G}_{\bar{e}f}^{-1}(i\omega)-\sum_l U_l\left(\delta_{\bar{e},\bar{a}_l}\delta_{f,b_l}\alpha_{11}^l+\delta_{\bar{e},\bar{c}_l}\delta_{f,d_l}\alpha_{00}^l-\delta_{\bar{e},\bar{a}_l}\delta_{f,d_l}\alpha_{10}^l-\delta_{\bar{e}\bar{c}_l}\delta_{fb_l}\alpha_{01}^l\right)\]

The path-integral of \(Z\) has to be adjusted accordingly

\[Z=Z_0\sum_{k=0}^\infty\frac{(-1)^k}{k!}\sum_{\mathbf{l}}\left(\prod_{i=1}^k U_{l_i}\right)\int_{\boldsymbol{\tau}}\det\hat{\mathcal{G}}_{\mathbf{l}}[\boldsymbol{\tau}],\]

with matrix

(5)\[\begin{split}\hat{\mathcal{G}}_{\mathbf{l}}[\boldsymbol{\tau}]=\begin{pmatrix}\tilde{G}_{\mathbf{b_l\bar{a}_l}}[\boldsymbol{\tau}]-\operatorname{diag}(\boldsymbol{\alpha_{00}^{\mathbf{l}}}) & & \tilde{G}_{\mathbf{b_l\bar{c}_l}}[\boldsymbol{\tau}]-\operatorname{diag}(\boldsymbol{\alpha}_{\mathbf{10}}^{\mathbf{l}})\\ \tilde{G}_{\mathbf{d_l\bar{a}_l}}[\boldsymbol{\tau}]-\operatorname{diag}(\boldsymbol{\alpha}_{\mathbf{01}}^{\mathbf{l}}) & & \tilde{G}_{\mathbf{d_l\bar{c}_l}}[\boldsymbol{\tau}]-\operatorname{diag}(\boldsymbol{\alpha_{11}^l}) \end{pmatrix}\end{split}\]

To see this, one has to rewrite an expectation value of the kind \(\langle\tilde{n}_{a_{l_1}}(\tau_1)\dots\tilde{n}_{b_{l_k}}(\tau_k)\rangle\) into a determinant. This procedure is again based on the Wick theorem, and leads to the matrix in Eq. (5), as can be easily proven, starting from the expression \(\langle n_{a_{l_k}}(\tau_1)\dots n_{b_{l_k}}(\tau_k)\rangle\), by successive replacements \(n\rightarrow\tilde{n}\). We then get

\[\Xi\equiv\sum_{k=0}^\infty\sum_{\mathbf{l}}\int_{\boldsymbol{\tau}}A_{\mathbf{l}}\times\det\hat{\mathcal{G}}_{\mathbf{l}}[\boldsymbol{\tau}]\]

and

\[A_{\mathbf{l}}=\frac{(-1)^k}{k!}\left(\prod_{i=1}^k U_{l_i}\right)\]

We solve the following set of equations self-consistently

\[\begin{split}\tilde{\mathcal{G}}_{b_l\bar{a}_l}[\hat{\alpha}^{\mathbf{l}}] & =\alpha_{00}^l\\ \tilde{\mathcal{G}}_{b_l\bar{c}_l}[\hat{\alpha}^{\mathbf{l}}] & =\alpha_{10}^l\\ \tilde{\mathcal{G}}_{d_l\bar{a}_l}[\hat{\alpha}^{\mathbf{l}}] & =\alpha_{01}^l\\ \tilde{\mathcal{G}}_{d_l\bar{d}_l}[\hat{\alpha}^{\mathbf{l}}] & =\alpha_{11}^l\end{split}\]

in order to determine alpha based off of the self-consistent Hartree Fock solution. We finally introduce a small term-wise assymetry to cope with the trivial part of the sign problem

\[\begin{split}\hat{\alpha}_l\rightarrow\hat{\alpha}_l+\begin{pmatrix}\delta & \delta\\ \delta & -\delta \end{pmatrix}\end{split}\]

General dynamic interactions

Let us now consider an interaction of the most general form

\[\begin{split}\mathcal{S}_{\mathrm{int}} & =\mathop{\sum\int}_{\bar{\alpha}\beta\bar{\gamma}\delta}U_{\bar{\alpha}\beta\bar{\gamma}\delta}\,\bar{c}_\alpha c_\beta \bar{c}_\gamma c_\delta\\ & =\sum_{\bar{a}b\bar{c}d}\int_{\boldsymbol{\tau}}U_{\bar{a}b\bar{c}d}(\tau_a,\tau_b,\tau_c,\tau_d)\times\bar{c}_a(\tau_a)c_b(\tau_b)\bar{c}_c(\tau_c)c_d(\tau_d).\end{split}\]

In this case the expansion of the partition function takes the form

\[Z=Z_0\sum_{k=0}^\infty\frac{(-1)^k}{k!}\mathop{\sum\int}_{\bar{\boldsymbol{\alpha}}\boldsymbol{\beta}\bar{\boldsymbol{\gamma}}\boldsymbol{\delta}}\left(\prod_{i=1}^k U_{\bar{\alpha}_i\beta_i\bar{\gamma}_i\delta_i}\right)\langle\bar{c}_{\alpha_1}c_{\beta_1}\bar{c}_{\gamma_1}c_{\delta_1}\ldots\bar{c}_{\alpha_k}c_{\beta_k}\bar{c}_{\gamma_k}c_{\delta_k}\rangle_0,\]

which again, using the Wick theorem, takes the compact form

\[Z=Z_0\sum_{k=0}^\infty\frac{(-1)^k}{k!}\mathop{\sum\int}_{\bar{\boldsymbol{\alpha}}\boldsymbol{\beta}\bar{\boldsymbol{\gamma}}\boldsymbol{\delta}}\left(\prod_{i=1}^k U_{\bar{\alpha}_i\beta_i\bar{\gamma}_i\delta_i}\right)\det\hat{\mathcal{G}}_k,\]

with a matrix

\[[\hat{\mathcal{G}}_k]_{ij}=\mathcal{G}_{x_i\bar{y}_j},\]

that now has indices

\[\boldsymbol{x}[\boldsymbol{\beta},\boldsymbol{\delta}]\equiv[{\beta_1}{\delta_1}\ldots{\beta_k}{\delta_k}],\qquad\bar{\boldsymbol{y}}[\bar{\boldsymbol{\alpha}},\bar{\boldsymbol{\gamma}}]\equiv[\bar{\alpha}_1\bar{\gamma}_1\ldots\bar{\alpha}_k\bar{\gamma}_k].\]

In summary we find that

(6)\[\Xi=\frac{Z}{Z_0}=\sum_{k=0}^\infty\mathop{\sum\int}_{\bar{\boldsymbol{\alpha}}\boldsymbol{\beta}\bar{\boldsymbol{\gamma}}\boldsymbol{\delta}}A_k[\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma},\boldsymbol{\delta}]\times\det\hat{\mathcal{G}}_k[\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma},\boldsymbol{\delta},\boldsymbol{s}]\]

with

\[A_k[\boldsymbol{\alpha},\boldsymbol{\beta},\boldsymbol{\gamma},\boldsymbol{\delta}]=\frac{(-1)^k}{k!}\left(\prod_{i=1}^k U_{\bar{\alpha}_i\beta_i\bar{\gamma}_i\delta_i}\right).\]

Eq. (6) can be understood as a sum over weights \(w_{\mathcal{C}}\) for all configurations \(\mathcal{C}\)

\[\sum_{\mathcal{C}}w_{\mathcal{C}}\]

Dynamic density-density and spin-exchange interactions

Let us now consider retarded interactions of both the density-density (\(\mathcal{D}\)) and the spin-exchange (\(\mathcal{J}^\perp\)) type. In this case the interacting part of the action reads

\[\begin{split}\mathcal{S}_{\mathrm{int}} &= \frac{1}{N_s}\iint_{\tau,\tau'}\sum_{abs}\mathcal{D}_{ab}(\tau-\tau')\tilde{n}_{as}(\tau)\tilde{n}_{bs}(\tau')\\ & +\frac{1}{2}\iint_{\tau,\tau'}\sum_{uv}\mathcal{J}_{uv}^\perp(\tau-\tau')\left[s_u^+(\tau)s_v^-(\tau')+s_u^-(\tau)s_v^+(\tau')\right]\end{split}\]

with

\[\begin{split}s_u^+(\tau) & =\bar{c}_{u,\uparrow}(\tau)c_{u,\downarrow}(\tau),\\ s_u^-(\tau) & =\bar{c}_{u,\downarrow}(\tau)c_{u,\uparrow}(\tau).\end{split}\]

As done already in Sec. Action with static density-density interaction, we have to adjust the non-interacting propagator due to the use of the shifted densities \(\tilde{n}\)

(7)\[\tilde{\mathcal{G}}_{\bar{a}b}^{-1}(i\omega)\equiv\mathcal{G}_{\bar{a}b}^{-1}(i\omega)-\delta_{\bar{a},b}\sum_c\frac{1}{N_s}\left[\mathcal{D}_{bc}(i\Omega=0)+\mathcal{D}_{cb}(i\Omega=0)\right]\sum_s\alpha_{cs}.\]

Monte-Carlo algorithm

The sum in Eq. (4) is evaluated using the Metropolis-Hastings algorithm. It can be understood as a sum over configurations \(\mathcal{C}\equiv(k,\mathbf{a},\mathbf{b},\boldsymbol{\tau})\) (and possibly \({\bf s}\)) with different probabilities \(|w_{\mathcal{C}}|\), i.e.

\[\Xi=\sum_{\mathcal{C}}w_{\mathcal{C}}=\sum_{\mathcal{C}}|w_{\mathcal{C}}|\operatorname{sign}(\mathcal{C}).\]

The weights \(|w_{\mathcal{C}}|\) are used to construct a Markov chain of configurations satisfying detailed balance and ergodicity. Under these conditions, the expectation value of an observable \(\hat{O}\) can be computed as:

(8)\[\langle\hat{O}\rangle=\frac{\sum_{\mathcal{C}}|w_{\mathcal{C}}|\operatorname{sign}(\mathcal{C})\hat{O}_{\mathcal{C}}}{\sum_{\mathcal{C}}|w_{\mathcal{C}}|\operatorname{sign}(\mathcal{C})}=\frac{\langle\operatorname{sign}\hat{O}\rangle_{\mathrm{MC}}}{\langle\operatorname{sign}\rangle_{\mathrm{MC}}}\]

where we have defined

(9)\[\langle\hat{A}\rangle_{\mathrm{MC}}\equiv\sum_{\mathcal{C}}|w_{\mathcal{C}}|\hat{A}_{\mathcal{C}}.\]

Here, \(\hat{A}_{\mathcal{C}}\) is the value of observable \(\hat{A}\) when evaluated in the MC configuration \(\mathcal{C}\).

Correlation functions and measurements

In this section we discuss the measurement of different correlation functions in CTInt. Note that in the following we will not explicitly denote that the propagator has been shifted using a \(\tilde{\mathcal{G}}\), but we will rather use the symbol \(\mathcal{G}\).

Single-particle Green’s function \(G\)

From Eq. (1) we see that \(\frac{\partial Z}{\partial\mathcal{G}_{\bar{\beta}\alpha}^{-1}}=Z\langle\bar{c}_{\bar{\beta}}c_\alpha\rangle\), i.e.

(10)\[G_{\alpha\bar{\beta}}=\frac{1}{Z}\frac{\partial Z}{\partial\mathcal{G}_{\bar{\beta}\alpha}^{-1}}\]

Using Eqs. (10), (3) and (4) this can be further evaluated to

(11)\[G_{\alpha\bar{\beta}}=\frac{1}{Z}\left[\frac{\partial Z_0}{\partial\mathcal{G}_{\bar{\beta}\alpha}^{-1}}\Xi+Z_0\frac{\partial\Xi}{\partial\mathcal{G}_{\bar{\beta}\alpha}^{-1}}\right]=\frac{1}{Z}\left[\mathcal{G}_{\alpha\bar{\beta}}Z_0\Xi-Z_0\Xi\mathcal{G}_{\alpha\bar{\delta}}\frac{1}{\Xi}\frac{\partial\Xi}{\partial\mathcal{G}_{\gamma\bar{\delta}}}\mathcal{G}_{\gamma\bar{\beta}}\right].\]

Hence

(12)\[G_{\alpha\bar{\beta}}=\mathcal{G}_{\alpha\bar{\beta}}+\mathcal{G}_{\alpha\bar{\delta}}M_{\bar{\delta}\gamma}\mathcal{G}_{\gamma\bar{\beta}},\]

with

(13)\[M_{\bar{\alpha}\beta}\equiv-\frac{1}{\Xi}\frac{\partial\Xi}{\partial\mathcal{G}_{\beta\bar{\alpha}}}.\]

Instead of measuring \(G\) directly, one typically measures the quantity \(M_{\bar{\alpha}\beta}\) instead, as is done in the CTINT code of TRIQS. Let us see how it can be written by means of a Monte-Carlo average. Eq. (4) gives

\[M_{\bar{\alpha}\beta}=-\frac{1}{\Xi}\sum_{k=0}^\infty\sum_{\mathbf{a},\mathbf{b},\mathbf{s}}\int_{\boldsymbol{\tau}}A_k\frac{\partial\det\hat{\mathcal{G}}_k}{\partial\mathcal{G}_{\beta\bar{\alpha}}}.\]

Using the chain rule

\[\frac{\partial\det\hat{\mathcal{G}}_k}{\partial\mathcal{G}_{\beta\bar{\alpha}}}=\sum_{ij}\frac{\partial\det\hat{\mathcal{G}}_k}{\partial[\mathcal{\hat{G}}_k]_{ij}}\frac{\partial[\hat{\mathcal{G}}_k]_{ij}}{\partial\mathcal{G}_{\beta\bar{\alpha}}}=\det\hat{\mathcal{G}}_k\sum_{ij}[\hat{\mathcal{G}}_k^{-1}]_{ji}\delta_{x_i\beta}\delta_{\bar{y}_i\bar{\alpha}}\]

we obtain

(14)\[\begin{split}M_{\bar{\alpha}\beta} & =-\frac{1}{\Xi}\sum_{k=0}^\infty\sum_{\mathbf{a},\mathbf{b},\mathbf{s}}\int_{\boldsymbol{\tau}}A_k\frac{\partial\det\hat{\mathcal{G}}_k}{\partial\mathcal{G}_{\beta\bar{\alpha}}}\nonumber \\ & =-\frac{1}{\Xi}\sum_{\mathcal{C}}w_{\mathcal{C}}\left\{ \sum_{ij}\delta_{\bar{y}_i\bar{\alpha}}\delta_{x_i\beta}[\hat{\mathcal{G}}_k^{-1}]_{ji}\right\} \nonumber \\ & =\frac{\Big\langle-\sum_{ij}\delta_{\bar{y}_i\bar{\alpha}}\delta_{x_i\beta}[\hat{\mathcal{G}}_k^{-1}]_{ji}\,\operatorname{sign}(\mathcal{C})\Big\rangle_{\mathrm{MC}}}{\Big\langle\operatorname{sign}(\mathcal{C})\Big\rangle_{\mathrm{MC}}}.\end{split}\]

Imaginary time measurement

Just like the Green function, \(M_{\bar{\alpha}\beta}\) depends only on the time-difference \(\tau=\tau_{\bar{\alpha}}-\tau_\beta\), which, in practice, is chosen from a discrete set of points \(\{\tau_l\}\). They are determined by binning the interval \([0,\beta)\) equidistantly, such that grid points are seperated by \(\Delta\tau\). We thus define

\[M_{\bar{a}b}^{\Delta\tau}(\tau_l)=-\frac{1}{\beta\langle\operatorname{sign}(\mathcal{C})\rangle_{\mathrm{MC}}}\left\langle \sum_{ij}\delta_{\bar{c}_j\bar{a}}\delta_{c_i b}\left[\delta_{\Delta\tau}(\tau_l-\tau_j+\tau_i)-\delta_{\Delta\tau}(\tau_l-\beta-\tau_j+\tau_i)\right][\hat{\mathcal{G}}_k^{-1}]_{ji}\operatorname{sign}(\mathcal{C})\right\rangle _{\mathrm{MC}}\]

where \(\delta_{\Delta\tau}\) is the broadened delta function

\[\begin{split}\delta_{\Delta\tau}(\tau)\equiv\begin{cases} \frac{1}{\Delta\tau} & 0<\tau<\Delta\tau\\ 0 & \tau\leq0\quad\mathrm{or}\quad\tau\geq\Delta\tau, \end{cases}\end{split}\]

and a factor \(1/\beta\) should be introduced when making use of time-translational invariance \((\tau_{\bar{\alpha}},\tau_\beta)\rightarrow\tau=\tau_{\bar{\alpha}}-\tau_\beta\). The second \(\delta\) function with the minus sign takes care of shifting negative time-differences \(\tau_{\bar{\alpha}}-\tau_\beta\) to the interval \([0,\beta)\) using the antiperiodicity of \(M\). As \(M_{\bar{u}v}(\tau_l)=\lim_{\Delta\tau\rightarrow0}{M}_{\bar{u}v}^{\Delta\tau}(\tau_l)\) we can approximate \(M_{\bar{u}v}(\tau_l)\approx M_{\bar{u}v}^{\Delta\tau}(\tau_l)\) for a sufficiently small \(\Delta\tau\).

Frequency measurement

Alternatively, we can measure \(M\) directly in frequency space, i.e.

\[M_{\bar{a}b}(i\omega_n)=-\frac{1}{\beta\langle\operatorname{sign}(\mathcal{C})\rangle_{\mathrm{MC}}}\left\langle \sum_{ij}\delta_{\bar{c}_j\bar{a}}\delta_{c_i b}e^{i\omega_n(\tau_j-\tau_i)}[\hat{\mathcal{G}}_k^{-1}]_{ji}\operatorname{sign}(\mathcal{C})\right\rangle _{\mathrm{MC}}\]

Here, we can avoid the binning, while a factor \(1/\beta\) has to be introduced as before. To evaluate this measure, we can make use of the NFFT library. To make this more obvious we rewrite

\[e^{i\omega_n(\tau_j-\tau_i)}=e^{i\pi(2n+1)\left(\frac{\tau_j-\tau_i}{\beta}-\frac{1}{2}+\frac{1}{2}\right)}=e^{2\pi inx}e^{i\pi\left(n+x+\frac{1}{2}\right)}\]

with \(x\equiv\left(\frac{\tau_j-\tau_i}{\beta}-\frac{1}{2}\right)\in[-0.5,0.5)\) and \(n\in\mathbb{Z}\), as is required by the interface of the library.

Two-particle Green’s function \(\chi^4\)

The same idea can be pursued when considering two-particle Green functions. From Eq. (1) we can directly see that \(\frac{\partial Z}{\partial\mathcal{G}_{\bar{\alpha}\beta}^{-1}\partial\mathcal{G}_{\bar{\gamma}\delta}^{-1}}=Z\langle\bar{c}_\alpha c_\beta \bar{c}_\gamma c_\delta\rangle\) i.e

(15)\[\chi_{\bar{\alpha}\beta\bar{\gamma}\delta}^4=\langle\bar{c}_\alpha c_\beta \bar{c}_\gamma c_\delta\rangle=\frac{1}{Z}\frac{\partial Z}{\partial\mathcal{G}_{\bar{\alpha}\beta}^{-1}\partial\mathcal{G}_{\bar{\gamma}\delta}^{-1}}\]

Starting from Eqs. (15) and (10) we first note:

(16)\[\chi_{\bar{\alpha}\beta\bar{\eta}\kappa}^4=\frac{1}{Z}\frac{\partial(ZG_{\beta\bar{\alpha}})}{\partial\mathcal{G}_{\bar{\eta}\kappa}^{-1}}=-\frac{1}{Z}\mathcal{G}_{\mu\bar{\eta}}\frac{\partial(ZG_{\beta\bar{\alpha}})}{\partial\mathcal{G}_{\mu\bar{\nu}}}\mathcal{G}_{\kappa\bar{\nu}}\]

Let us now evaluate \(\frac{\partial(ZG_{\beta\bar{\alpha}})}{\partial\mathcal{G}_{\mu\bar{\nu}}}\) using (11). We have

(17)\[\begin{split}\frac{\partial(ZG_{\beta\bar{\alpha}})}{\partial\mathcal{G}_{\mu\bar{\nu}}} & =\frac{\partial}{\partial\mathcal{G}_{\mu\bar{\nu}}}\left[\mathcal{G}_{\beta\bar{\alpha}}Z_0\Xi-Z_0\mathcal{G}_{\gamma\bar{\alpha}}\frac{\partial\Xi}{\partial\mathcal{G}_{\gamma\bar{\delta}}}\mathcal{G}_{\beta\bar{\delta}}\right]\nonumber \\ & =\delta_{\beta\mu}\delta_{\bar{\nu}\bar{\alpha}}Z_0\Xi+\mathcal{G}_{\beta\bar{\alpha}}\frac{\partial Z_0}{\partial\mathcal{G}_{\mu\bar{\nu}}}\Xi+\mathcal{G}_{\beta\bar{\alpha}}Z_0\frac{\partial\Xi}{\partial\mathcal{G}_{\mu\bar{\nu}}}\nonumber \\ & \;\;-\frac{\partial Z_0}{\partial\mathcal{G}_{\mu\bar{\nu}}}\mathcal{G}_{\gamma\bar{\alpha}}\frac{\partial\Xi}{\partial\mathcal{G}_{\gamma\bar{\delta}}}\mathcal{G}_{\beta\bar{\delta}}-Z_0\delta_{\mu\gamma}\delta_{\bar{\alpha}\bar{\nu}}\frac{\partial\Xi}{\partial\mathcal{G}_{\gamma\bar{\delta}}}\mathcal{G}_{\beta\bar{\delta}}-Z_0\mathcal{G}_{\gamma\bar{\alpha}}\frac{\partial\Xi}{\partial\mathcal{G}_{\mu\bar{\nu}}\partial\mathcal{G}_{\gamma\bar{\delta}}}\mathcal{G}_{\beta\bar{\delta}}-Z_0\mathcal{G}_{\gamma\bar{\alpha}}\frac{\partial\Xi}{\partial\mathcal{G}_{\gamma\bar{\delta}}}\delta_{\beta\mu}\delta_{\bar{\delta}\bar{\nu}}\nonumber \\ & =\delta_{\beta\mu}\delta_{\bar{\nu}\bar{\alpha}}Z_0\Xi-\mathcal{G}_{\beta\bar{\alpha}}\mathcal{G}_{\bar{\nu}\mu}^{-1}Z_0\Xi-\mathcal{G}_{\beta\bar{\alpha}}Z_0\Xi M_{\bar{\nu}\mu}\nonumber \\ & \;\;-\mathcal{G}_{\bar{\nu}\mu}^{-1}Z_0\mathcal{G}_{\gamma\bar{\alpha}}\Xi M_{\bar{\delta}\gamma}\mathcal{G}_{\beta\bar{\delta}}+\Xi Z_0\delta_{\mu\gamma}\delta_{\bar{\alpha}\bar{\nu}}M_{\bar{\delta}\gamma}\mathcal{G}_{\beta\bar{\delta}}-Z_0\mathcal{G}_{\gamma\bar{\alpha}}\Xi M_{\bar{\nu}\mu\bar{\delta}\gamma}^4\mathcal{G}_{\beta\bar{\delta}}+Z_0\mathcal{G}_{\gamma\bar{\alpha}}\Xi M_{\bar{\delta}\gamma}\delta_{\beta\mu}\delta_{\bar{\delta}\bar{\nu}}\nonumber \\ & =Z\Big\{\delta_{\beta\mu}\delta_{\bar{\nu}\bar{\alpha}}-\mathcal{G}_{\beta\bar{\alpha}}\mathcal{G}_{\bar{\nu}\mu}^{-1}-\mathcal{G}_{\beta\bar{\alpha}}M_{\bar{\nu}\mu}-\mathcal{G}_{\bar{\nu}\mu}^{-1}\mathcal{G}_{\gamma\bar{\alpha}}M_{\bar{\delta}\gamma}\mathcal{G}_{\beta\bar{\delta}}\nonumber \\ & \;\;\;+\delta_{\bar{\alpha}\bar{\nu}}M_{\bar{\delta}\mu}\mathcal{G}_{\beta\bar{\delta}}-\mathcal{G}_{\gamma\bar{\alpha}}M_{\mu\bar{\nu}\gamma\bar{\delta}}^4\mathcal{G}_{\beta\bar{\delta}}+\mathcal{G}_{\gamma\bar{\alpha}}M_{\bar{\nu}\gamma}\delta_{\beta\mu}\Big\}\end{split}\]

where we have defined

(18)\[M_{\mu\bar{\nu}\gamma\bar{\delta}}^4\equiv\frac{1}{\Xi}\frac{\partial\Xi}{\partial\mathcal{G}_{\mu\bar{\nu}}\partial\mathcal{G}_{\gamma\bar{\delta}}}.\]

Plugging (17) into (16), we obtain

\[\begin{split}\chi_{\bar{\alpha}\beta\bar{\eta}\kappa}^4 & =-\mathcal{G}_{\beta\bar{\eta}}\mathcal{G}_{\kappa\bar{\alpha}}+\mathcal{G}_{\kappa\bar{\eta}}\mathcal{G}_{\beta\bar{\alpha}}+\mathcal{G}_{\beta\bar{\alpha}}\mathcal{G}_{\kappa\bar{\nu}}M_{\bar{\nu}\mu}\mathcal{G}_{\mu\bar{\eta}}+\mathcal{G}_{\beta\bar{\delta}}M_{\bar{\delta}\gamma}\mathcal{G}_{\gamma\bar{\alpha}}\mathcal{G}_{\kappa\bar{\eta}}\\ & \;\;-\mathcal{G}_{\beta\bar{\delta}}M_{\bar{\delta}\mu}\mathcal{G}_{\mu\bar{\eta}}\mathcal{G}_{\kappa\bar{\alpha}}+\mathcal{G}_{\mu\bar{\alpha}}\mathcal{G}_{\beta\bar{\nu}}M_{\mu\bar{\nu}\gamma\bar{\delta}}^4\mathcal{G}_{\gamma\bar{\eta}}\mathcal{G}_{\kappa\bar{\delta}}-\mathcal{G}_{\beta\bar{\eta}}\mathcal{G}_{\kappa\bar{\nu}}M_{\bar{\nu}\gamma}\mathcal{G}_{\gamma\bar{\alpha}}\end{split}\]

i.e, after relabeling and reordering,

(19)\[\begin{split}\chi_{\bar{\alpha}\beta\bar{\gamma}\delta}^4 & =\mathcal{G}_{\mu\bar{\alpha}}\mathcal{G}_{\beta\bar{\nu}}M_{\mu\bar{\nu}\kappa\bar{\eta}}^4\mathcal{G}_{\kappa\bar{\gamma}}\mathcal{G}_{\delta\bar{\eta}}\nonumber \\ & \;\;+\mathcal{G}_{\beta\bar{\alpha}}\mathcal{G}_{\delta\bar{\gamma}}+[\mathcal{G}M\mathcal{G}]_{\beta\bar{\alpha}}\mathcal{G}_{\delta\bar{\gamma}}+\mathcal{G}_{\beta\bar{\alpha}}[\mathcal{G}M\mathcal{G}]_{\delta\bar{\gamma}}\nonumber \\ & \;\;-\mathcal{G}_{\beta\bar{\gamma}}\mathcal{G}_{\delta\bar{\alpha}}-[\mathcal{G}M\mathcal{G}]_{\beta\bar{\gamma}}\mathcal{G}_{\delta\bar{\alpha}}-\mathcal{G}_{\beta\bar{\gamma}}[\mathcal{G}M\mathcal{G}]_{\delta\bar{\alpha}}.\end{split}\]

Thus, in order to compute \(\chi^4\), one needs to compute \(M\) and \(M^4\).

Connected four-point function

By inspection of Eq.:eq:chi4_final we see that the connected part of the two-particle Green function now reads

(20)\[\chi_{\bar{\alpha}\beta\bar{\gamma}\delta}^{4,\mathrm{conn}}=\mathcal{G}_{\mu\bar{\alpha}}\mathcal{G}_{\beta\bar{\nu}}M_{\mu\bar{\nu}\kappa\bar{\eta}}^4\mathcal{G}_{\kappa\bar{\gamma}}\mathcal{G}_{\delta\bar{\eta}}-\left[\mathcal{G}M\mathcal{G}\right]_{\alpha\bar{\beta}}\left[\mathcal{G}M\mathcal{G}\right]_{\gamma\bar{\delta}}+\left[\mathcal{G}M\mathcal{G}\right]_{\alpha\bar{\delta}}\left[\mathcal{G}M\mathcal{G}\right]_{\gamma\bar{\beta}}.\]

Instead, we can write

\[\chi_{\bar{\alpha}\beta\bar{\gamma}\delta}^{4,\mathrm{conn}}=\mathcal{G}_{\mu\bar{\alpha}}\mathcal{G}_{\beta\bar{\nu}}M_{\mu\bar{\nu}\kappa\bar{\eta}}^{4,\mathrm{conn}}\mathcal{G}_{\kappa\bar{\gamma}}\mathcal{G}_{\delta\bar{\eta}}\]

with

\[M_{\mu\bar{\nu}\kappa\bar{\eta}}^{4,\mathrm{conn}}\equiv M_{\mu\bar{\nu}\kappa\bar{\eta}}^4-M_{\bar{\nu}\mu}M_{\bar{\eta}\kappa}+M_{\bar{\nu}\kappa}M_{\bar{\eta}\mu}.\]

Fully reducible vertex

Using the definition of the fully reducible vertex function \(F\) (see vertex conventions), we get

\[F_{\alpha\bar{\beta}\gamma\bar{\delta}}=\bigl\{ G_{\bar{\beta}\beta}^{-1}\mathcal{G}_{\beta\bar{\nu}}\bigr\}\left\{ G_{\bar{\delta}\delta}^{-1}\mathcal{G}_{\delta\bar{\eta}}\right\} M_{\mu\bar{\nu}\kappa\bar{\eta}}^{4,\mathrm{conn}}\left\{ \mathcal{G}_{\mu\bar{\alpha}}G_{\bar{\alpha}\alpha}^{-1}\right\} \left\{ \mathcal{G}_{\kappa\bar{\gamma}}G_{\bar{\gamma}\gamma}^{-1}\right\} .\]

With Eq. (12) we can calculate \(F\) directly from \(M^4\) and \(M\)

\[F_{\alpha\bar{\beta}\gamma\bar{\delta}}=\left\{ 1+M\mathcal{G}\right\} _{\bar{\beta}\bar{\nu}}^{-1}\left\{ 1+M\mathcal{G}\right\} _{\bar{\delta}\bar{\eta}}^{-1}M_{\mu\bar{\nu}\kappa\bar{\eta}}^{4,\mathrm{conn}}\left\{ 1+M\mathcal{G}\right\} _{\mu\alpha}^{-1}\left\{ 1+M\mathcal{G}\right\} _{\kappa\gamma}^{-1}.\]

From this expression, one can see that \(F\) and \(M^{4,\mathrm{conn}}\) have the same asymptotics.

Imaginary time measurement of \(M^4\)

In order to write \(M^4\) by means of a Monte-Carlo average we have to evaluate

\[M_{\alpha\bar{\beta}\gamma\bar{\delta}}^4=\frac{1}{\Xi}\frac{\partial\Xi}{\partial\mathcal{G}_{\alpha\bar{\beta}}\partial\mathcal{G}_{\gamma\bar{\delta}}}=-\frac{1}{\Xi}\sum_{k=0}^\infty\sum_{\mathbf{a},\mathbf{b},\mathbf{s}}\int_{\boldsymbol{\tau}}A_k\,\frac{\partial\det\hat{\mathcal{G}}_k}{\partial\mathcal{G}_{\alpha\bar{\beta}}\partial\mathcal{G}_{\gamma\bar{\delta}}}.\]

This derivative evaluates to (neglecting the \(k\) index for now)

\[\begin{split}\frac{\partial\det\hat{\mathcal{G}}}{\partial\mathcal{G}_{\alpha\bar{\beta}}\partial\mathcal{G}_{\gamma\bar{\delta}}} & =\frac{\partial}{\partial\mathcal{G}_{\alpha\bar{\beta}}}\det\hat{\mathcal{G}}\,\sum_{kl}[\hat{\mathcal{G}}^{-1}]_{lk}\delta_{x_k\gamma}\delta_{\bar{y}_l\bar{\delta}}\\ & =\det\hat{\mathcal{G}}\,\sum_{ijkl}\left\{ [\hat{\mathcal{G}}^{-1}]_{ji}[\hat{\mathcal{G}}^{-1}]_{lk}-[\hat{\mathcal{G}}^{-1}]_{li}[\hat{\mathcal{G}}^{-1}]_{jk}\right\} \delta_{x_i\alpha}\delta_{\bar{y}_i\bar{\beta}}\delta_{x_k\gamma}\delta_{\bar{y}_l\bar{\delta}},\end{split}\]

where we had to use that

\[\phantom{=}\frac{\partial}{\partial\mathcal{G}_{\alpha\bar{\beta}}}\sum_{kl}[\hat{\mathcal{G}}^{-1}]_{lk}\delta_{x_k\gamma}\delta_{\bar{y}_l\bar{\delta}}=-\sum_{ijkl}[\hat{\mathcal{G}}^{-1}]_{li}\underbrace{\frac{\partial\hat{\mathcal{G}}_{ij}}{\partial\mathcal{G}_{\alpha\bar{\beta}}}}_{\delta_{x_i\alpha}\delta_{\bar{y}_i\bar{\beta}}}[\hat{\mathcal{G}}^{-1}]_{jk}\delta_{x_k\gamma}\delta_{\bar{y}_l\bar{\delta}}.\]

The Monte-Carlo average for \(M^4\) thus takes the form

(21)\[M_{\alpha\bar{\beta}\gamma\bar{\delta}}^4=\frac{1}{\langle\operatorname{sign}(\mathcal{C})\rangle_{\mathrm{MC}}}\Bigg\langle\sum_{ijkl}\left\{ [\hat{\mathcal{G}}^{-1}]_{ji}[\hat{\mathcal{G}}^{-1}]_{lk}-[\hat{\mathcal{G}}^{-1}]_{li}[\hat{\mathcal{G}}^{-1}]_{jk}\right\} \delta_{x_i\alpha}\delta_{\bar{y}_j\bar{\beta}}\delta_{x_k\gamma}\delta_{\bar{y}_l\bar{\delta}}\operatorname{sign}(\mathcal{C})\Bigg\rangle_{\mathrm{MC}}.\]

Using the intermediate scattering matrix

\[\bar{M}_{\bar{\beta}\alpha}\equiv\sum_{ij}[\hat{\mathcal{G}}^{-1}]_{ji}\delta_{x_i\alpha}\delta_{\bar{y}_j\bar{\beta}}\]

it takes the compact form

\[M_{\alpha\bar{\beta}\gamma\bar{\delta}}^4=\frac{1}{\langle\operatorname{sign}(\mathcal{C})\rangle_{\mathrm{MC}}}\Big\langle\left(\bar{M}_{\bar{\beta}\alpha}\bar{M}_{\bar{\delta}\gamma}-\bar{M}_{\bar{\delta}\alpha}\bar{M}_{\bar{\beta}\gamma}\right)\operatorname{sign}(\mathcal{C})\Big\rangle_{\mathrm{MC}}.\]

This factorization is directly used in our implementation.

Frequency measurement of \(M^4\)

Just like \(M\), we can measure \(M^4\) also directly in Matsubara frequencies [1]

\[\begin{split}\begin{split}M_{a\bar{b}c\bar{d}}^4(i\omega_a,i\omega_{\bar{b}},i\omega_c) & =\frac{1}{\beta\langle\operatorname{sign}(\mathcal{C})\rangle_{\mathrm{MC}}}\Bigg\langle\left(\bar{M}_{\bar{\beta}\alpha}\bar{M}_{\bar{\delta}\gamma}-\bar{M}_{\bar{\delta}\alpha}\bar{M}_{\bar{\beta}\gamma}\right)\\ & \hspace{2cm}\times e^{-i\omega_a(\tau_a-\tau_d)}\,e^{i\omega_{\bar{b}}(\tau_b-\tau_d)}\,e^{-i\omega_c(\tau_c-\tau_d)}\operatorname{sign}(\mathcal{C})\Bigg\rangle_{\mathrm{MC}}\\ & =\frac{1}{\beta\langle\operatorname{sign}(\mathcal{C})\rangle_{\mathrm{MC}}}\Bigg\langle\Big[\bar{M}_{\bar{b}a}(\omega_b,\omega_a)\bar{M}_{\bar{d}c}(\omega_a+\omega_c-\omega_b,\omega_c)\\ & \hspace{2cm}-\bar{M}_{\bar{d}a}(\underbrace{\omega_a+\omega_c-\omega_b}_{\omega_d},\omega_a)\bar{M}_{\bar{b}c}(\omega_b,\omega_c)\Big]\operatorname{sign}(\mathcal{C})\Bigg\rangle_{\mathrm{MC}}. \end{split}\end{split}\]

Here we have defined the intermediate scattering matrix in the Matsubara domain

\[\bar{M}_{\bar{b}a}(\omega_b,\omega_a)\equiv\sum_{ij}[\hat{\mathcal{G}}^{-1}]_{ji}\delta_{c_i a}\delta_{\bar{c}_j\bar{b}}\times\,e^{i\omega_b\tau_j}e^{-i\omega_a\tau_i}.\]

In the implementation we precompute this intermediate scattering matrix, which allows us to reduce the effort from \(\mathcal{O}(k^4\log(k))\) down to \(\mathcal{O}(k^2\log(k))\).

Equal-time correlation functions \(\chi^3\)

The equal-time correlation functions \(\chi^{3,r}\) can be directly calculated from the intermediate quantity \(M^{3,r}\) defined in the following. In the first step we define three extensions of the intermediate scattering matrix

\[\overline{\mathcal{G}M}_{b\alpha}\equiv\overline{\mathcal{G}M}_{ba}(\tau_a)\equiv\sum_{ij}\mathcal{G}_{b\bar{c}_j}(-\tau_j)[\hat{\mathcal{G}}^{-1}]_{ji}\delta_{x_i\alpha}=-\sum_{ij}\mathcal{G}_{b\bar{c}_j}(\beta-\tau_j)[\hat{\mathcal{G}}^{-1}]_{ji}\delta_{x_i\alpha},\]

\[\overline{M\mathcal{G}}_{\bar{\beta}\bar{a}}\equiv\overline{M\mathcal{G}}_{\bar{b}\bar{a}}(\tau_b)\equiv\sum_{ij}\delta_{y_j\beta}[\hat{\mathcal{G}}^{-1}]_{ji}\mathcal{G}_{c_i\bar{a}}(\tau_i),\]

\[\overline{\mathcal{G}M\mathcal{G}}_{b\bar{a}}\equiv\sum_{ij}\mathcal{G}_{b\bar{c}_j}(-\tau_j)[\hat{\mathcal{G}}^{-1}]_{ji}\mathcal{G}_{c_i\bar{a}}(\tau_i)=-\sum_{ij}\mathcal{G}_{b\bar{c}_j}(\beta-\tau_j)[\hat{\mathcal{G}}^{-1}]_{ji}\mathcal{G}_{c_i\bar{a}}(\tau_i).\]

Using these, \(M^{3,r}\) is measured as

\[M_{abcd}^{3,pp}(\tau_a,\tau_c)=\frac{1}{\langle\operatorname{sign}(\mathcal{C})\rangle_{\mathrm{MC}}}\Bigg\langle\Big[\overline{\mathcal{G}M}_{ba}(\tau_a)\overline{\mathcal{G}M}_{dc}(\tau_c)-\overline{\mathcal{G}M}_{da}(\tau_a)\overline{\mathcal{G}M}_{bc}(\tau_c)\Big]\operatorname{sign}(\mathcal{C})\Bigg\rangle_{\mathrm{MC}},\]

\[M_{a\bar{b}\bar{c}d}^{3,ph}(\tau_a,\tau_{\bar{b}})=\frac{1}{\langle\operatorname{sign}(\mathcal{C})\rangle_{\mathrm{MC}}}\Bigg\langle\Big[\bar{M}_{\bar{\beta}\alpha}\overline{\mathcal{G}M\mathcal{G}}_{d\bar{c}}-\overline{\mathcal{G}M}_{d\alpha}\overline{M\mathcal{G}}_{\bar{\beta}\bar{c}}\Big]\operatorname{sign}(\mathcal{C})\Bigg\rangle_{\mathrm{MC}}.\]

Frequency measurement of \(M^{3,r}\)

Again, we can measure \(M^{3,r}\) also directly in Matsubara frequencies. For this we need the scattering matrices

\[\overline{\mathcal{G}M}_{ba}(\omega_a)\equiv\sum_{ij}\mathcal{G}_{b\bar{c}_j}(-\tau_j)[\hat{\mathcal{G}}^{-1}]_{ji}e^{-i\omega_a\tau_i}\delta_{c_i a}=\sum_{ij}\mathcal{G}_{b\bar{c}_j}(\beta-\tau_j)[\hat{\mathcal{G}}^{-1}]_{ji}e^{i\omega_a(\beta-\tau_i)}\delta_{c_i a},\]

\[\overline{M\mathcal{G}}_{\bar{b}\bar{a}}(\omega_{\bar{b}})\equiv\sum_{ij}e^{i\omega_{\bar{b}}\tau_j}\delta_{\bar{c}_j\bar{b}}[\hat{\mathcal{G}}^{-1}]_{ji}\mathcal{G}_{c_i\bar{a}}(\tau_i),\]

\[\overline{\mathcal{G}M\mathcal{G}}_{b\bar{a}}\equiv\sum_{ij}\mathcal{G}_{b\bar{c}_j}(-\tau_j)[\hat{\mathcal{G}}^{-1}]_{ji}\mathcal{G}_{c_i\bar{a}}(\tau_i)=-\sum_{ij}\mathcal{G}_{b\bar{c}_j}(\beta-\tau_j)[\hat{\mathcal{G}}^{-1}]_{ji}\mathcal{G}_{c_i\bar{a}}(\tau_i),\]

Then we can measure

\[M_{abcd}^{3,pp}(\omega_a,\omega_c)=\frac{1}{\langle\operatorname{sign}(\mathcal{C})\rangle_{\mathrm{MC}}}\Bigg\langle\Big[\overline{\mathcal{G}M}_{ba}(\omega_a)\overline{\mathcal{G}M}_{dc}(\omega_c)-\overline{\mathcal{G}M}_{da}(\omega_a)\overline{\mathcal{G}M}_{bc}(\omega_c)\Big]\operatorname{sign}(\mathcal{C})\Bigg\rangle_{\mathrm{MC}},\]

\[M_{a\bar{b}\bar{c}d}^{3,ph}(\omega_a,\omega_{\bar{b}})=\frac{1}{\langle\operatorname{sign}(\mathcal{C})\rangle_{\mathrm{MC}}}\Bigg\langle\Big[\bar{M}_{\bar{b}a}(\omega_{\bar{b}},\omega_a)\overline{\mathcal{G}M\mathcal{G}}_{c\bar{d}}-\overline{\mathcal{G}M}_{ca}(\omega_a)\overline{M\mathcal{G}}_{\bar{b}\bar{d}}(\omega_{\bar{b}})\Big]\operatorname{sign}(\mathcal{C})\Bigg\rangle_{\mathrm{MC}}.\]

Constructing \(\chi^{3,r}\) from \(M^{3,r}\)

With

\[M_{abcd}^{3,pp,conn}(\omega_a,\omega_c)=M_{abcd}^{3,pp}(\omega_a,\omega_c)-\overline{\mathcal{G}M}_{ba}(\omega_a)\overline{\mathcal{G}M}_{dc}(\omega_c)+\overline{\mathcal{G}M}_{da}(\omega_a)\overline{\mathcal{G}M}_{bc}(\omega_c),\]

\[M_{a\bar{b}\bar{c}d}^{3,ph,conn}(\omega_a,\omega_{\bar{b}})=M_{a\bar{b}\bar{c}d}^{3,ph}(\omega_a,\omega_{\bar{b}})-\beta\delta_{\omega_a,\omega_{\bar{b}}}\bar{M}_{\bar{b}a}(\omega_a)\overline{\mathcal{G}M\mathcal{G}}_{c\bar{d}}+\overline{\mathcal{G}M}_{ca}(\omega_a)\overline{M\mathcal{G}}_{\bar{b}\bar{d}}(\omega_{\bar{b}}),\]

we have that

\[\chi_{\bar{a}b\bar{c}d}^{3,pp,conn}(\omega_a,\omega_c)=\mathcal{G}_{\bar{a}i}(\omega_a)\mathcal{G}_{\bar{c}j}(\omega_c)M_{ibjd}^{3,pp,conn}(\omega_a,\omega_c),\]

\[\chi_{\bar{a}b\bar{c}d}^{3,ph,conn}(\omega_{\bar{a}},\omega_b)=\mathcal{G}_{i\bar{a}}(\omega_{\bar{a}})\mathcal{G}_{b\bar{j}}(\omega_b)M_{i\bar{j}\bar{c}d}^{3,ph,conn}(\omega_{\bar{a}},\omega_b).\]

Note that \(M_{abcd}^{3,pp}(\omega_a,\omega_c)\) and \(M_{a\bar{b}\bar{c}d}^{3,ph}(\omega_a,\omega_{\bar{b}})\) can be either measured directly in frequencies, or calculated via a Fourier transforms

\[M_{abcd}^{3,pp}(\omega_a,\omega_c)=\int_{\tau_a,\tau_c}e^{-i\tau_a\omega_a}e^{-i\tau_c\omega_c}M_{abcd}^{3,pp}(\tau_a,\tau_c)\]

\[M_{a\bar{b}\bar{c}d}^{3,ph}(\omega_a,\omega_{\bar{b}})=\int_{\tau_a,\tau_{\bar{b}}}e^{-i\tau_a\omega_a}e^{i\tau_{\bar{b}}\omega_{\bar{b}}}M_{a\bar{b}\bar{c}d}^{3,ph}(\tau_a,\tau_{\bar{b}})\]

Mixed fermion-boson notation

Instead of using the frequencies of the fermionic operators, it is often convenient to work in a mixed fermion-boson notation. It is defined as

\[\begin{split}\tilde{M}_{abcd}^{3,pp}(\omega,\Omega)=M_{abcd}^{3,pp}(\omega,\Omega-\omega) & =\int_{\tau_a,\tau_c}e^{i(\tau_c-\tau_a)\omega}e^{-i\tau_c\Omega}M_{abcd}^{3,pp}(\tau_a,\tau_c)\\ & =\int_{-\beta}^0d\tau'\int_0^\beta d\tau_a e^{-i(\tau'+\tau_a)\omega}e^{i\tau'\Omega}M_{abcd}^{3,pp}(\tau_a,-\tau')\\ & =-\int_0^\beta d\tau'\int_0^\beta d\tau_a e^{-i(\tau'+\tau_a)\omega}e^{i\tau'\Omega}M_{abcd}^{3,pp}(\tau_a,\beta-\tau')\\ & =\int_0^\beta d\tau'\int_{-\tau'}^{-\beta-\tau'}d\tau e^{i\tau\omega}e^{i\tau'\Omega}M_{abcd}^{3,pp}(-\tau-\tau',\beta-\tau')\\ & =\int_0^\beta d\tau'\int_{-\tau'}^{\beta-\tau'}d\tau e^{i\tau\omega}e^{i\tau'\Omega}M_{abcd}^{3,pp}(\beta-\tau-\tau',\beta-\tau')\\\end{split}\]

\[\tilde{M}_{a\bar{b}\bar{c}d}^{3,ph}(\omega,\Omega)=M_{a\bar{b}\bar{c}d}^{3,ph}(\omega,\Omega+\omega)=\int_{\tau_a,\tau_{\bar{b}}}e^{i(\tau_{\bar{b}}-\tau_a)\omega}e^{i\tau_{\bar{b}}\Omega}M_{a\bar{b}\bar{c}d}^{3,ph}(\tau_a,\tau_{\bar{b}}),\]

and accordingly for \(\chi^{3,r}\).

Equal-time correlation functions \(\chi^2\)

We now consider the case that two pairs of times are equal. In this case we measure in imaginary time, and a Fourier transform is performed in a post-processing step. We require again an intermediate object

\[\overline{\mathcal{G}M\mathcal{G}}_{\beta\bar{\alpha}}\equiv\overline{\mathcal{G}M\mathcal{G}}_{b\bar{a}}(\tau_b,\tau_a)\equiv\sum_{ij}\mathcal{G}_{b\bar{c}_j}(\tau_b-\tau_j)[\hat{\mathcal{G}}^{-1}]_{ji}\mathcal{G}_{c_i\bar{a}}(\tau_i-\tau_{\bar{a}}).\]

The measurements then read

\[M_{\bar{a}b\bar{c}d}^{2,pp}(\tau)=\frac{1}{\langle\operatorname{sign}(\mathcal{C})\rangle_{\mathrm{MC}}}\Bigg\langle\Big[\overline{\mathcal{G}M\mathcal{G}}_{b\bar{a}}(0^+,\tau^+)\overline{\mathcal{G}M\mathcal{G}}_{d\bar{c}}(0,\tau)-\overline{\mathcal{G}M\mathcal{G}}_{d\bar{a}}(0,\tau^+)\overline{\mathcal{G}M\mathcal{G}}_{b\bar{c}}(0^+,\tau)\Big]\operatorname{sign}(\mathcal{C})\Bigg\rangle_{\mathrm{MC}},\]

\[M_{\bar{a}b\bar{c}d}^{2,ph}(\tau)=\frac{1}{\langle\operatorname{sign}(\mathcal{C})\rangle_{\mathrm{MC}}}\Bigg\langle\Big[\overline{\mathcal{G}M\mathcal{G}}_{b\bar{a}}(\tau,\tau^+)\overline{\mathcal{G}M\mathcal{G}}_{d\bar{c}}(0,0^+)-\overline{\mathcal{G}M\mathcal{G}}_{d\bar{a}}(0,\tau^+)\overline{\mathcal{G}M\mathcal{G}}_{b\bar{c}}(\tau,0^+)\Big]\operatorname{sign}(\mathcal{C})\Bigg\rangle_{\mathrm{MC}}.\]

Constructing \(\chi^{2,r}\) from \(M^{2,r}\)

We have that

\[\chi_{\bar{a}b\bar{c}d}^{2,pp,conn}(\tau)=M_{\bar{a}b\bar{c}d}^{2,pp}(\tau)-\overline{\mathcal{G}M\mathcal{G}}_{b\bar{a}}(\beta-\tau)\overline{\mathcal{G}M\mathcal{G}}_{d\bar{c}}(\beta-\tau)+\overline{\mathcal{G}M\mathcal{G}}_{d\bar{a}}(\beta-\tau)\overline{\mathcal{G}M\mathcal{G}}_{b\bar{c}}(\beta-\tau),\]

\[\chi_{\bar{a}b\bar{c}d}^{2,ph,conn}(\tau)=M_{\bar{a}b\bar{c}d}^{2,ph}(\tau)-\overline{\mathcal{G}M\mathcal{G}}_{b\bar{a}}(\beta^-)\overline{\mathcal{G}M\mathcal{G}}_{d\bar{c}}(\beta^-)-\overline{\mathcal{G}M\mathcal{G}}_{d\bar{a}}(\beta-\tau)\overline{\mathcal{G}M\mathcal{G}}_{b\bar{c}}(\tau).\]

[1]

Careful: The Fourier conventions for \(M\)-objects are indeed reversed, meaning that barred indices are transformed to frequencies with \(\int_\tau e^{iw\tau}\).