Brillouin-zone integration in scuff-em
Many codes in the scuff-em 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 Brillouin-zone integrations include
the Casimir force per unit imaginary frequency on an extended object in scuff-cas3d
the Casimir-Polder potential per unit imaginary frequency experienced by a polarizable particle near an extended surface in scuff-caspol
the local density of states at a given angular frequency at user-specified evaluation points in scuff-ldos.
In general, Brillouin-zone integrations are evaluated by numerical cubature---that 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 scuff-em calculation at a fixed Bloch wavevector. The scuff-em workflow offers two options for evaluating such cubatures:
- You can design and implement your own cubature scheme
involving your own custom-chosen weights and points
. In this case, you will use
--byOmegakBlochcommand-line option to instruct a scuff-em application code to report values of the quantity at each of your points (this output will typically be written to file with extension
.byXikBloch), then compute the weighted sums yourself in e.g. julia.
- Alternatively, you can ask scuff-em to perform the
BZ integration internally, using one of several
built-in cubature schemes. In this case the BZ-integrated
quantities will typically be written to
an output file with extension
.byXi. You will also get an output file named
.byXikBlochthat reports the Bloch-vector-resolved integrand samples chosen internally by the scuff-em BZ integrator.
Command-line options for customizing internal BZ integration
If you choose the second option above, you may specify various command-line options to customize the algorithm used by scuff-em 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).
Sets the order (accuracy parameter) of the integration
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.
This option lets you tell scuff-em that your integrand function is invariant under 2, 4, or 8-fold 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 command-line options above, this section describes the various integration algorithms available and how they are affected by the command-line parameters.
Integration methods for 1D Brillouin zones
For one-dimensional Brillouin zones, there is only one
BZ integration method available---namely,
--BZIMethod CC, the default)
with the number of integrand samples either fixed
or chosen adaptively until user-specified error tolerances
More specifically, you may say either
--BZIOrder [11 | 13 | 15 | ... | 97 | 99]
The former option selects fixed-order 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
using this algorithm.
In this case the number of sample points will be chosen
automatically subject to the values you select for the
Symmetry factors for 1D Brillouin zones
For 1D Brillouin zones, the only allowed
value of the
--BZSymmetryFactor option is
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
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.)
Clenshaw-Curtis cubature (
Nested 2D fixed-order or adaptive Clenshaw-Curtis cubature.
For fixed-order nested CC cubature with
NNsample points per dimension, say
--BZIOrder NN. Here
NNmust be an odd integer between 11 and 99 inclusive.
For adaptive 2D CC cubature (subject to your specified values of
Triangle cubature (
This algorithm divides the Brillouin zone into 8 triangles and applies a fixed-order triangle cubature scheme to each triangle, omitting repetition of triangles that are symmetry-equivalent given the value you specified for
--BZSymmetryFactor(see pictures below).
For this algorithm, the allowed values of
--BZIOrderare 1, 2, 4, 5, 7, 9, 13, 14, 16, 20, or 25.
Polar cubature (
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 Clenshaw-Curtis quadrature, and the quadrature is evaluated using rectangular-rule 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
--BZIOrderin the form where and are each odd integers between 11 and 99 inclusive.
Thus, for example,
--BZIOrder 3321specifies that the integral is to be evaluated via 33-point CC cubature, while the integral is to be evaluated via 21-point rectangular-rule cubature.
Polar cubature with change of variables (
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 variables
These modifications are useful for cases in which the free-space 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 solid-state physics, and changing variables to introduces a Jacobian factor that neutralizes these singularities to yield a better-behaved integrand.
Special values for rotationally-invariant integrands
As discussed above, for the
Polar2 integration methods
the value of the
--BZIOrder option is interpreted as the composite
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 rotationally-symmetric geometries in which the BZ integrand is independent of , you can specify to indicate that the integral is to be evaluated by a 1-point 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
take the value 2, 4, or 8, specifying that the Brillouin-zone
integrand obeys symmetries as follows:
We have so the BZ integration may be restricted to just the right half of the BZ.
We have so the BZ integration may be restricted to just the upper-right 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 Brillouin-zone integration with various values of the command-line 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