Mode Decomposition
Meep contains a feature to decompose arbitrary fields into a superposition of the harmonic modes of a given structure using the eigenmode solver MPB. This section provides a description of the analytical as well as implementation details of this feature. A tutorial example is provided in Tutorials/Mode Decomposition.
Theoretical Background
The analytical theory for waveguide mode decomposition is described in Chapter 31 ("Modal methods for Maxwell's equations") of Optical Waveguide Theory by Snyder and Love.
Consider a waveguide with propagation axis along the direction and constant cross section in the transverse direction . For a given angular frequency ω we can solve for the eigenmodes of the structure. Thus, arbitrary fields of the form and can be decomposed into a basis of these eigenmodes:
β are the propagation wavevectors and α are the basis coefficients. Mode decomposition involves solving for these unknown quantities. The following steps are involved in the computation:

In Meep, compute the Fouriertransformed fields and on a surface that is transverse to the waveguide and stored in a
dft_flux
object. 
In MPB, compute the eigenmodes and as well as the propagation wavevectors β for the same crosssectional structure.

Compute the coefficients α for any number of eigenmodes
This is all done automatically in Meep using the get_eigenmode_coefficients
routine.
Function Description
The modedecomposition feature is available via the meep::fields::get_eigenmode_coefficients
function callable from Python or C++. This function makes use of several lowerlevel functions which are described in more detail below. The C++ prototype for this routine is:
std::vector<cdouble>
fields::get_eigenmode_coefficients(dft_flux *flux,
direction d,
const volume &where,
std::vector<int> bands,
kpoint_func k_func=0,
void *user_data=0)
The following are the parameters:

