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 the --byOmegakBloch command-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 .byOmegakBloch or .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 .byOmega or .byXi. You will also get an output file named .byOmegakBloch or .byXikBloch that 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).

--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 [2|4|8]

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, Clenshaw-Curtis quadrature (--BZIMethod CC, the default) with the number of integrand samples either fixed or chosen adaptively until user-specified error tolerances are achieved.

More specifically, you may say either

  • --BZIOrder [11 | 13 | 15 | ... | 97 | 99]


  • --BZIOrder 0

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 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 command-line 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.)

  • Clenshaw-Curtis cubature (--BZIMethod CC)

    Nested 2D fixed-order or adaptive Clenshaw-Curtis cubature.

    For fixed-order nested CC cubature with NN sample points per dimension, say --BZIOrder NN. Here NN must be an odd integer between 11 and 99 inclusive.

    For adaptive 2D CC cubature (subject to your specified values of --BZIRelTol, --BZIAbsTol, and --BZIMaxEvals) say --BZIOrder 0.


  • Triangle cubature (--BZIMethod TC)

    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 --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 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 --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 33-point CC cubature, while the integral is to be evaluated via 21-point rectangular-rule 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 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 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 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 --BZSymmetryFactor may take the value 2, 4, or 8, specifying that the Brillouin-zone 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 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