Dyadic Green's Functions and LDOS above a material halfspace
In this example, we exploit scuffem's support for 2D periodic geometries by using scuffldos to compute dyadic Green's functions (DGFs) and electromagnetic local densities of states (LDOS) at evaluation points lying above an infinite planar dielectric interface.
Because this geometry happens to be amenable to an analytical treatment, we will also check the numerical results of scuffem against the results of analytical calculations of the LDOS.
The files for this example may be found in the
share/scuffem/examples/HalfSpaceLDOS
subdirectory
of your scuffem installation.
Review: What we are computing and why we need Brillouinzone integration
For the benefit of users who may be wondering why numerical Brillouinzone integration is needed for this calculation, this section briefly reviews some of the theoretical background. (This material is discussed in more detail in the memo Computation of Green's Functions and LDOS in scuffem). If you just want to get started running the calculations in scuffem, feel free to skip to the next section.
Green's functions and LDOS
The scattering dyadic Green's function of a material geometry is a matrix whose entry is proportional to the component of the scattered electric field at due to a directed point electric dipole source at , all quantities having time dependence . (We could also consider the magnetic Green's function , which describes the scattered magnetic fields arising from a point magnetic dipole excitation.) The subscripts "D" and "S" label the "destination" and "source" points; if these coincide , we are talking about onepoint DGFs; otherwise we have twopoint DGFs.
The onepoint scattering DGFs for a material geometry contain information on how the bodies in the geometry modify the electromagnetic local density of states (LDOS) in the vicinity of the body. More specifically, the full LDOS at a point may be decomposed as a sum of electric and magnetic contributions, which themselves may be further decomposed into vacuum and scattering contributions:
Computing Green's functions
The electric and magnetic DGFs for a given source point (and any number of destination points) collectively describe the results of 6 scattering problems (with incident fields produced by electric and magnetic dipoles oriented in each of the 3 cartesian directions), andfor compact geometriescould be computed in scuffem simply by assembling 6 righthandside vectors, solving the BEM system for the corresponding surface currents, and computing the scattered fields at all desired destination points.
However, the situation is more complicated for infinitely extended geometries. Recall that such geometries are described in scuffem by a unitcell geometry together with a basis of vectors for a (one or twodimensional) lattice . By restricting our attention to scattering problems in which the incident fields (and thus all fields and currents) satisfies the Blochperiodicity condition for some Bloch wavevector , we can get away with solving only for the surface currents within the unit cell, the currents everywhere else being then fixed by the Bloch periodicity. This is the basic magic trick that allows infinite geometries to be modeled on finite computers.
Of course, the restriction to scattering problems involving only Blochperiodic incident fields is in fact highly restrictive; indeed, it excludes the particular problem we must solve to compute DGFs and LDOS, in which the incident field emanates from a point source at and is certainly not Blochperiodic. But here we can avail to the linearity of Maxwell's equations and the magic of Fourier analysis. Instead of solving a scattering problem involving the field of a point source at , we consider the field of an infinite phased array of point sources, located at points with phases , for all . This is a Blochperiodic field, so we can take it as the incident field in a scuffem scattering problem, and the resulting scattered fields at define the Blochperiodic scattering DGFs To recover the nonperiodic DGF from , we now imagine solving and summing the results of many scattering problems involving our infinite array of point sources phased with many different Bloch vectors , chosen cleverly to ensure that the contributions of the single point source at add constructively, while the contributions of all other point sources in the array add destructively and cancel. The formal expression of this procedure is a sort of inverse Fourier transform relating to the average of over the Brillouin zone (BZ) of the reciprocal lattice: or, evaluating the integral approximately by numerical cubature, where are the points and weights in a cubature rule for integration over the Brillouin zone. Thus, getting the LDOS at a single point in a periodic geometry requires summing the results of many scattering calculations at various Bloch vectors lying in the Brillouin zone of the reciprocal lattice.
The question of how to choose the cubature rulethat is, the set of Bloch vectors at which the integrand is sampled and the weights with which the samples are combinedto yield accurate integral estimates without exorbitant computational cost is a tricky one, and one whose answer varies from problem to problem. There are two ways to proceed:
 You can design your own cubature scheme for integration over
