API

greens!(
    G::AbstractArray{C,D},
    a::Int, b::Int, unit_cell::UnitCell{D}, lattice::Lattice{D},
    G_τ0::AbstractMatrix{T},
    sgn::T=one(T)
) where {D, C<:Number, T<:Number}

Measure the unequal time Green's function averaged over translation symmetry

\[G_{\sigma,\mathbf{r}}^{a,b}(\tau)=\frac{1}{N}\sum_{\mathbf{i}}G_{\sigma,\mathbf{i}+\mathbf{r},\mathbf{i}}^{a,b}(\tau,0) =\frac{1}{N}\sum_{\mathbf{i}}\langle\hat{\mathcal{T}}\hat{a}_{\sigma,\mathbf{i}+\mathbf{r}}^{\phantom{\dagger}}(\tau)\hat{b}_{\sigma,\mathbf{i}}^{\dagger}(0)\rangle,\]

with the result being added to G.

Fields

• G::AbstractArray{C,D}: Array the green's function $G_{\sigma,\mathbf{r}}^{a,b}(\tau)$ is written to.
• a::Int: Index specifying an orbital species in the unit cell.
• b::Int: Index specifying an orbital species in the unit cell.
• unit_cell::UnitCell{D}: Defines unit cell.
• lattice::Lattice{D}: Specifies size of finite lattice.
• G_τ0::AbstractMatrix{T}: The matrix $G(\tau,0).$
• sgn::T=one(T): The sign of the weight appearing in a DQMC simulation.

density_correlation!(
    DD::AbstractArray{C,D},
    a::Int, b::Int, unit_cell::UnitCell{D}, lattice::Lattice{D},
    Gup_τ0::AbstractMatrix{T}, Gup_0τ::AbstractMatrix{T},
    Gup_ττ::AbstractMatrix{T}, Gup_00::AbstractMatrix{T},
    Gdn_τ0::AbstractMatrix{T}, Gdn_0τ::AbstractMatrix{T},
    Gdn_ττ::AbstractMatrix{T}, Gdn_00::AbstractMatrix{T},
    sgn::T=one(T)
) where {D, C<:Number, T<:Number}

Calculate the unequal-time density-density (charge) correlation function

\[\begin{align*} \mathcal{D}_{\mathbf{r}}^{a,b}(\tau) & = \frac{1}{N}\sum_{\mathbf{i}} \mathcal{D}_{\mathbf{i}+\mathbf{r},\mathbf{i}}^{a,b}(\tau,0)\\ & = \frac{1}{N}\sum_{\mathbf{i}} \langle \hat{n}_{a,\mathbf{i} + \mathbf{r}}(\tau)\hat{n}_{b,\mathbf{i}}(0) \rangle, \end{align*}\]

where $\hat{n}_{b,\mathbf{i}} = (\hat{n}_{\uparrow, b, \mathbf{i}} + \hat{n}_{\downarrow, b, \mathbf{i}})$ and $\hat{n}_{\sigma, b,\mathbf{i}} = \hat{b}^\dagger_{\sigma, \mathbf{i}} \hat{b}_{\sigma, \mathbf{i}}$ is the number operator for an electron with spin $\sigma$ on orbital $b$ in unit cell $\mathbf{i}$, with the result being added to the array DD.

Fields

• DD::AbstractArray{C,D}: Array the density correlation function $\mathcal{D}_{\mathbf{r}}^{a,b}(\tau)$ is added to.
• a::Int: Index specifying an orbital species in the unit cell.
• b::Int: Index specifying an orbital species in the unit cell.
• unit_cell::UnitCell{D}: Defines unit cell.
• lattice::Lattice{D}: Specifies size of finite lattice.
• Gup_τ0::AbstractMatrix{T}: The matrix $G_{\uparrow}(\tau,0).$
• Gup_0τ::AbstractMatrix{T}: The matrix $G_{\uparrow}(0,\tau).$
• Gup_ττ::AbstractMatrix{T}: The matrix $G_{\uparrow}(\tau,\tau).$
• Gup_00::AbstractMatrix{T}: The matrix $G_{\uparrow}(0,0).$
• Gdn_τ0::AbstractMatrix{T}: The matrix $G_{\downarrow}(\tau,0).$
• Gdn_0τ::AbstractMatrix{T}: The matrix $G_{\downarrow}(0,\tau).$
• Gdn_ττ::AbstractMatrix{T}: The matrix $G_{\downarrow}(\tau,\tau).$
• Gdn_00::AbstractMatrix{T}: The matrix $G_{\downarrow}(0,0).$
• sgn::T=one(T): The sign of the weight appearing in a DQMC simulation.

density_correlation!(
    DD::AbstractArray{C,D},
    a::Int, b::Int, unit_cell::UnitCell{D}, lattice::Lattice{D},
    Gσ_τ0::AbstractMatrix{T}, Gσ_0τ::AbstractMatrix{T},
    Gσ_ττ::AbstractMatrix{T}, Gσ′_00::AbstractMatrix{T},
    σ::Int, σ′::Int, sgn::T=one(T)
) where {D, C<:Number, T<:Number}

Calculate the spin-resolved unequal-time density correlation function

