Brillouinzone integration in scuffem
Many codes in the scuffem suite require evaluating integrals over the Brillouin zone (BZ) of a 1D or 2D reciprocal lattice, i.e.
where we generically use the overlined symbol to denote the contribution of Bloch vector to quantity . Examples of calculations that require Brillouinzone integrations include

the Casimir force per unit imaginary frequency on an extended object in scuffcas3d

the CasimirPolder potential per unit imaginary frequency experienced by a polarizable particle near an extended surface in scuffcaspol

the local density of states at a given angular frequency at userspecified evaluation points in scuffldos.
In general, Brillouinzone integrations are evaluated by numerical cubaturethat is, as weighted sums of integrand samples: where are the weights and points in a cubature rule for the Brillouin zone of your reciprocal lattice, and where each integrand sample is computed by performing a single scuffem calculation at a fixed Bloch wavevector. The scuffem workflow offers two options for evaluating such cubatures:
 You can design and implement your own cubature scheme
involving your own customchosen weights and points
. In this case, you will use
the
byOmegakBloch
commandline option to instruct a scuffem application code to report values of the quantity at each of your points (this output will typically be written to file with extension.byOmegakBloch
or.byXikBloch
), then compute the weighted sums yourself in e.g. julia.
 Alternatively, you can ask scuffem to perform the
BZ integration internally, using one of several
builtin cubature schemes. In this case the BZintegrated
quantities will typically be written to
an output file with extension
.byOmega
or.byXi
. You will also get an output file named.byOmegakBloch
or.byXikBloch
that reports the Blochvectorresolved integrand samples chosen internally by the scuffem BZ integrator.
Commandline options for customizing internal BZ integration
If you choose the second option above, you may specify various commandline options to customize the algorithm used by scuffem to select the cubature points and weights . The options are listed here and discussed in more detail below.
BZIMethod [CC  TC  Polar  Polar2]
Selects the integration algorithm (see below for details).
BZIOrder NN
Sets the order (accuracy parameter) of the integration
algorithm to NN
. The allowed values of NN
here depend
on the integration algorithm you chose (see below).
BZIRelTol xx
BZIAbsTol xx
BZIMaxEvals NN
For adaptive integration algorithms in which the order is determined internally (see below), these options allow you to specify relative and absolute error tolerances and an upper limit on the number of integrand samples that will be used.
BZSymmetryFactor [248]
This option lets you tell scuffem that your integrand function is invariant under 2, 4, or 8fold rotational symmetry transformations applied to . See below for more details on what this means.
Understanding the internal BZ integration algorithms
To help you understand how to configure the various commandline options above, this section describes the various integration algorithms available and how they are affected by the commandline parameters.
Integration methods for 1D Brillouin zones
For onedimensional Brillouin zones, there is only one
BZ integration method availablenamely,
ClenshawCurtis quadrature
(BZIMethod CC
, the default)
with the number of integrand samples either fixed
or chosen adaptively until userspecified error tolerances
are achieved.
More specifically, you may say either
BZIOrder [11  13  15  ...  97  99]
or
BZIOrder 0
The former option selects fixedorder CC cubature with 11, 13, ... 99 sample points. (This number must be an odd integer between 11 and 99 inclusive.)
The latter option selects adaptive
CC cubature
using this algorithm.
In this case the number of sample points will be chosen
automatically subject to the values you select for the
BZIRelTol
, BZIAbsTol
, and BZIMaxEvals
commandline options.
Symmetry factors for 1D Brillouin zones
For 1D Brillouin zones, the only allowed
value of the BZSymmetryFactor
option is 2
,
indicating that your integrand is symmetric under
sign flip of the Bloch wavevector, i.e.
.
In this case the BZ integration may be restricted to
the range
Integration methods for 2D Brillouin zones
The following integration methods are implemented for 2D Brillouin zones. (See below for pictures of where the various methods place their sample points.)

ClenshawCurtis cubature (
BZIMethod CC
)Nested 2D fixedorder or adaptive ClenshawCurtis cubature.
For fixedorder nested CC cubature with
NN
sample points per dimension, sayBZIOrder NN
. HereNN
must be an odd integer between 11 and 99 inclusive.For adaptive 2D CC cubature (subject to your specified values of
BZIRelTol
,BZIAbsTol
, andBZIMaxEvals
) sayBZIOrder 0.

Triangle cubature (
BZIMethod TC
)This algorithm divides the Brillouin zone into 8 triangles and applies a fixedorder triangle cubature scheme to each triangle, omitting repetition of triangles that are symmetryequivalent given the value you specified for
BZSymmetryFactor
(see pictures below).For this algorithm, the allowed values of
BZIOrder
are 1, 2, 4, 5, 7, 9, 13, 14, 16, 20, or 25.

