# Examples

The following are several examples of how to use the algorithms and integrands provided by AutoBZCore.jl. For background on defining integrals see the Problem definitions page

## Green's function

A common integral appearing in Dynamical mean-field theory is that of the local Green's function:

\[G(\omega) = \int d^d \vec{k}\ \operatorname{Tr} \left[ \left( \omega - H \left( \vec{k} \right) - \Sigma(\omega) \right)^{-1} \right].\]

For simplicity, we take $\Sigma(\omega) = -i\eta$. We can define the integrand as a function of $\vec{k}$ and $H$ and (required) parameters $\eta, \omega$.

```
using LinearAlgebra
gloc_integrand(k, h; η, ω) = inv(complex(ω,η)*I-h(k))
```

Commonly, $H(\vec{k})$ is evaluated using Wannier interpolation, i.e. as a Fourier series. For a simple tight-binding model, the integer lattice, the Hamiltonian is given by

\[H(k) = \cos(2\pi k) = \frac{1}{2} \left( e^{2\pi ik} + e^{-2\pi ik} \right)\]

We can use the built-in function `cos`

to evaluate this, however, for more complex Fourier series it becomes easier to use the representation in terms of Fourier coefficients. Using the package FourierSeriesEvaluators.jl, we can define $H(k) = \cos(2\pi k)$ by the following:

```
using FourierSeriesEvaluators
h = FourierSeries([0.5, 0.0, 0.5]; period=1, offset=-2)
```

The coefficient values of $1/2$ can be determined from Euler's formula, as used in the expansion of $cos$ above, and the value of `offset`

is chosen to offset the coefficient array indices, `1:3`

since Julia has 1-based indexing, to correspond to values of $n$ in the phase factors $e^{2\pi i n k}$ used in the Fourier series above, i.e. `-1:1`

. Now we proceed to the define the integral problem

```
using AutoBZCore
using AutoBZCore: IntegralProblem
integrand = ParameterIntegrand(gloc_integrand, h, η=0.1)
prob = IntegralProblem(integrand, 0, 1)
```

Here we wrapped our function with two of its arguments, `h, η`

as a `ParameterIntegrand`

that allows us to provide partial arguments so that we can solve the integral as a function of the remaining parameters, in this case `ω`

. We also created an `AutoBZCore.IntegralProblem`

to integrate our function over its period $[0,1]$. To solve this problem, we pick any of the package's Integral algorithms and the tolerance to which we would like the solution. Then we make an `IntegralSolver`

to evaluate $G(\omega)$ as a function.

```
alg = QuadGKJL()
gloc = IntegralSolver(prob, alg; abstol=1e-3)
gloc(ω=0.0) # -2.7755575615628914e-17 - 0.9950375451895513im
```

Now we proceed to multi-dimensional integrals. In this case, Wannier interpolation is much more efficient when Fourier series are evaluated one variable at a time. To understand, this suppose we have a series defined by $M \times M$ coefficients (i.e. a 2d series) that we want to evaluate on an $N \times N$ grid. Naively evaluating the series at each grid point will require $\mathcal{O}(M^{2} N^{2})$ operations, however, we can reduce the complexity by pre-evaluating certain coefficients as follows

\[f(x, y) = \sum_{m,n=1}^{M} f_{nm} e^{i(nx + my)} = \sum_{n=1}^{M} e^{inx} \left( \sum_{m=1}^{M} f_{nm} e^{imy} \right) = \sum_{n=1}^{M} e^{inx} \tilde{f}_{n}(y)\]

This means we can evaluate the series on the grid in $\mathcal{O}(M N^2 + M^2 N)$ operations. When $N \gg M$, this is $\mathcal{O}(M N^{2})$ operations, which is comparable to the computational complexity of a multi-dimensional FFT. Since the constants of a FFT may not be trivial, this scheme is competitive.

This capability is provided by `FourierIntegrand`

. Let's use this with a Fourier series corresponding to $H(\vec{k}) = \cos(2\pi k_{x}) + \cos(2\pi k_{y})$

```
h = FourierSeries([0.0; 0.5; 0.0;; 0.5; 0.0; 0.5;; 0.0; 0.5; 0.0]; period=1, offset=-2)
integrand = FourierIntegrand(gloc_integrand, h, η=0.1)
```

However, since `FourierIntegrand`

evaluates $H(k)$ for us and gives it as a `FourierValue`

together with $k$, we need to define another method for our integrand to comply with the interface

`gloc_integrand(h_k::FourierValue; η, ω) = inv(complex(ω,η)*I-h_k.s)`

Similar to before, we construct an `AutoBZCore.IntegralProblem`

and this time we take the integration domain to correspond to the full Brillouin zone of a square lattice with lattice vectors `2pi*I(2)`

. (See the Reference for more details on constructing BZs.)

```
bz = load_bz(FBZ(2), 2pi*I(2))
prob = IntegralProblem(integrand, bz)
```

This package provides several BZ-specific integral algorithms that we can use to solve the multidimensional integral.

```
alg = IAI()
gloc = IntegralSolver(prob, alg; abstol=1e-3)
gloc(ω=0.0) # 1.5265566588595902e-16 - 1.3941704019631334im
```

## Density of states

The repo's demo on density of states provides a complete example of how to compute and interpolate an integral as a function of its parameters.