the Brillouin zone, and simply ask scuffldos to give you
values of the integrand
at specific Bloch vectors .
In this case you will pass the
OmegakBlochFile
commandline argument to scuffldos, and you will get back a data file (with extension.byOmegakBloch
) reporting values of at each point you specified. (See below for examples).
 Alternatively, you can use the builtin Brillouinzone integrator
in scuffldos to perform the BZ integral automatically. In this
case you will get back a data file with extension
.LDOS
reporting values of and at each frequency and evaluation point you specified, plus a data file with extension.byOmegakBloch
reporting the Bloch vectors at which the internal integrator sampled the integrand and the values of at each sample point. In this case you will probably want to set various commandline options to optimize the integral evaluation; this is illustrated below and discussed in more detail on the page Brillouinzone integration in scuffem.
gmsh geometry files for unitcell geometries
To compute DGFs and LDOS in a halfspace geometry, we will describe the interface between the upper (vacuum) and lower (dielectric) regions of our geometry as the infinite periodic replication of a unitcell mesh over the sites of a 2D square lattice. Because this particular geometry has infinite translational symmetry, we are free to choose the lattice constant however we like; here we will consider the two values nm and nm.
The gmsh geometry file Square_L_N.geo
describes a square with useradjustable parameters
L
and N
that may be configured using the setnumber
option on the gmsh command line to define the square
side length and the number of triangle edges per side length
(the meshing fineness). To produce unitcell meshes for
the two desired lattice constants, I go like this:
% gmsh 2 setnumber L 0.1 setnumber N 4 Square_L_N.geo o Square_L0P1.msh
% RenameMesh Square_L0P1.msh
% gmsh 2 setnumber L 0.01 setnumber N 4 Square_L_N.geo o Square_L0P01.msh
% RenameMesh Square_L0P01.msh
(Here RenameMesh
is a simple bash
script
that uses scuffanalyze to count the number of interior
edges in a surface mesh and rename the mesh file accordingly.)
These commands produce the files
Square_L0P1_40.msh
and Square_L0P01_40.msh.
I also repeat the
process with setnumber N 8
to yield finer meshes (with 176 instead
of 40 interior edges).
These meshes may be visualized in gmsh by saying e.g.
% gmsh Square_L0P1.40.msh
% gmsh Square_L0P01.176.msh
Note the following:

For 2D periodic geometries in scuffem, the lattice vectors must lie in the plane.

For surfaces that straddle the unitcell boundaries (as is the case here), each triangle edge that lies on any edge of the unit cell must have an identical image edge on the opposite side of the unit cell. An easy way to achieve this is to use extrusions in gmsh, as in the
.geo
file above. 
Among the various possible ways to discretize a square intro triangles, I have chosen one that preserves as much symmetry as possiblespecifically, the mesh is symmetric under the transformations , , and . This is important for exploiting symmetries to reduce the cost of Brillouinzone integration (see below).
The .msh
files needed to run the calculations in this
tutorial may be found in the mshFiles/
subfolder
of the HalfSpaceLDOS
example distributed with
scuffem.
scuffem geometry files
For each unitcell mesh I create scuffem geometry files describing three material geometries:

an infinitearea perfectly electrically conducting (PEC) plate at height ;

an infinite lossless dielectric () filling space for , with vacuum for ;