Polar cubature (
BZIMethod Polar
)This algorithm uses a polar decomposition to evaluate the BZ integral as two nested 1D integrals, one over and the second over . (The integral is evaluated via ClenshawCurtis quadrature, and the quadrature is evaluated using rectangularrule quadrature, not necessarily of the same order).
This algorithm is useful for integrands that are strongly peaked near the origin of the Brillouin zone and highly attenuated near the boundaries, so that most of the integral comes from the region near the origin.
For this algorithm, the number of integration points used for the and integrals (call these numbers and ) are encoded into the value passed to
BZIOrder
in the form where and are each odd integers between 11 and 99 inclusive.Thus, for example,
BZIOrder 3321
specifies that the integral is to be evaluated via 33point CC cubature, while the integral is to be evaluated via 21point rectangularrule cubature.

Polar cubature with change of variables (
BZIMethod Polar2
)This is the same as
BZIMethod Polar
, but with two modifications: (a) The integral is split into two separate integrals covering the ranges and (where is the vacuum photon wavenumber at the frequency in question). (b) In each of the two integrals we make the change of variablesThese modifications are useful for cases in which the freespace wavevector falls within the Brillouin zone. In these cases, the convergence of the integral is degraded by the phenomena known as "Wood anomalies" in optics or "van Hove singularities" in solidstate physics, and changing variables to introduces a Jacobian factor that neutralizes these singularities to yield a betterbehaved integrand.
Special values for rotationallyinvariant integrands
As discussed above, for the Polar
and Polar2
integration methods
the value of the BZIOrder
option is interpreted as the composite
quantity ,
where and are odd integers in the range
specifying the number of cubature points used for the and
integrals (or set to 0 to request adaptive quadrature).
For fully rotationallysymmetric geometries in which the BZ integrand is independent of , you can specify to indicate that the integral is to be evaluated by a 1point cubature with the single sample taken at .
(Of course, no geometry discretized into triangles can actually be fully rotationally invariant, but pretending so may be a reasonable approximation in some cases, such as this one.)
Symmetry factors for 2D Brillouin zones
For 2D geometries, the option BZSymmetryFactor
may
take the value 2, 4, or 8, specifying that the Brillouinzone
integrand obeys symmetries as follows:

BZSymmetryFactor 2
:We have so the BZ integration may be restricted to just the right half of the BZ.

BZSymmetryFactor 4
:We have so the BZ integration may be restricted to just the upperright quadrant of the BZ.
BZSymmetryFactor 8
: In addition to symmetry under sign changes, the integrand is symmetric under , so the BZ integration may be restricted to the triangular region (This is only possible for square lattices.)
Locations of quadrature points for 2D Brillouin zones
Here are some diagrams indicating the Bloch wavevectors that will be sampled by the internal algorithms for Brillouinzone integration with various values of the commandline parameters:
BZIMethod CC BZIOrder 21 BZSymmetryFactor 1
BZIMethod CC BZIOrder 21 BZSymmetryFactor 2
BZIMethod CC BZIOrder 21 BZSymmetryFactor 4
BZIMethod CC BZIOrder 21 BZSymmetryFactor 8
BZIMethod TC BZIOrder 20 BZSymmetryFactor 1
BZIMethod TC BZIOrder 20 BZSymmetryFactor 2
BZIMethod TC BZIOrder 20 BZSymmetryFactor 4
BZIMethod TC BZIOrder 20 BZSymmetryFactor 8
BZIMethod Polar BZIOrder 3111 BZSymmetryFactor 1
BZIMethod Polar BZIOrder 3111 BZSymmetryFactor 2
BZIMethod Polar BZIOrder 3111 BZSymmetryFactor 4
BZIMethod Polar BZIOrder 3111 BZSymmetryFactor 8
Note: The next four diagrams are for a frequency of
. The accumulation
of points near
is noticeable for the Polar2
integration strategy.
BZIMethod Polar2 BZIOrder 1111 BZSymmetryFactor 1
BZIMethod Polar2 BZIOrder 1111 BZSymmetryFactor 2
BZIMethod Polar2 BZIOrder 1111 BZSymmetryFactor 4
BZIMethod Polar2 BZIOrder 1111 BZSymmetryFactor 8
Locations of quadrature points for 1D Brillouin zones
BZIMethod CC BZIOrder 31 BZSymmetryFactor 1
BZIMethod CC BZIOrder 31 BZSymmetryFactor 2