The CTSEG double expansion algorithm

Double expansion

The segment-picture continuous-time quantum Monte Carlo algorithm (CTSEG), in its double expansion version, is based on an expansion of the partition function in powers of the hybridization function \(\Delta(\tau)\) and perpendicular spin-spin interaction \(J_{\perp} (\tau)\). Details can be found in [1] [2]. The partition function is \(Z = \int [Dc][D\overline{c}]e^{-\mathcal{S}[c,\overline{c}]}\), with the action \(\mathcal{S}\) given by

\[\begin{split}\begin{split} \mathcal{S} &= \iint_0^{\beta} \mathrm{d} \tau \mathrm{d} \tau' \sum_{a,b} \left\{ \overline{c}_{a\sigma} (\tau) \left( (\partial_{\tau} + \epsilon_{a\sigma} - \mu)\delta_{ab}^{\sigma \sigma'} \delta_{\tau - \tau'} + \Delta_{ab}^{\sigma \sigma'}(\tau - \tau')\right) c_{b\sigma'}(\tau') \right\} \\ &+ \frac{1}{2} \iint_0^{\beta} \mathrm{d} \tau \mathrm{d} \tau' \sum_{a,b} \mathcal{U}_{ab}(\tau - \tau') n_a(\tau) n_b(\tau') + \frac{1}{2} \iint_0^{\beta} \mathrm{d} \tau \mathrm{d} \tau' \sum_{a, \xi = x, y, z} s_a^{\xi}(\tau) \mathcal{J}_a^{\xi}(\tau - \tau') s_a^{\xi} (\tau') \end{split}\end{split}\]

Here \(\beta\) is the inverse temperature, \(a\) denote orbital indices, \(\sigma\) spin indices (\(\sigma = \uparrow, \downarrow\)), \(n_a \equiv \sum_{\sigma} n_{a\sigma}\), \(s_a^{\xi} \equiv \frac{1}{2} \sum_{\sigma \sigma'} \overline{c}_{a\sigma} \sigma_{\sigma \sigma'}^{\xi} c_{a \sigma'}\) and \(\sigma^{\xi}\) are the Pauli matrices. \(\overline{c}_{a\sigma}(\tau)\) and \(c_{a\sigma}(\tau)\) are the \(\beta\)-antiperiodic Grassmann fields corresponding to the fermion creation and annihilation operators on the impurity, respectively. The \(\epsilon_{a\sigma}\) are orbital energies and \(mu\) is the chemical potential.

This action can be recast as

\[\begin{split}\begin{split} \mathcal{S} &= \iint_0^{\beta} \mathrm{d} \tau \mathrm{d} \tau' \sum_{a,b} \left\{ \overline{c}_{a\sigma} (\tau) \left( (\partial_{\tau} + \epsilon_{a\sigma} - \mu)\delta_{ab}^{\sigma \sigma'} \delta_{\tau - \tau'} + \Delta_{ab}^{\sigma \sigma'}(\tau - \tau')\right) c_{b\sigma'}(\tau') \right\} \\ &+ \frac{1}{2} \iint_0^{\beta} \mathrm{d} \tau \mathrm{d} \tau' \sum_{u,v} \mathcal{U}_{uv}(\tau - \tau') n_u(\tau) n_v(\tau') + \frac{1}{2} \iint_0^{\beta} \mathrm{d} \tau \mathrm{d} \tau' \sum_{a} \mathcal{J}_a^{\perp}(\tau - \tau') s_a^{+}(\tau) s_a^{-} (\tau') \end{split}\end{split}\]

where

\[\mathcal{U}_{uv}(\tau - \tau') = \mathcal{U}_{ab}(\tau - \tau') + (-1)^{\sigma \sigma'} \frac{1}{4} \mathcal{J}_a^z(\tau) \delta_{ab},\]

\(s^{\pm} = S_x \pm i s_y\) and \(\mathcal{J}^{\perp} \equiv \mathcal{J}^x = \mathcal{J}^y\). The CTSEG solver stochastically explores the terms (or configurations) generated by the expansion of \(\mathcal{S}\) in powers of \(\Delta(\tau)\) and \(\mathcal{J}^{\perp}(\tau)\) and samples the observables of interest (e.g. the Green’s function) every few configurations.

Note

Our CTSEG implementation supports the \(\mathcal{J}^{\perp}(\tau)\) expansion only for a single orbital. For multiple orbitals, an expansion only in \(\Delta(\tau)\) is possible. For a single orbital, it is also possible to carry out an expansion in \(\mathcal{J}^{\perp}(\tau)\) only (i.e., with \(\Delta(\tau) = 0\)).

Configuration

A term of the double expansion takes the general form

\[w = \mathrm{det} [\Delta] \prod_{i = 1}^{k} \left[ - \frac{\mathcal{J}_{\perp}(q(\alpha_i) - \beta_i) }{2} \right] w_{\rm loc},\]

where \(w_{\rm loc}\) is the local trace factor:

\[w_{\rm loc} = \langle \psi \vert e^{-S_{\rm loc}} c(v_1)\overline{c}(u_1) \dots c(v_n) \overline{c}(u_n) s^+(\alpha_1)s^-(\beta_1) \dots s^+(\alpha_k) s^-(\beta_k) \vert \psi \rangle.\]

We have denoted \(\vert \Psi \rangle\) the basis states of the isolated impurity and \(S_{\rm loc}\) the corresponding action. \([\Delta]\) is the matrix with elements \([\Delta]_{ij} = \Delta(u_i - v_j)\) and \(q\) is a permutation of \((1, \dots, k)\). We chose to regroup all the possible permutations of the “hybridized” operators into a determinant because it allows for numerically fast updates, while there is no advantage to regrouping the permutations of the spin operators into a permanent.

We will designate by “color” the set of indices carried by an operator in addition to its time (typically, spin and orbital). Under the applicability conditions of CTSEG (the interaction commutes with the density operator) the configuration can be represented as a set of imaginary time segments for each color: within a segment, the color is occupied. Then, the local trace factor supports a simple expression in terms of the overlaps of those segments between different colors (up to a sign, discussed below):

\[w_{\rm loc} = w_{\rm st}(\tilde U_{ab}, \tilde \mu_a) \cdot \exp \left[ \sum_{ij} s_i s_j K_{a(i)b(j)} (\tau_i - \tau_j) \right],\]

where \(s_i = 1\) for creation operators and -1 otherwise. We split the interaction into a static and a dynamic part: \(\mathcal{U}_{ab}(\tau) = U_{ab} + D_{ab}(\tau)\). The dynamical interaction kernel \(K_{ab}\) is defined by \(K_{ab}''(\tau) = D_{ab}(\tau)\) and \(K_{ab}(\tau = 0^+) = K_{ab}(\tau = \beta^-) = 0\). The static trace factor is

\[w_{\rm st}(U_{ab}, \mu_a) = \exp \left[ -\sum_{a \neq b} U_{ab} O_{ab} + \sum_a \mu_a \ell_a \right].\]

Here \(O_{ab}\) is the overlap between the segments of colors \(a\) and \(b\), and \(\ell_a\) is the occupied length in color a. In the presence of dynamical interactions, the static trace factor is computed with a renormalized interaction and chemical potential: \(\tilde U_{ab} = U_{ab} - 2 K'_{ab}(0)\) and \(\tilde \mu_{a} = \mu - \epsilon_a + K'_{aa}(0)\).

Sign

The sign of the Monte Carlo weight depends on the initial ordering of the operators in the expansion. The total sign is the product of the sign of the determinant and the sign of the trace, and the ordering of the spin operators does not influence the sign. The trace is positive if all the hybridized operators (whatever their nature) are ordered by decreasing (from left to right) time and color. In practice, we define the matrix \([\Delta]\) in such a way that the times corresponding to its rows (creation operators) and columns (annihilation operators) are both in increasing order, regardless of color: this fixes the sign of the determinant. We then compute the sign of the permutation that takes the sequence \([c c^{\dagger} \dots c c^{\dagger} ]\) with the \(c\) and \(c^{\dagger}\) increasing-time-ordered to the decreasing-time-and-color-ordered sequence of operators.