same as previous, but now with lossy metal (aluminum) instead of lossless dielectric.
The resulting .scuffgeo
files may be found in the
scuffgeoFiles/
subfolder of the HalfSpaceLDOS
example directory.
For example, the PECplate geometry at the coarser resolution
with the longer lattice constant is
PECPlate_L0P1_40.scuffgeo
,
while the aluminum halfspace geometry at the finer resolution
with the shorter lattice constant is
AlHalfSpace_L0P01_176.scuffgeo
.
The latter file looks like this:
# this comes from Phys Rev B **68** 245405
MATERIAL ALUMINUM
wp = 1.747e16;
gamma = 7.596e13;
Eps(w) = 1  wp^2 / (w * (w + i*gamma));
ENDMATERIAL
LATTICE
VECTOR 0.01 0.00
VECTOR 0.00 0.01
ENDLATTICE
REGION UpperHalfSpace MATERIAL Vacuum
REGION Exterior MATERIAL Aluminum
SURFACE Plate
MESHFILE Square_L0P01_176.msh
REGIONS Exterior UpperHalfSpace
ENDSURFACE
List of evaluation points
We'll compute onepoint DGFs and LDOS at four points,
located at distances of 1, 10, 100, and 1000 nm above
the surface on the axis (). Thus our evaluationpoint
file (EPFile
) looks like this:
0.0 0.0 0.001
0.0 0.0 0.010
0.0 0.0 0.100
0.0 0.0 1.000
Studying the integrand to select the optimal cubature rule
Ultimately we will want to launch a full Brillouinzoneintegrated calculation to get the total LDOS at our desired evaluation points and frequencies. However, before doing this we should acquire some insight into the behavior of the integrand so that we can make intelligent choices for the various options controlling Brillouinzone integration.
Checking Brillouinzone symmetry
A first thing to check is the symmetry of the integrand within the Brillouin zone (BZ). As discussed here, the builtin cubature routines for 2D BZ integration in scuffem can exploit 2, 4, or 8fold rotational symmetry. In this case, because

(1) the underlying (realspace) geometry is symmetric under , , and ,

(2) we have chosen a mesh discretization that preserves these symmetries, and

(3) we are computing a onepoint DGF,
we expect the BZ integrand to have 8fold rotational symmetry, and thus we should be able to restrict the domain of the BZ integration to just the irreducible Brillouin zone (a triangular wedge having 1/8 the total area of the BZ, as shown here) and multiply by 8 to get the full BZ integral.
To test this, we'll run a quick first calculation in which we compute just the integrand at 8 points related by symmetry. Specifically, for a square lattice with lattice constant (reciprocal lattice constant ) we'll consider the point and the 7 other points obtained by flipping signs and interchanging components. Here's a picture showing the locations of these 8 points within in the BZ:
Working at angular frequency rad/sec,
we create a little text file named EightPoints
listing the
values at which we want to sample the BZ integrand:
10 +25 +17
10 +17 +25
10 17 +25
10 25 +17
10 25 17
10 17 25
10 +17 25
10 +25 17
Now run scuffldos:
% scuffldos geometry scuffgeoFiles/AlHalfSpace_L0P1_40.scuffgeo OmegakBlochFile EightPoints EPFile EPFile
This produces the file AlHalfSpace_L0P1_40.byOmegakBloch
,
in which each individual line reports values
for a single evaluation point, frequency, and Bloch vector.
Looking at just the portion of the file corresponding to the
evaluation point at nm, we see this:
This shows that, as expected, the Brillouinzone integrands for
the electric and magnetic LDOS are
invariantup to small numerical errors on the order
of around 1%, which are reduced to around 0.1% when we rerun
the calculation on the finer [N=176] meshunder the 8fold
symmetry relating the various Brillouinzone points here.
Later we will exploit this symmetry by specifying the commandline
option BZSymmetryFactor 8
to reduce the cost of
evaluating the full BZ integral.
Incidentally, notice that the individual components of the DGFs (such as ) do not obey the full set of symmetries obeyed by the LDOS. You will need to take this into account if your goal is to compute DGFs and not just LDOS.
Observing the variation of the integrand with
Another diagnostic test is to inspect the variation of the integrand with the magnitude of the Bloch vector. For this purpose we'll compute the integrand at points lying on the axis between the origin and the edge of the BZ, again working with the nm lattice at . We create a file called XAxisPoints
specifying the () points we want to sample:
10 00.0 0.0
10 1.57 0.0
10 3.14 0.0
...
10 29.8 0.0
10 31.4 0.0
Run scuffldos:
% scuffldos geometry scuffgeoFiles/AlHalfSpace_L0P1_40.scuffgeo OmegakBlochFile XAxisPoints EPFile EPFile
This produces the file AlHalfSpace_L0P1_40.byOmegakBloch.
We plot the BZ integrand for the electric LDOS
versus at nm:
% gnuplot
gnuplot> FILE='AlHalfSpace_L0P1_40.byOmegakBloch'
gnuplot> set yrange [100:100]
gnuplot> plot FILE u ($3==0.1 ? $6 : 1/0):8 w lp pt 7 ps 1
Whoops! There seems to be some sort of blowup happening near , i.e. near the point where the Bloch wavelength coincides with the freespace photon wavelength. Zooming in on the region near , the data look like this (where we have additionally plotted data for evaluation points at a couple of other heights off the surface):
The important things to read off from this plot are:

