API

greens!(
    G::AbstractArray{C,D},
    a::Int, b::Int, unit_cell::UnitCell{D}, lattice::Lattice{D},
    G_τ0::AbstractMatrix{T},
    sgn=one(C)
) 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=one(C): 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=one(C)
) 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=one(C): 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=one(C)
) 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=one(C): 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=one(C)
) 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}= & \frac{1}{2}(\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)\\ = & \frac{1}{2}\left(\hat{a}_{\uparrow,\mathbf{i}}^{\dagger}\hat{a}_{\downarrow,\mathbf{i}}+\hat{a}_{\downarrow,\mathbf{i}}^{\dagger}\hat{a}_{\uparrow,\mathbf{i}}\right). \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=one(C): 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=one(C)
) 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}}= & frac{1}{2} (\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)\\ = & \frac{1}{2} \left( \hat{n}_{\uparrow,a,\mathbf{i}}-\hat{n}_{\downarrow,a,\mathbf{i}} \right). \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=one(C): 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=one(C)
) 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=one(C): 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=one(C)
) 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=one(C): 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=one(C)
) 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=one(C): 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=one(C)
) 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 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}} - t_{\sigma,\mathbf{i}, \mathbf{i}+\mathbf{r}}^{b,a} \hat{b}_{\sigma,\mathbf{i}}^{\dagger}\hat{a}_{\sigma,\mathbf{i}+\mathbf{r}})\\ \end{align*}\]

and $t_{\sigma,\mathbf{i}, \mathbf{i}+\mathbf{r}}^{b,a} = \big( t_{\sigma,\mathbf{i}+\mathbf{r},\mathbf{i}}^{a,b} \big)^*$.

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=one(C): 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=one(C)
) 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}} - t_{\sigma,\mathbf{i}, \mathbf{i}+\mathbf{r}}^{b,a} \hat{b}_{\sigma,\mathbf{i}}^{\dagger}\hat{a}_{\sigma,\mathbf{i}+\mathbf{r}})\\ \end{align*}\]

and $t_{\sigma,\mathbf{i}, \mathbf{i}+\mathbf{r}}^{b,a} = \big( t_{\sigma,\mathbf{i}+\mathbf{r},\mathbf{i}}^{a,b} \big)^*$.

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=one(C): 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.

cubic_spline_τ_to_ωn!(
    Cn::AbstractVector{Complex{E}},
    Cτ::AbstractVector{T},
    β::E, Δτ::E;
    spline_type::String = "C2",
    M1::T = NaN,
    M2::T = NaN
) where {E<:AbstractFloat, T<:Number}

Calculate the Matsubara frequency representation

\[\begin{align*} C(\mathrm{i}\omega_n) = \int_0^\beta d\tau \ e^{\mathrm{i}\omega_n\tau} C(\tau) \end{align*}\]

of an imaginary-time correlation function $C(\tau)$ defined on a regular grid of $L + 1$ imaginary-time points $\tau \in \{ 0, \Delta\tau, 2\Delta\tau, \ldots, (\beta-\Delta\tau), \beta \}.$ This is done by fitting a cubic spline through the $C(\tau)$ points, and then transforming the spline. Here, the Matsubara correlations $C(\mathrm{i}\omega_n)$ are written to the vector Cn, and the imaginary-time correlations $C(\tau)$ are stored in the vector Cτ. It is assumed that C[1] and C[end] correspond to $C(\tau = 0)$ and $C(\tau = \beta)$ respectively. It is also assumed that $L = \beta/\Delta\tau.$

Note that the length of the vector Cn determines whether the correlation function is assumed to be fermionic or bosonic.

For fermoinic correlation functions, mod(M, 2L) == 0, where M = length(Cn) and L = length(Cτ) - 1. Put another way, it must be that length(Cn) == 2*n*(length(Cτ)-1) for some positive integer n ≥ 1. In this case the vector Cn contains Matsubara frequency correlation $C(\mathrm{i}\omega_n)$ for $n \in [-N, (N-1)]$, where N = M÷2 and $\omega_n = (2n+1)\pi/\beta$.

For bosonic correlation functions mod(M+1, 2L) == 0, where M = length(Cn) and L = length(Cτ) - 1. Put another way, it must be that length(Cn) == 2*n*(length(Cτ)-1) - 1 for some positive integer n ≥ 1. In this case the vector Cn contains Matsubara frequency correlation $C(\mathrm{i}\omega_n)$ for $n \in [-N, N]$, where N = (M-1)÷2 and $\omega_n = 2n\pi/\beta$.

The spline_type argument can be set to "C2", "akima", or "makima" to specify the type of spline to use. The default argument (and best for most cases) is "C2", which refers to a C2 cubic spline with continuous first and second derivatives. In this case, boundary conditions are imposed such that the second derivative of the spline at τ=0 and τ=β are set equal to the second derivative as calculated using a second-order forward and backward finite difference respectively.

The M1 and M2 arguments can be used to specify the first and second moments of the spectral function if known. By default the are set to NaN, which results in them being calculated internally based on the imaginary-time correlation function data points and corresponding spline fit.