flux
is adft_flux
object containing the frequencydomain fields on a crosssectional slice perpendicular to the waveguide axis 
d
is the direction of the waveguide axis 
where
is avolume
describing the crosssectional slice 
bands
is an array of integers corresponding to the mode indices 
user_func
is an optional function you supply to provide an initial guess of the wavevector of a mode with given frequency and band index having the following prototype:
vec (*kpoint_func)(void user_data, double freq, int band);
The return value of get_mode_coefficients
is an array of type std::complex<double>
(shortened to vector<cdouble>
) of length 2*num_freqs*num_bands
where num_freqs
is the number of frequencies stored in the flux
object (equivalent to flux>Nfreq
) and num_bands
is the length of the bands
input array. The expansion coefficients for the mode with frequency nf
and band index nb
are stored sequentially as α, α starting at slot 2*nb*num_freqs+nf
of this array:
std::vector<cdouble> coeffs=f.get_eigenmode_coefficient(...)
fields::get_eigenmode_coefficients(dft_flux *flux,
direction d,
const volume &where,
std::vector<int> bands,
kpoint_func k_func=0,
void *user_data=0);
int num_bands = bands.size();
int num_freqs = Flux>Nfreq;
for(int nb=0; nb<num_bands; nb++)
for(int nf=0; nf<num_freqs++; nf++)
{
// get coefficients of forward and backwardtraveling
// waves in eigenmode bands[nb] at frequency #nf
cdouble AlphaPlus = coeffs[2*nb*num_freqs + nf + 0];
cdouble AlphaMinus = coeffs[2*nb*num_freqs + nf + 1];
...
Normalization
The coefficients computed by get_eigenmode_coefficients
are normalized to ensure that their squared magnitude equals the
power carried by the corresponding eigenmode, i.e.
where is the power carried by traveling eigenmode . This is discussed in more detail below.
Related Computational Routines
Besides get_eigenmode_coefficients,
there are a few
computational routines in libmeep
that you may find useful
for problems like those considered above.
Computing MPB Eigenmodes
void *fields::get_eigenmode(double &omega,
direction d, const volume &where,
const volume &eig_vol,
int band_num,
const vec &kpoint, bool match_frequency,
int parity,
double resolution,
double eigensolver_tol);
Calls MPB to compute the band_num
th eigenmode at frequency omega
for the portion of your geometry lying in where
which is typically a crosssectional slice of a waveguide. kpoint
is an initial starting guess for what the propagation vector of the waveguide mode will be. This is implemented in mpb.cpp.
Working with MPB Eigenmodes
The return value of get_eigenmode
is an opaque pointer to a data structure storing information about the computed eigenmode, which may be passed to the following routines:
// get a single component of the eigenmode field at a given point in space
std::complex<double> eigenmode_amplitude(const vec &p, void *vedata, component c);
// get the group velocity of the eigenmode
double get_group_velocity(void *vedata);
// free all memory associated with the eigenmode
void destroy_eigenmode_data(void *vedata);
These are implemented in mpb.cpp.
Exporting FrequencyDomain Fields
void output_flux_fields(dft_flux *flux, const volume where,
const char *HDF5FileName);
void output_mode_fields(void *mode_data, dft_flux *flux,
const volume where,
const char *HDF5FileName);
output_flux_fields
exports the components of the frequencydomain fields stored in flux
to an HDF5 file with the given filename. where
is the volume
passed to the flux
constructor. In general, flux
will store data for fields at multiple frequencies.
output_mode_fields
is similar, but instead exports the components of the eigenmode described by mode_data
which should be the return value of a call to get_eigenmode
.
These are implemented in dft.cpp.
Computing Overlap Integrals
std::complex<double> get_mode_flux_overlap(void *mode_data,
dft_flux *flux,
int num_freq,
const volume where);
std::complex<double> get_mode_mode_overlap(void *mode1_data,
void *mode2_data,
dft_flux *flux,
const volume where);
get_mode_flux_overlap
computes the overlap integral between the eigenmode described by mode_data
and the fields stored in flux
for the num_freq
th stored frequency, where num_freq
ranges from 0 to flux>Nfreq1
. mode_data
should be the return value of a previous call to get_eigenmode.
get_mode_mode_overlap
is similar, but computes the overlap integral between two eigenmodes. mode1_data
and mode2_data
may be identical, in which case you get the inner product of the mode with itself. This should equal the group velocity of the mode based on the MPB's normalization convention.
These are implemented in dft.cpp.
How Mode Decomposition Works
The theoretical basis of the modedecomposition algorithm is the orthogonality relation satisfied by the normal modes:
where the inner product involves an integration over transverse coordinates:
where is any surface transverse to the direction of propagation and is the unit normal vector to (i.e. just in the case considered above). The normalization constant is a matter of convention, but in MPB it is taken to be the group velocity of the mode, , times the area of the crosssectional surface : .
Now consider a Meep calculation in which we have accumulated frequencydomain and fields on a dftflux
object located on a crosssectional surface . Invoking the eigenmode expansion and choosing the origin of the axis to be the position of the crosssectional plane, the tangential components of the frequencydomain Meep fields take the form:
We have used the wellknown relations between the tangential components of the forwardtraveling and backwardtraveling field modes:
Taking the inner product of both equations with the and fields of each eigenmode, we find
Thus, by evaluating the integrals on the LHS of these equations — numerically, using the MPBcomputed eigenmode fields and the Meepcomputed fields as tabulated on the computational grid — and combining the results appropriately, we can extract the coefficients . This calculation is carried out by the routine meep::fields::get_mode_flux_overlap
. Although simple in principle, the implementation is complicated by the fact that, in multiprocessor calculations, the Meep fields needed to evaluate the integrals are generally not all present on any one processor, but are instead distributed over multiple processors, requiring some interprocess communication to evaluate the full integral.
The Poynting flux carried by the Meep fields may be expressed in the form:
Thus, the power carried by a given forward or backwardtraveling eigenmode is given by:
or
$$ \textit{power} = \alpha_n^\pm^2
where we have defined .
These coefficients are the quantities reported by
get_eigenmode_coefficients.