As from below, the integrand blows up like

On the other side of , the integrand decays exponentially at a dependent rate.
This behavior is expected for cases in which the photon wavenumber lies within the Brillouin zone; the squareroot blowup in the LDOS is just the usual van Hove singularity, while the exponential decay for larger arises because plane waves with transverse wavelength shorter than the freespace photon wavelength are evanescent (not propagating) waves that decay exponentially with distance from the surface.
(BTW, here is the gnuplot script that produces the two plots above.)
Running full Brillouinzoneintegrated calculations
Now it's time to launch runs for full BZintegrated LDOS. For
this purpose we will want to set appropriate values for the
commandline options controlling Brillouinzone integration.
We have already determined that the integrand has
8fold symmetry, so we will be specifying BZSymmetryFactor 8.
The next task is to choose a cubature
rule ("BZ integration method"), i.e. a strategy for sampling
the integrand at appropriatelychosen
points to estimate the integral accurately with as few
samples as possible.
The presence of the singularity
in the Brillouin zone makes this a good candidate for
the Polar2
strategy, which performs the BZ integral
using a polar decomposition
,
with the integrals over the regions
and handled separately, and with the
change of variables
to neutralize the squareroot singularity at
In this case the integrals
over and are performed using
and evaluation points,
respectively, where and
are combined into a single value specified as
the BZIOrder
commandline option, i.e.
For example, to use 31 points for the radial integral
and 11 points for the integral, we would
say BZIOrder 3111
.
(Actually, in this particular case, we expect to have full rotational symmetry, i.e. the integrand should be independent of , and thus we can specify the special value to request that the integral be evaluated by a onepoint cubature sampling just the point .)
Choosing an 11point integration rule for the radial integral and the abovediscussed 1point cubature for the integral, our scuffldos commandline looks like this:
% scuffldos geometry scuffgeoFiles/AlHalfSpace_L0P1_40.scuffgeo OmegaFile OmegaFile EPFile EPFile BZIMethod Polar2 BZIOrder 1102
As a comparison, I will also run calculations with
. This selects an adaptive integration
algorithm for the radial integral, which determines
the number and location of the sample points
automatically to achieve userspecified relative
and absolute error tolerances, which are set with the
BZIRelTol
and BZIAbsTol
commandline options:
% scuffldos geometry scuffgeoFiles/AlHalfSpace_L0P1_40.scuffgeo OmegaFile OmegaFile EPFile EPFile BZIMethod Polar2 BZIOrder 0002 BZIRelTol 1.0e2 BZIRelTol 1.0e8
I also run calculations on the finerresolution mesh (the _176.scuffgeo
files) and the
mesh with lattice constant (the L0P01
files).
Here's a plot of the LDOS enhancement at :
The green dots are for the calculation with all numerical resolution parameters set to their maximum (most accurate and most expensive) values: finest mesh, smallest lattice constant, and adaptive BZ cubature. These points agree well with the analytical prediction up to the highest frequencies, where they start to deviate.
The other data symbols are the results of calculations in which one fineness parameter is relaxed to a lessexpensive value:
 The blue triangles are for the geometry with lattice constant nm (but still using a fine mesh with adaptive cubature). Evidently this degrades the quality of the numerical data completely. That makes sense, since the evaluation point is at a height nm above the surface. For larger values the data values for the two lattice constants are in better agreement.
 The orange triangles are for the geometry with lattice constant nm, but now meshed at coarser resolution (resulting in triangle edges instead of for the finer case). Evidently this doesn't degrade the quality of the data at all, so for this case we can run calculations with the (significantly cheaper) coarser mesh at no cost in accuracy.
 The yellow stars are for the nm geometry with the fine () mesh resolution, but with the BZ integration performed using a fixedorder numerical integration ( samples) instead of adaptive quadrature for the radial part of the BZ integral. Evidently this yields data that agree with the more accurate adaptivequadrature data at most frequencies, but deviate over a range of frequencies in the range and again at high frequencies.