jackknife_cubic_spline_τ_to_ωn!(
    Cn::AbstractVector{Complex{E}},
    Cτ::AbstractMatrix{T},
    β::E, Δτ::E;
    # KEYWORD ARGUMENTS
    spline_type::String = "C2",
    M1::T = NaN,
    M2::T = NaN,
    bias_corrected = false,
    return_covariance_matrix = false
) where {E<:AbstractFloat, T<:Number}

Perform a jackknife error analysis on the cubic spline Matsubara transformation. This function takes as input a matrix Cτ of imaginary-time correlation functions, where each column corresponds to a separate sample, and the rows correspond to the imaginary-time points. Refer to cubic_spline_τ_to_ωn! for details on the Matsubara transformation itself, and more information on the meaning of the spline_type, M1, and M2 keyword arguments. In addition to filling in the vector Cn with the mean Matsubara correlation function, this function also returns the corresponding standard error for each Cn value if return_covariance_matrix = false, and the full covariance matrix if return_covariance_matrix = true.

fourier_transform!(
    C::AbstractArray{Complex{T}},
    a::Int,
    b::Int,
    dims,
    unit_cell::UnitCell{D,T},
    lattice::Lattice{D},
    fftplan!::Union{F, Nothing} = nothing
) where {D, T<:AbstractFloat, F<:AbstractFFTs.Plan}

fourier_transform!(
    C::AbstractArray{Complex{T}},
    r::AbstractVector{T},
    dims,
    unit_cell::UnitCell{D,T},
    lattice::Lattice{D},
    fftplan!::Union{F, Nothing} = nothing
) where {D, T<:AbstractFloat, F<:AbstractFFTs.Plan}

fourier_transform!(
    C::AbstractArray{Complex{T},D},
    a::Int,
    b::Int,
    unit_cell::UnitCell{D,T},
    lattice::Lattice{D},
    fftplan!::Union{F, Nothing} = nothing
) where {D, T<:AbstractFloat, F<:AbstractFFTs.Plan}

fourier_transform!(
    C::AbstractArray{Complex{T},D},
    r::AbstractVector{T},
    unit_cell::UnitCell{D,T},
    lattice::Lattice{D},
    fftplan!::Union{F, Nothing} = nothing
) where {D, T<:AbstractFloat, F<:AbstractFFTs.Plan}

Calculate the fourier transform from position to momentum space

\[\begin{align*} C_{\mathbf{K}} = & \sum_{\mathbf{R}}e^{{\rm -i}\mathbf{K}\cdot(\mathbf{R}+\mathbf{r})}C_{\mathbf{R}} \end{align*}\]

where $\mathbf{r}$ is a constant displacement vector. If orbitals $a$ and $b$ are passed, then $\mathbf{r} = \mathbf{r}_a - \mathbf{r}_b$ where $\mathbf{r}_a$ and $\mathbf{r}_b$ are the basis vectors for each orbital 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. It is also possible to optionally pass an FFT plan to accelerate the fourier transformation, though it need to be an FFT plan that operates in-place on an array of complexes with the dimension of C.

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 Δτ.

jackknife(
    # ARGUMENTS
    g::Function,
    samples...;
    # KEYWORD ARGUMENTS
    bias_corrected = false,
    jackknife_sample_means = similar.(samples),
    jackknife_g = similar(first(samples)),
)

Propagate errors through the evaluation of a function g given the binned samples, returning both the mean and standard error as a tuple (mean, std_error).

Arguments

• g::Function: Function to evaluate, taking the same number of arguments as samples
• samples...: Vectors of individual observations (binned samples)

Keyword Arguments

• bias_corrected = false: If true, apply $\mathcal{O}(1/N)$ bias correction
• jackknife_sample_means: Preallocated arrays to avoid temporary allocations.
• jackknife_g: Preallocated array for jackknife function evaluations to avoid temporary allocations.

Returns

• mean: By default, the full sample estimate g(mean.(samples)...). If bias_corrected=true, returns the bias-corrected estimate.
• std_error: Jackknife standard error estimate.

Notes

The standard error is computed from the jackknife replicates. By default, the returned mean is the full sample estimate rather than the jackknife mean, as this is typically the best point estimate from the available data.

average_Gr0(
    G::AbstractMatrix{T},
    r::Bond{D}, unit_cell::UnitCell{D,E}, lattice::Lattice{D},
    sgn=one(C)
) 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=one(C)
) 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=one(C)
) 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=one(C)
) 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=one(C)
) 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=one(C)
) 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=one(C), conj_η₂::Bool = false, conj_η₁::Bool = false
) 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=one(C)
) 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=one(C), conj_η₂::Bool = false, conj_η₁::Bool = false
) 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=one(C)
) 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=one(C), conj_η₂::Bool = false, conj_η₁::Bool = false
) 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=one(C)
) 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=one(C), conj_η₂::Bool = false, conj_η₁::Bool = false
) 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=one(C)
) 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=one(C), conj_η₂::Bool = false, conj_η₁::Bool = false
) 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=one(C)
) 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}.$

function 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=one(C), conj_η₂::Bool = false, conj_η₁::Bool = false
) 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.