\[\begin{align*} \mathcal{D}_{\mathbf{r}}^{(a,\sigma),(b,\sigma')}(\tau) = \frac{1}{N}\sum_{\mathbf{i}} & \mathcal{D}_{\mathbf{i}+\mathbf{r},\mathbf{i}}^{(a,\sigma),(b,\sigma')}(\tau,0) \\ = \frac{1}{N}\sum_{\mathbf{i}} & \langle\hat{n}_{\sigma,a,\mathbf{i}+\mathbf{r}}(\tau)\hat{n}_{\sigma',b,\mathbf{r}}(0)\rangle, \end{align*}\]

where $\hat{n}_{\sigma', b,\mathbf{i}} = \hat{b}^\dagger_{\sigma', \mathbf{i}} \hat{b}_{\sigma, \mathbf{i}}$ is the number operator for an electron with spin $\sigma'$ on orbital $b$ in unit cell $\mathbf{i}$, with the result being added to the array DD.

Fields

• DD::AbstractArray{C,D}: The array the spin-resolved density correlation $\mathcal{D}_{\mathbf{r}}^{(a,\sigma),(b,\sigma')}(\tau)$ is added to.
• a::Int: Index specifying an orbital species in the unit cell.
• b::Int: Index specifying an orbital species in the unit cell.
• unit_cell::UnitCell{D}: Defines unit cell.
• lattice::Lattice{D}: Specifies size of finite lattice.
• Gσ_τ0::AbstractMatrix{T}: The matrix $G_{\sigma}(\tau,0).$
• Gσ_0τ::AbstractMatrix{T}: The matrix $G_{\sigma}(0,\tau).$
• Gσ_ττ::AbstractMatrix{T}: The matrix $G_{\sigma}(\tau,\tau).$
• Gσ′_00::AbstractMatrix{T}: The matrix $G_{\sigma'}(0,0).$
• σ::Int: The electron spin appearing in the $\langle \hat{n} \rangle_{\sigma,a,\mathbf{i}+\mathbf{r}}$ density operator.
• σ′::Int: The electron spin appearing in the $\langle \hat{n} \rangle_{\sigma',b,\mathbf{i}}$ density operator.
• sgn::T=one(T): The sign of the weight appearing in a DQMC simulation.

spin_x_correlation!(
    SxSx::AbstractArray{C,D}, a::Int, b::Int, unit_cell::UnitCell{D}, lattice::Lattice{D},
    Gup_τ0::AbstractMatrix{T}, Gup_0τ::AbstractMatrix{T},
    Gdn_τ0::AbstractMatrix{T}, Gdn_0τ::AbstractMatrix{T},
    sgn::T=one(T)
) where {D, C<:Number, T<:Number}

Calculate the unequal-time spin-spin correlation function in the $\hat{x}$ direction, given by

\[\mathcal{S}_{x,\mathbf{r}}^{a,b}(\tau)=\frac{1}{N}\sum_{\mathbf{i}}\mathcal{S}_{x,\mathbf{i}+\mathbf{r},\mathbf{i}}^{ab}(\tau,0) =\frac{1}{N}\sum_{\mathbf{i}}\big\langle\hat{S}_{x,a,\mathbf{i}+\mathbf{r}}(\tau)\hat{S}_{x,b,\mathbf{i}}(0)\big\rangle,\]

where the spin-$\hat{x}$ operator is given by

\[\begin{align*} \hat{S}_{x,\mathbf{i},a}= & (\hat{a}_{\uparrow,\mathbf{i}}^{\dagger},\hat{a}_{\downarrow,\mathbf{i}}^{\dagger})\left[\begin{array}{cc} 0 & 1\\ 1 & 0 \end{array}\right]\left(\begin{array}{c} \hat{a}_{\uparrow,\mathbf{i}}\\ \hat{a}_{\downarrow,\mathbf{i}} \end{array}\right)\\ = & \hat{a}_{\uparrow,\mathbf{i}}^{\dagger}\hat{a}_{\downarrow,\mathbf{i}}+\hat{a}_{\downarrow,\mathbf{i}}^{\dagger}\hat{a}_{\uparrow,\mathbf{i}}. \end{align*}\]

Fields

• SxSx::AbstractArray{C,D}: Array the spin-$x$ correlation function $\mathcal{S}_{x,\mathbf{r}}^{a,b}(\tau)$ is added to.
• a::Int: Index specifying an orbital species in the unit cell.
• b::Int: Index specifying an orbital species in the unit cell.
• unit_cell::UnitCell{D}: Defines unit cell.
• lattice::Lattice{D}: Specifies size of finite lattice.
• Gup_τ0::AbstractMatrix{T}: The matrix $G_{\uparrow}(\tau,0).$
• Gup_0τ::AbstractMatrix{T}: The matrix $G_{\uparrow}(0,\tau).$
• Gdn_τ0::AbstractMatrix{T}: The matrix $G_{\downarrow}(\tau,0).$
• Gdn_0τ::AbstractMatrix{T}: The matrix $G_{\downarrow}(0,\tau).$
• sgn::T=one(T): The sign of the weight appearing in a DQMC simulation.

```

spin_z_correlation!(
    SzSz::AbstractArray{C,D},
    a::Int, b::Int, unit_cell::UnitCell{D}, lattice::Lattice{D},
    Gup_τ0::AbstractMatrix{T}, Gup_0τ::AbstractMatrix{T},
    Gup_ττ::AbstractMatrix{T}, Gup_00::AbstractMatrix{T},
    Gdn_τ0::AbstractMatrix{T}, Gdn_0τ::AbstractMatrix{T},
    Gdn_ττ::AbstractMatrix{T}, Gdn_00::AbstractMatrix{T},
    sgn::T=one(T)
) where {D, C<:Complex, T<:Number}

Calculate the unequal-time spin-spin correlation function in the $\hat{z}$ direction, given by

\[\begin{align*} \mathcal{S}_{z,\mathbf{r}}^{a,b}(\tau)=\frac{1}{N}\sum_{\mathbf{i}}\mathcal{S}_{z,\mathbf{i}+\mathbf{r},\mathbf{i}}^{ab}(\tau,0) = & \frac{1}{N}\sum_{\mathbf{i}}\big\langle\hat{S}_{z,a,\mathbf{i}+\mathbf{r}}(\tau)\hat{S}_{z,b,\mathbf{i}}(0)\big\rangle, \end{align*}\]

where the spin-$\hat{z}$ operator is given by

\[\begin{align*} \hat{S}_{z,a,\mathbf{i}}= & (\hat{a}_{\uparrow,\mathbf{i}}^{\dagger},\hat{a}_{\downarrow,\mathbf{i}}^{\dagger})\left[\begin{array}{cc} 1 & 0\\ 0 & -1 \end{array}\right]\left(\begin{array}{c} \hat{a}_{\uparrow,\mathbf{i}}\\ \hat{a}_{\downarrow,\mathbf{i}} \end{array}\right)\\ = & \hat{n}_{\uparrow,a,\mathbf{i}}-\hat{n}_{\downarrow,a,\mathbf{i}}. \end{align*}\]

Fields

• SzSz::AbstractArray{C,D}: Array the spin-$z$ correlation function $\mathcal{S}_{z,\mathbf{r}}^{a,b}(\tau)$ is added to.
• a::Int: Index specifying an orbital species in the unit cell.
• b::Int: Index specifying an orbital species in the unit cell.
• unit_cell::UnitCell{D}: Defines unit cell.
• lattice::Lattice{D}: Specifies size of finite lattice.
• Gup_τ0::AbstractMatrix{T}: The matrix $G_{\uparrow}(\tau,0).$
• Gup_0τ::AbstractMatrix{T}: The matrix $G_{\uparrow}(0,\tau).$
• Gup_ττ::AbstractMatrix{T}: The matrix $G_{\uparrow}(\tau,\tau).$
• Gup_00::AbstractMatrix{T}: The matrix $G_{\uparrow}(0,0).$
• Gdn_τ0::AbstractMatrix{T}: The matrix $G_{\downarrow}(\tau,0).$
• Gdn_0τ::AbstractMatrix{T}: The matrix $G_{\downarrow}(0,\tau).$
• Gdn_ττ::AbstractMatrix{T}: The matrix $G_{\downarrow}(\tau,\tau).$
• Gdn_00::AbstractMatrix{T}: The matrix $G_{\downarrow}(0,0).$
• sgn::T=one(T): The sign of the weight appearing in a DQMC simulation.

pair_correlation!(
    P::AbstractArray{C,D},
    b″::Bond{D}, b′::Bond{D}, unit_cell::UnitCell{D}, lattice::Lattice{D},
    Gup_τ0::AbstractMatrix{T}, Gdn_τ0::AbstractMatrix{T},
    sgn::T=one(T)
) where {D, C<:Number, T<:Number}

Calculate the unequal-time pair correlation function

\[\mathcal{P}_{\mathbf{r}}^{(a,b,r''),(c,d,r')}(\tau)=\frac{1}{N}\sum_{\mathbf{i}}\mathcal{P}_{\mathbf{i}+\mathbf{r},\mathbf{i}}^{(a,b,r''),(c,d,r')}(\tau,0) =\frac{1}{N}\sum_{\mathbf{i}}\langle\hat{\Delta}_{\mathbf{i}+\mathbf{r},a,b,\mathbf{r}''}(\tau)\hat{\Delta}_{\mathbf{i},c,d,\mathbf{r}'}^{\dagger}(0)\rangle,\]

where the bond b″ defines the pair creation operator

\[\hat{\Delta}_{\mathbf{i},a,b,\mathbf{r}''}^{\dagger}=\hat{a}_{\uparrow,\mathbf{i}+\mathbf{r}''}^{\dagger}\hat{b}_{\downarrow,\mathbf{i}}^{\dagger},\]

and the bond b′ defines the pair creation operator

\[\hat{\Delta}_{\mathbf{i},c,d,\mathbf{r}'}^{\dagger}=\hat{c}_{\uparrow,\mathbf{i}+\mathbf{r}'}^{\dagger}\hat{d}_{\downarrow,\mathbf{i}}^{\dagger}.\]

Fields

• P::AbstractArray{C,D}: Array the pair correlation function $\mathcal{P}_{\mathbf{r}}^{(a,b,r''),(c,d,r')}(\tau)$ is added to.
• b″::Bond{D}: Bond defining pair annihilation operator appearing in pair correlation function.
• b′::Bond{D}: Bond defining pair creation operator appearing in pair correlation function.
• unit_cell::UnitCell{D}: Defines unit cell.
• lattice::Lattice{D}: Specifies size of finite lattice.
• Gup_τ0::AbstractMatrix{T}: The matrix $G_{\uparrow}(\tau,0).$
• Gdn_τ0::AbstractMatrix{T}: The matrix $G_{\downarrow}(\tau,0).$
• sgn::T=one(T): The sign of the weight appearing in a DQMC simulation.

bond_correlation!(
    BB::AbstractArray{C,D},
    b′::Bond{D}, b″::Bond{D}, unit_cell::UnitCell{D}, lattice::Lattice{D},
    Gup_τ0::AbstractMatrix{T}, Gup_0τ::AbstractMatrix{T},
    Gup_ττ::AbstractMatrix{T}, Gup_00::AbstractMatrix{T},
    Gdn_τ0::AbstractMatrix{T}, Gdn_0τ::AbstractMatrix{T},
    Gdn_ττ::AbstractMatrix{T}, Gdn_00::AbstractMatrix{T},
    sgn::T=one(T)
) where {D, C<:Number, T<:Number}

Calculate the uneqaul-time bond-bond correlation function

\[\begin{align*} \mathcal{B}_{\mathbf{r}}^{(\mathbf{r}',a,b),(\mathbf{r}'',c,d)}(\tau) = & \frac{1}{N}\sum_{\mathbf{i}} \langle[\hat{B}_{\uparrow,\mathbf{i}+\mathbf{r},(\mathbf{r}',a,b)}(\tau)+\hat{B}_{\downarrow,\mathbf{i}+\mathbf{r},(\mathbf{r}',a,b)}(\tau)] \cdot[\hat{B}_{\uparrow,\mathbf{i},(\mathbf{r}'',c,d)}(0)+\hat{B}_{\downarrow,\mathbf{i},(\mathbf{r}'',c,d)}(0)]\rangle \end{align*}\]

where the

\[\hat{B}_{\sigma,\mathbf{i},(\mathbf{r},a,b)} = \hat{a}_{\sigma,\mathbf{i}+\mathbf{r}}^{\dagger}\hat{b}_{\sigma,\mathbf{i}}^{\phantom{\dagger}} + \hat{b}_{\sigma,\mathbf{i}}^{\dagger}\hat{a}_{\sigma,\mathbf{i}+\mathbf{r}}^{\phantom{\dagger}}\]

is the bond operator.

Fields

• BB::AbstractArray{C,D}: Array the bond correlation function $\mathcal{B}_{\mathbf{r}}^{(\mathbf{r}',a,b),(\mathbf{r}'',c,d)}(\tau)$ is added to.
• b′::Bond{D}: Bond defining the bond operator appearing on the left side of the bond correlation function.
• b″::Bond{D}: Bond defining the bond operator appearing on the right side of the bond correlation function.
• unit_cell::UnitCell{D}: Defines unit cell.
• lattice::Lattice{D}: Specifies size of finite lattice.
• Gup_τ0::AbstractMatrix{T}: The matrix $G_{\uparrow}(\tau,0).$
• Gup_0τ::AbstractMatrix{T}: The matrix $G_{\uparrow}(0,\tau).$
• Gup_ττ::AbstractMatrix{T}: The matrix $G_{\uparrow}(\tau,\tau).$
• Gup_00::AbstractMatrix{T}: The matrix $G_{\uparrow}(0,0).$
• Gdn_τ0::AbstractMatrix{T}: The matrix $G_{\downarrow}(\tau,0).$
• Gdn_0τ::AbstractMatrix{T}: The matrix $G_{\downarrow}(0,\tau).$
• Gdn_ττ::AbstractMatrix{T}: The matrix $G_{\downarrow}(\tau,\tau).$
• Gdn_00::AbstractMatrix{T}: The matrix $G_{\downarrow}(0,0).$
• sgn::T=one(T): The sign of the weight appearing in a DQMC simulation.

bond_correlation!(
    BB::AbstractArray{C,D},
    b′::Bond{D}, b″::Bond{D}, unit_cell::UnitCell{D}, lattice::Lattice{D},
    Gσ′_τ0::AbstractMatrix{T}, Gσ′_0τ::AbstractMatrix{T},
    Gσ′_ττ::AbstractMatrix{T}, Gσ″_00::AbstractMatrix{T},
    σ′::Int, σ″::Int, sgn::T=one(T)
) where {D, C<:Number, T<:Number}

Calculate the spin-resolved uneqaul-time bond-bond correlation function

\[\begin{align*} \mathcal{B}_{\mathbf{r}}^{(\mathbf{r}',a,b,\sigma'),(\mathbf{r}'',c,d,\sigma'')}(\tau) = \frac{1}{N}\sum_{\mathbf{i}} & \mathcal{B}_{\mathbf{i}+\mathbf{r},\mathbf{i}}^{(\mathbf{r}',a,b,\sigma'),(\mathbf{r}'',c,d,\sigma'')}(\tau,0)\\ =\frac{1}{N}\sum_{\mathbf{i}} & \langle\hat{B}_{\sigma',\mathbf{i}+\mathbf{r},(\mathbf{r}',a,b)}(\tau)\hat{B}_{\sigma'',\mathbf{i},(\mathbf{r}'',c,d)}(0)\rangle, \end{align*}\]

where

\[\begin{align*} \hat{B}_{\sigma,\mathbf{i},(\mathbf{r},a,b)} & = \hat{a}_{\sigma,\mathbf{i}+\mathbf{r}}^{\dagger}\hat{b}_{\sigma,\mathbf{i}}^{\phantom{\dagger}}+\hat{b}_{\sigma,\mathbf{i}}^{\dagger}\hat{a}_{\sigma,\mathbf{i}+\mathbf{r}}^{\phantom{\dagger}}\\ & = -\hat{a}_{\sigma,\mathbf{i}+\mathbf{r}}^{\phantom{\dagger}}\hat{b}_{\sigma,\mathbf{i}}^{\dagger}-\hat{b}_{\sigma,\mathbf{i}}^{\phantom{\dagger}}\hat{a}_{\sigma,\mathbf{i}+\mathbf{r}}^{\dagger} \end{align*}\]

is the bond operator.

Fields

• BB::AbstractArray{C,D}: Array the spin-resolved bond correlation function $\mathcal{B}_{\mathbf{r}}^{(\mathbf{r}',a,b,\sigma'),(\mathbf{r}'',c,d,\sigma'')}(\tau)$ is added to.
• b′::Bond{D}: Bond defining the bond operator appearing on the left side of the bond correlation function.
• b″::Bond{D}: Bond defining the bond operator appearing on the right side of the bond correlation function.
• unit_cell::UnitCell{D}: Defines unit cell.
• lattice::Lattice{D}: Specifies size of finite lattice.
• Gσ′_τ0::AbstractMatrix{T}: The matrix $G_{\sigma'}(\tau,0).$
• Gσ′_0τ::AbstractMatrix{T}: The matrix $G_{\sigma'}(0,\tau).$
• Gσ′_ττ::AbstractMatrix{T}: The matrix $G_{\sigma'}(\tau,\tau).$
• Gσ″_00::AbstractMatrix{T}: The matrix $G_{\sigma''}(0,0).$
• σ′::Int: The electron spin appearing in the $\hat{B}_{\sigma',\mathbf{i}+\mathbf{r},(\mathbf{r}',a,b)}$ bond operator.
• σ″::Int: The electron spin appearing in the $\hat{B}_{\sigma'',\mathbf{i},(\mathbf{r}'',c,d)}$ bond operator.
• sgn::T=one(T): The sign of the weight appearing in a DQMC simulation.

current_correlation!(
    CC::AbstractArray{C,D},
    b′::Bond{D}, b″::Bond{D},
    tup′::AbstractArray{T,D}, tup″::AbstractArray{T,D},
    tdn′::AbstractArray{T,D}, tdn″::AbstractArray{T,D},
    unit_cell::UnitCell{D}, lattice::Lattice{D},
    Gup_τ0::AbstractMatrix{T}, Gup_0τ::AbstractMatrix{T},
    Gup_ττ::AbstractMatrix{T}, Gup_00::AbstractMatrix{T},
    Gdn_τ0::AbstractMatrix{T}, Gdn_0τ::AbstractMatrix{T},
    Gdn_ττ::AbstractMatrix{T}, Gdn_00::AbstractMatrix{T},
    sgn::T=one(T)
) where {D, C<:Number, T<:Number}

Calculate the uneqaul-time current-current correlation function

\[\mathcal{J}_{\mathbf{r}}^{(\mathbf{r}',a,b),(\mathbf{r}'',c,d)}(\tau) = \frac{1}{N}\sum_{\mathbf{i}} \langle[\hat{J}_{\uparrow,\mathbf{i}+\mathbf{r},(\mathbf{r}',a,b)}(\tau)+\hat{J}_{\downarrow,\mathbf{i}+\mathbf{r},(\mathbf{r}',a,b)}(\tau)] \cdot[\hat{J}_{\uparrow,\mathbf{i},(\mathbf{r}'',c,d)}(0)+\hat{J}_{\downarrow,\mathbf{i},(\mathbf{r}'',c,d)}(0)]\rangle,\]

where the current operator is given by

\[\begin{align*} \hat{J}_{\sigma,\mathbf{i},(\mathbf{r},a,b)}= & -{\rm i}t_{\mathbf{i}+\mathbf{r},\mathbf{i}}^{a,b}(\hat{a}_{\sigma,\mathbf{i}+\mathbf{r}}^{\dagger}\hat{b}_{\sigma,\mathbf{i}}-\hat{b}_{\sigma,\mathbf{i}}^{\dagger}\hat{a}_{\sigma,\mathbf{i}+\mathbf{r}})\\ = & -{\rm i}t_{\mathbf{i}+\mathbf{r},\mathbf{i}}^{a,b}(\hat{a}_{\sigma,\mathbf{i}+\mathbf{r}}\hat{b}_{\sigma,\mathbf{i}}^{\dagger}-\hat{b}_{\sigma,\mathbf{i}}\hat{a}_{\sigma,\mathbf{i}+\mathbf{r}}^{\dagger}). \end{align*}\]

Fields

• CC::AbstractArray{C,D}: Array the current correlation function $\mathcal{J}_{\mathbf{r}}^{(\mathbf{r}',a,b),(\mathbf{r}'',c,d)}(\tau)$ is added to.
• b′::Bond{D}: Bond defining the current operator appearing on the left side of the current correlation function.
• b″::Bond{D}: Bond defining the current operator appearing on the right side of the current correlation function.
• tup′::AbstractArray{T,D}: Spin up position and imaginary time dependent hopping amplitudes corresponding to bond b′.
• tup″::AbstractArray{T,D}: Spin up position and imaginary time dependent hopping amplitudes corresponding to bond b″.
• tdn′::AbstractArray{T,D}: Spin down position and imaginary time dependent hopping amplitudes corresponding to bond b′.
• tdn″::AbstractArray{T,D}: Spin down position and imaginary time dependent hopping amplitudes corresponding to bond b″.
• unit_cell::UnitCell{D}: Defines unit cell.
• lattice::Lattice{D}: Specifies size of finite lattice.
• Gup_τ0::AbstractMatrix{T}: The matrix $G_{\uparrow}(\tau,0).$
• Gup_0τ::AbstractMatrix{T}: The matrix $G_{\uparrow}(0,\tau).$
• Gup_ττ::AbstractMatrix{T}: The matrix $G_{\uparrow}(\tau,\tau).$
• Gup_00::AbstractMatrix{T}: The matrix $G_{\uparrow}(0,0).$
• Gdn_τ0::AbstractMatrix{T}: The matrix $G_{\downarrow}(\tau,0).$
• Gdn_0τ::AbstractMatrix{T}: The matrix $G_{\downarrow}(0,\tau).$
• Gdn_ττ::AbstractMatrix{T}: The matrix $G_{\downarrow}(\tau,\tau).$
• Gdn_00::AbstractMatrix{T}: The matrix $G_{\downarrow}(0,0).$
• sgn::T=one(T): The sign of the weight appearing in a DQMC simulation.

current_correlation!(
    CC::AbstractArray{C,D},
    b′::Bond{D}, b″::Bond{D},
    t′::AbstractArray{T,D}, t″::AbstractArray{T,D},
    unit_cell::UnitCell{D}, lattice::Lattice{D},
    Gσ′_τ0::AbstractMatrix{T}, Gσ′_0τ::AbstractMatrix{T},
    Gσ′_ττ::AbstractMatrix{T}, Gσ″_00::AbstractMatrix{T},
    σ′::Int, σ″::Int, sgn::T=one(T)
) where {D, C<:Number, T<:Number}

Calculate the spin-resolved uneqaul-time current-current correlation function

\[\begin{align*} \mathcal{J}_{\mathbf{r}}^{(\mathbf{r}',a,b,\sigma'),(\mathbf{r}'',c,d,\sigma'')}(\tau) & = \frac{1}{N}\sum_{\mathbf{i}}\mathcal{J}_{\mathbf{i}+\mathbf{r},\mathbf{i}}^{(\mathbf{r}',a,b,\sigma'),(\mathbf{r}'',c,d,\sigma'')}(\tau,0)\\ & = \frac{1}{N}\sum_{\mathbf{i}}\langle\hat{J}_{\sigma',\mathbf{i}+\mathbf{r},(\mathbf{r}',a,b)}(\tau)\hat{J}_{\sigma'',\mathbf{i},(\mathbf{r}'',c,d)}(0)\rangle, \end{align*}\]

where the spin-resolved current operator is given by

\[\begin{align*} \hat{J}_{\sigma,\mathbf{i},(\mathbf{r},a,b)} & = -{\rm i}t_{\sigma,\mathbf{i}+\mathbf{r},\mathbf{i}}^{a,b}(\hat{a}_{\sigma,\mathbf{i}+\mathbf{r}}^{\dagger}\hat{b}_{\sigma,\mathbf{i}}-\hat{b}_{\sigma,\mathbf{i}}^{\dagger}\hat{a}_{\sigma,\mathbf{i}+\mathbf{r}})\\ & = -{\rm i}t_{\sigma,\mathbf{i}+\mathbf{r},\mathbf{i}}^{a,b}(\hat{a}_{\sigma,\mathbf{i}+\mathbf{r}}\hat{b}_{\sigma,\mathbf{i}}^{\dagger}-\hat{b}_{\sigma,\mathbf{i}}\hat{a}_{\sigma,\mathbf{i}+\mathbf{r}}^{\dagger}). \end{align*}\]

Fields

• CC::AbstractArray{C,D}: Array the spin-resolved current correlation function $\mathcal{J}_{\mathbf{r}}^{(\mathbf{r}',a,b,\sigma'),(\mathbf{r}'',c,d,\sigma'')}(\tau)$ is added to.
• b′::Bond{D}: Bond defining the current operator appearing on the left side of the current correlation function.
• b″::Bond{D}: Bond defining the current operator appearing on the right side of the current correlation function.
• t′::AbstractArray{T,D}: Position and imaginary time dependent hopping amplitude associated with bond b′.
• t″::AbstractArray{T,D}: Position and imaginary time dependent hopping amplitude associated with bond b″.
• unit_cell::UnitCell{D}: Defines unit cell.
• lattice::Lattice{D}: Specifies size of finite lattice.
• Gσ′_τ0::AbstractMatrix{T}: The matrix $G_{\sigma'}(\tau,0).$
• Gσ′_0τ::AbstractMatrix{T}: The matrix $G_{\sigma'}(0,\tau).$
• Gσ′_ττ::AbstractMatrix{T}: The matrix $G_{\sigma'}(\tau,\tau).$
• Gσ″_00::AbstractMatrix{T}: The matrix $G_{\sigma''}(0,0).$
• σ′::Int: The electron spin appearing in the left current operator.
• σ″::Int: The electron spin appearing in the right current operator.
• sgn::T=one(T): The sign of the weight appearing in a DQMC simulation.

measure_n(G::AbstractMatrix{T}) where {T}

Measure the average density $\langle \hat{n}_\sigma \rangle$ given the equal-time Green's function matrix $G_\sigma(\tau,\tau).$

measure_n(G::AbstractMatrix{T}, a::Int, unit_cell::UnitCell) where {T}

Measure the average density $\langle \hat{n}_{\sigma,a} \rangle$ for orbital species $a,$ given the equal-time Green's function matrix $G_\sigma(\tau,\tau).$

measure_double_occ(Gup::AbstractMatrix{T}, Gdn::AbstractMatrix{T}) where {T}

Measure the double-occupancy $\langle \hat{n}_\uparrow \hat{n}_\downarrow \rangle$ given both the spin-up and spin-down equal-time Green's function matrices $G_\uparrow(\tau,\tau)$ and $G_\downarrow(\tau,\tau)$ respectively.

measure_double_occ(Gup::AbstractMatrix{T}, Gdn::AbstractMatrix{T}, a::Int, unit_cell::UnitCell) where {T}

Measure the double-occupancy $\langle \hat{n}_{\uparrow,a} \hat{n}_{\downarrow,a} \rangle$ for orbital species $a,$ given both the spin-up and spin-down equal-time Green's function matrices $G_\uparrow(\tau,\tau)$ and $G_\downarrow(\tau,\tau)$ respectively.

measure_N(G::AbstractMatrix{T}) where {T}

Measure the total particle number $\langle \hat{N}_\sigma \rangle$ given an equal-time Green's function matrix $G_\sigma(\tau,\tau).$

measure_N(G::AbstractMatrix{T}, a::Int, unit_cell::UnitCell) where {T}

Measure the total particle number $\langle \hat{N}_{\sigma,a} \rangle$ in orbital species $a,$ given an equal-time Green's function matrix $G_\sigma(\tau,\tau).$

measure_Nsqrd(Gup::AbstractMatrix{T}, Gdn::AbstractMatrix{T}) where {T}

Measure the expectation value of the total particle number squared $\langle \hat{N}^2 \rangle$ given both the spin-up and spin-down equal-time Green's function matrices $G_\uparrow(\tau,\tau)$ and $G_\downarrow(\tau,\tau)$ respectively.

fourier_transform!(
    C::AbstractArray{Complex{T}},
    a::Int,
    b::Int,
    dims,
    unit_cell::UnitCell{D,T},
    lattice::Lattice{D}
) where {D, T<:AbstractFloat}

fourier_transform!(C::AbstractArray{Complex{T}},
    a::Int,
    b::Int,
    unit_cell::UnitCell{D,T},
    lattice::Lattice{D}
) where {D, T<:AbstractFloat}

Calculate the fourier transform from position to momentum space

\[\begin{align*} C_{\mathbf{k}}^{a,b}= & \sum_{\mathbf{r}}e^{{\rm -i}\mathbf{k}\cdot(\mathbf{r}+\mathbf{r}_{a}-\mathbf{r}_{b})}C_{\mathbf{r}}^{a,b} \end{align*}\]

where $a$ and $b$ specify orbital species in the unit cell. Note that the array C is modified in-place. If dims is passed, iterate over these dimensions of the array, performing a fourier transform on each slice.

susceptibility!(χ::AbstractArray{T}, S::AbstractArray{T}, Δτ::E, dim::Int) where {T<:Number, E<:AbstractFloat}

Calculate the susceptibilities

\[\chi_\mathbf{n} = \int_0^\beta S_\mathbf{n}(\tau) d\tau,\]

where the $\chi_\mathbf{n}$ susceptibilities are written to χ, and S contains the $S_\mathbf{n}(\tau)$ correlations that need to be integrated over. The parameter Δτ is the discretization in imaginary time $\tau,$ and is the step size used in Simpson's method to numerically evaluate the integral over imaginary time. The argument dim specifies which dimension of S corresponds to imaginary time, and needs to be integrated over. Accordingly,

ndim(χ)+1 == ndim(S)

and

size(χ) == size(selectdim(S, dim, 1))

must both be true.

susceptibility(S::AbstractVector{T}, Δτ::E) where {T<:Number, E<:AbstractFloat}

Calculate the suceptibility

\[\chi = \int_0^\beta S(\tau) d\tau,\]

where the correlation data is stored in S. The integration is performed using Simpson's method using a step size of Δτ.

average_Gr0(
    G::AbstractMatrix{T},
    r::Bond{D}, unit_cell::UnitCell{D,E}, lattice::Lattice{D},
    sgn::T=one(T)
) where {D, T<:Number, E<:AbstractFloat}

Calculate the translationally averaged Green's function

G_\mathbf{r}^(a,b) = \frac{1}{N} \sum_{\mathbf{i}} G_{\mathbf{i}+\mathbf{r},\mathbf{i}}^{a,b}

defined by the bond r, which corresponds to the displacement $\mathbf{r} + (\mathbf{u}_a-\mathbf{u}_b).$

average_ηGr0(
    G::AbstractMatrix{T}, η::AbstractArray{T,D},
    r::Bond{D}, unit_cell::UnitCell{D,E}, lattice::Lattice{D},
    sgn::T=one(T)
) where {D, T<:Number, E<:AbstractFloat}

Calculate and return the sum

\[\frac{1}{N} \sum_{\mathbf{i}} \eta_{\mathbf{i}} G_{\mathbf{i}+\mathbf{r},\mathbf{i}}^{a,b},\]

where the bond r defines the displacement $\mathbf{r} + (\mathbf{u}_a-\mathbf{u}_b).$

contract_G00!(
    S::AbstractArray{C}, G::AbstractMatrix{T},
    a::Int, b::Int, α::Int,
    unit_cell::UnitCell{D,E}, lattice::Lattice{D},
    sgn::T=one(T)
) where {D, C<:Complex, T<:Number, E<:AbstractFloat}

Evaluate the sum

\[\begin{align*} S_{\mathbf{r}} := S_{\mathbf{r}} + \frac{\alpha}{N}\sum_{\mathbf{i}}G_{\sigma,\mathbf{i},\mathbf{i}}^{a,b}(\tau,0) \end{align*}\]

for all $\mathbf{r}.$

contract_Gr0!(
    S::AbstractArray{C,D}, G::AbstractMatrix{T},
    r′::Bond{D}, α::Int,
    unit_cell::UnitCell{D,E}, lattice::Lattice{D},
    sgn::T=one(T)
) where {D, C<:Complex, T<:Number, E<:AbstractFloat}

Evaluate the sum

\[S_{\mathbf{r}}:=S_{\mathbf{r}}+\frac{\alpha}{N}\sum_{\mathbf{i}} G_{\sigma,\mathbf{i}+\mathbf{r}+\mathbf{r}_{1},\mathbf{i}}^{a,b}(\tau,0)\]

for all $\mathbf{r},$ where the bond r′ represents the static displacement $\mathbf{r}_1+(\mathbf{r}_a-\mathbf{r}_b).$

contract_Gr0!(
    S::AbstractArray{C,D}, G::AbstractMatrix{T},
    η::AbstractArray{T}, r′::Bond{D}, α::Int,
    unit_cell::UnitCell{D,E}, lattice::Lattice{D},
    sgn::T=one(T)
) where {D, C<:Complex, T<:Number, E<:AbstractFloat}

Evaluate the sum

\[S_{\mathbf{r}}:=S_{\mathbf{r}}+\frac{\alpha}{N}\sum_{\mathbf{i}} \eta_{\mathbf{i}} G_{\sigma,\mathbf{i}+\mathbf{r}+\mathbf{r}_{1},\mathbf{i}}^{a,b}(\tau,0)\]

for all $\mathbf{r},$ where the bond r′ represents the static displacement $\mathbf{r}_1+(\mathbf{r}_a-\mathbf{r}_b).$

contract_Grr_G00!(
    S::AbstractArray{C,D}, G₂::AbstractMatrix{T}, G₁::AbstractMatrix{T},
    b₂::Bond{D}, b₁::Bond{D},
    α::Int, unit_cell::UnitCell{D,E}, lattice::Lattice{D},
    sgn::T=one(T)
) where {D, C<:Complex, T<:Number, E<:AbstractFloat}

Evaluate the sum

\[S_{\mathbf{r}}:=S_{\mathbf{r}}+\frac{\alpha}{N}\sum_{\mathbf{i}}G_{\sigma_{2},\mathbf{i}+\mathbf{r}+\mathbf{r}_{2},\mathbf{i}+\mathbf{r}}^{a,b}(\tau_{2},0)\cdot G_{\sigma_{1},\mathbf{i}+\mathbf{r}_{1},\mathbf{i}}^{c,d}(\tau_{1},0)\]

for all $\mathbf{r},$ where the bond b₂ represents the static displacement $\mathbf{r}_2 + (\mathbf{r}_a - \mathbf{r}_b),$ and the bond b₁ represents the static displacement $\mathbf{r}_1 + (\mathbf{r}_c - \mathbf{r}_d).$

contract_Grr_G00!(
    S::AbstractArray{C,D}, G₂::AbstractMatrix{T}, G₁::AbstractMatrix{T},
    η₂::AbstractArray{T}, η₁::AbstractArray{T}, b₂::Bond{D}, b₁::Bond{D},
    α::Int, unit_cell::UnitCell{D,E}, lattice::Lattice{D},
    sgn::T=one(T)
) where {D, C<:Complex, T<:Number, E<:AbstractFloat}

Evaluate the sum

\[S_{\mathbf{r}}:=S_{\mathbf{r}}+\frac{\alpha}{N}\sum_{\mathbf{i}} \eta_2{\mathbf{i}+\mathbf{r}}\eta_1{\mathbf{i}} G_{\sigma_{2},\mathbf{i}+\mathbf{r}+\mathbf{r}_{2},\mathbf{i}+\mathbf{r}}^{a,b}(\tau_{2},0)\cdot G_{\sigma_{1},\mathbf{i}+\mathbf{r}_{1},\mathbf{i}}^{c,d}(\tau_{1},0)\]

for all $\mathbf{r},$ where the bond b₂ represents the static displacement $\mathbf{r}_2 + (\mathbf{r}_a - \mathbf{r}_b),$ and the bond b₁ represents the static displacement $\mathbf{r}_1 + (\mathbf{r}_c - \mathbf{r}_d).$

contract_Grr_G00!(
    S::AbstractArray{C,D}, G₂::AbstractMatrix{T}, G₁::AbstractMatrix{T},
    a::Int, b::Int, c::Int, d::Int,
    r₄::AbstractVector{Int}, r₃::AbstractVector{Int}, r₂::AbstractVector{Int}, r₁::AbstractVector{Int},
    α::Int, unit_cell::UnitCell{D,E}, lattice::Lattice{D},
    sgn::T=one(T)
) where {D, C<:Complex, T<:Number, E<:AbstractFloat}

Evaluate the sum

\[S_{\mathbf{r}}:=S_{\mathbf{r}}+\frac{\alpha}{N}\sum_{\mathbf{i}} G_{\sigma_{2},\mathbf{i}+\mathbf{r}+\mathbf{r}_{4},\mathbf{i}+\mathbf{r}+\mathbf{r}_{3}}^{a,b}(\tau_{2},0)\cdot G_{\sigma_{1},\mathbf{i}+\mathbf{r}_{2},\mathbf{i}+\mathbf{r}_{1}}^{c,d}(\tau_{1},0)\]

for all $\mathbf{r}.$

contract_Grr_G00!(
    S::AbstractArray{C,D}, G₂::AbstractMatrix{T}, G₁::AbstractMatrix{T},
    η₂::AbstractArray{T}, η₁::AbstractArray{T}, a::Int, b::Int, c::Int, d::Int,
    r₄::AbstractVector{Int}, r₃::AbstractVector{Int}, r₂::AbstractVector{Int}, r₁::AbstractVector{Int},
    α::Int, unit_cell::UnitCell{D,E}, lattice::Lattice{D},
    sgn::T=one(T)
) where {D, C<:Complex, T<:Number, E<:AbstractFloat}

Evaluate the sum

\[S_{\mathbf{r}}:=S_{\mathbf{r}}+\frac{\alpha}{N}\sum_{\mathbf{i}} \eta_2{\mathbf{i}+\mathbf{r}}\eta_1{\mathbf{i}} G_{\sigma_{2},\mathbf{i}+\mathbf{r}+\mathbf{r}_{4},\mathbf{i}+\mathbf{r}+\mathbf{r}_{3}}^{a,b}(\tau_{2},0)\cdot G_{\sigma_{1},\mathbf{i}+\mathbf{r}_{2},\mathbf{i}+\mathbf{r}_{1}}^{c,d}(\tau_{1},0)\]

for all $\mathbf{r}.$

contract_Gr0_Gr0!(
    S::AbstractArray{C,D}, G₂::AbstractMatrix{T}, G₁::AbstractMatrix{T}, b₂::Bond, b₁::Bond,
    α::Int, unit_cell::UnitCell{D,E}, lattice::Lattice{D},
    sgn::T=one(T)
) where {D, C<:Complex, T<:Number, E<:AbstractFloat}

Evaluate the sum

\[S_{\mathbf{r}}:=S_{\mathbf{r}}+\frac{\alpha}{N}\sum_{\mathbf{i}}G_{\sigma_{2},\mathbf{i}+\mathbf{r}+\mathbf{r}_{2},\mathbf{i}+\mathbf{r}_{1}}^{a,c}(\tau_{2},0)\cdot G_{\sigma_{1},\mathbf{i}+\mathbf{r},\mathbf{i}}^{b,d}(\tau_{1},0)\]

for all $\mathbf{r},$ where the bond b₂ represents the static displacement $\mathbf{r}_2 + (\mathbf{r}_a - \mathbf{r}_b),$ and the bond b₁ represents the static displacement $\mathbf{r}_1 + (\mathbf{r}_c - \mathbf{r}_d).$

contract_Gr0_Gr0!(
    S::AbstractArray{C,D}, G₂::AbstractMatrix{T}, G₁::AbstractMatrix{T},
    η₂::AbstractArray{T}, η₁::AbstractArray{T}, b₂::Bond, b₁::Bond,
    α::Int, unit_cell::UnitCell{D,E}, lattice::Lattice{D},
    sgn::T=one(T)
) where {D, C<:Complex, T<:Number, E<:AbstractFloat}

Evaluate the sum

\[S_{\mathbf{r}}:=S_{\mathbf{r}}+\frac{\alpha}{N}\sum_{\mathbf{i}} \eta_2(\mathbf{i}+\mathbf{r})\eta_1(\mathbf{i}) G_{\sigma_{2},\mathbf{i}+\mathbf{r}+\mathbf{r}_{2},\mathbf{i}+\mathbf{r}_{1}}^{a,c}(\tau_{2},0)\cdot G_{\sigma_{1},\mathbf{i}+\mathbf{r},\mathbf{i}}^{b,d}(\tau_{1},0)\]

for all $\mathbf{r},$ where the bond b₂ represents the static displacement $\mathbf{r}_2 + (\mathbf{r}_a - \mathbf{r}_b),$ and the bond b₁ represents the static displacement $\mathbf{r}_1 + (\mathbf{r}_c - \mathbf{r}_d).$

contract_Gr0_Gr0!(
    S::AbstractArray{C,D}, G₂::AbstractMatrix{T}, G₁::AbstractMatrix{T},
    a::Int, b::Int, c::Int, d::Int,
    r₄::AbstractVector{Int}, r₃::AbstractVector{Int}, r₂::AbstractVector{Int}, r₁::AbstractVector{Int},
    α::Int, unit_cell::UnitCell{D,E}, lattice::Lattice{D},
    sgn::T=one(T)
) where {D, C<:Complex, T<:Number, E<:AbstractFloat}

Evaluate the sum

\[S_{\mathbf{r}}:=S_{\mathbf{r}}+\frac{\alpha}{N}\sum_{\mathbf{i}}G_{\sigma_{2},\mathbf{i}+\mathbf{r}+\mathbf{r}_{4},\mathbf{i}+\mathbf{r}_{3}}^{a,b}(\tau_{2},0)\cdot G_{\sigma_{1},\mathbf{i}+\mathbf{r}+\mathbf{r}_{2},\mathbf{i}+\mathbf{r}_{1}}^{c,d}(\tau_{1},0)\]

for all $\mathbf{r}.$

contract_Gr0_Gr0!(
    S::AbstractArray{C,D}, G₂::AbstractMatrix{T}, G₁::AbstractMatrix{T},
    η₂::AbstractArray{T}, η₁::AbstractArray{T}, a::Int, b::Int, c::Int, d::Int,
    r₄::AbstractVector{Int}, r₃::AbstractVector{Int}, r₂::AbstractVector{Int}, r₁::AbstractVector{Int},
    α::Int, unit_cell::UnitCell{D,E}, lattice::Lattice{D},
    sgn::T=one(T)
) where {D, C<:Complex, T<:Number, E<:AbstractFloat}

Evaluate the sum

\[S_{\mathbf{r}}:=S_{\mathbf{r}}+\frac{\alpha}{N}\sum_{\mathbf{i}} \eta_2(\mathbf{i}+\mathbf{r})\eta_1(\mathbf{i}) G_{\sigma_{2},\mathbf{i}+\mathbf{r}+\mathbf{r}_{4},\mathbf{i}+\mathbf{r}_{3}}^{a,b}(\tau_{2},0)\cdot G_{\sigma_{1},\mathbf{i}+\mathbf{r}+\mathbf{r}_{2},\mathbf{i}+\mathbf{r}_{1}}^{c,d}(\tau_{1},0)\]

for all $\mathbf{r}.$

contract_G0r_Gr0!(
    S::AbstractArray{C,D}, G₂::AbstractMatrix{T}, G₁::AbstractMatrix{T},
    b₂::Bond, b₁::Bond, α::Int, unit_cell::UnitCell{D,E}, lattice::Lattice{D},
    sgn::T=one(T)
) where {D, C<:Complex, T<:Number, E<:AbstractFloat}

Evaluate the sum

\[S_{\mathbf{r}}:=S_{\mathbf{r}}+\frac{\alpha}{N}\sum_{\mathbf{i}}G_{\sigma_{2},\mathbf{i}+\mathbf{r}_{2},\mathbf{i}+\mathbf{r}}^{a,b}(\tau_{2},0)\cdot G_{\sigma_{1},\mathbf{i}+\mathbf{r}+\mathbf{r}_{1},\mathbf{i}}^{c,d}(\tau_{1},0)\]

for all $\mathbf{r},$ where the bond b₂ represents the static displacement $\mathbf{r}_2 + (\mathbf{r}_a - \mathbf{r}_b),$ and the bond b₁ represents the static displacement $\mathbf{r}_1 + (\mathbf{r}_c - \mathbf{r}_d).$

contract_G0r_Gr0!(
    S::AbstractArray{C,D}, G₂::AbstractMatrix{T}, G₁::AbstractMatrix{T},
    η₂::AbstractArray{T}, η₁::AbstractArray{T}, b₂::Bond, b₁::Bond,
    α::Int, unit_cell::UnitCell{D,E}, lattice::Lattice{D},
    sgn::T=one(T)
) where {D, C<:Complex, T<:Number, E<:AbstractFloat}

Evaluate the sum

\[S_{\mathbf{r}}:=S_{\mathbf{r}}+\frac{\alpha}{N}\sum_{\mathbf{i}} \eta_2(\mathbf{i}+\mathbf{r})\eta_1(\mathbf{i}) G_{\sigma_{2},\mathbf{i}+\mathbf{r}_{2},\mathbf{i}+\mathbf{r}}^{a,b}(\tau_{2},0)\cdot G_{\sigma_{1},\mathbf{i}+\mathbf{r}+\mathbf{r}_{1},\mathbf{i}}^{c,d}(\tau_{1},0)\]

for all $\mathbf{r},$ where the bond b₂ represents the static displacement $\mathbf{r}_2 + (\mathbf{r}_a - \mathbf{r}_b),$ and the bond b₁ represents the static displacement $\mathbf{r}_1 + (\mathbf{r}_c - \mathbf{r}_d).$

contract_G0r_Gr0!(
    S::AbstractArray{C,D}, G₂::AbstractMatrix{T}, G₁::AbstractMatrix{T},
    a::Int, b::Int, c::Int, d::Int,
    r₄::AbstractVector{Int}, r₃::AbstractVector{Int}, r₂::AbstractVector{Int}, r₁::AbstractVector{Int},
    α::Int, unit_cell::UnitCell{D,E}, lattice::Lattice{D},
    sgn::T=one(T)
) where {D, C<:Complex, T<:Number, E<:AbstractFloat}

Evaluate the sum

\[S_{\mathbf{r}}:=S_{\mathbf{r}}+\frac{\alpha}{N}\sum_{\mathbf{i}}G_{\sigma_{2},\mathbf{i}+\mathbf{r}_{4},\mathbf{i}+\mathbf{r}+\mathbf{r}_{3}}^{a,b}(\tau_{2},0)\cdot G_{\sigma_{1},\mathbf{i}+\mathbf{r}+\mathbf{r}_{2},\mathbf{i}+\mathbf{r}_{1}}^{c,d}(\tau_{1},0)\]

for all $\mathbf{r}.$

contract_G0r_Gr0!(
    S::AbstractArray{C,D}, G₂::AbstractMatrix{T}, G₁::AbstractMatrix{T},
    η₂::AbstractArray{T}, η₁::AbstractArray{T}, a::Int, b::Int, c::Int, d::Int,
    r₄::AbstractVector{Int}, r₃::AbstractVector{Int}, r₂::AbstractVector{Int}, r₁::AbstractVector{Int},
    α::Int, unit_cell::UnitCell{D,E}, lattice::Lattice{D},
    sgn::T=one(T)
) where {D, C<:Complex, T<:Number, E<:AbstractFloat}

Evaluate the sum

\[S_{\mathbf{r}}:=S_{\mathbf{r}}+\frac{\alpha}{N}\sum_{\mathbf{i}} \eta_2(\mathbf{i}+\mathbf{r})\eta_1(\mathbf{i}) G_{\sigma_{2},\mathbf{i}+\mathbf{r}_{4},\mathbf{i}+\mathbf{r}+\mathbf{r}_{3}}^{a,b}(\tau_{2},0)\cdot G_{\sigma_{1},\mathbf{i}+\mathbf{r}+\mathbf{r}_{2},\mathbf{i}+\mathbf{r}_{1}}^{c,d}(\tau_{1},0)\]

for all $\mathbf{r}.$

simpson(f::AbstractVector{T}, dx::E) where {T<:Number, E<:AbstractFloat}

Applying Simpson's rule, integrate over the vector f using a stepsize dx.