There is a newer version of the record available.

Published April 17, 2023 | Version 1
Software Open


  • 1. Zuse Institute Berlin



The software RPExpand is designed to provide easy and efficient access to modal expansions of optical quantities based on Riesz projections [1,2,3,4]. Furthermore, modal expansions based on normalized eigenmodes and associated expansion coefficients are supported [5,6,7].


Get a copy of the two class folders and the directory containing post processes for JCMsuite scattering simulations

  • @Scattering (interface to JCMsuite)
  • @RieszProjection
  • postprocesses


and add their parent directory to your matlab path or copy them into the directory from which your scripts are executed.


Information about the usage of RPExpand can be obtained using the help function. E.g., typing help RieszProjection in the Matlab terminal will print a short description of the class and will list its properties and methods. Furthermore, differnent examples are provided that can be run to get familiar with the functionalities and can serve as a guideline for setting up new projects. Additionally, the following section summarizes the prerequisites needed in order to use this code with the FEM solver JCMsuite.

The class 'RieszProjection' can be used independently of JCMsuite. An example is given in the script 'quantumExample.m'. This example runs very fast, does not require an installation of JCMsuite and covers the most important features of the code. Therefore it is recommended to run this script to get a first impression of RPExpand.


RPExpand has been tested with MATLAB R2019b under Linux but should be independent of the platform. The class 'RieszProjection' can be used if an interface to a software capable of solving linear systems is available. This interface must have the form of a callable object that is described in more detail in later. RPExpand includes the ready-to-use interface to JCMsuite 'Scattering'. A current release of JCMsuite can be downloaded from the homepage of JCMwave. Please follow the installation and activation instructions provided there. A 14 day trial licence is available free of charge. It does not allow for computing different scattering problems in parallel.

As a prerequisite for the usage of RPExpand with JCMsuite, its Matlab interface must be set up and a daemon has to be started. A minimal example for a small workstation would be:

% Set up thirdparty support
jcm_root = '/path/to/local/installation/JCMsuite.x.x.x'

% shutdown a possibly running daemon

% register a new resource
options = struct(...
                 'Hostname', 'localhost', ...
                 'Multiplicity', 1, ...
                 'NThreads', 2, ...

If you have access to a remote machine via ssh, you can pass hostname, username, and other parameters. For further details on the usage of the daemon with a queue or a cluster, please refer to the daemon command reference. If you have access to a license which allows for the parallel submission of different jobs, you can increase the multiplicity to a value larger than one. The parameter 'NThreads' refers to the number of threads per job.

A scattering project for JCMsuite compatible with RPExpand must contain the following files:

The trailing 't' at the file extentions marks template files, which are used for parameter substitution. Whereas in the last two cases it is optional, 'project.jcmpt' and 'sources.jcmt' must contain predefined parameters. A minimal project file is:

Project {
  Electromagnetics {
    TimeHarmonic {
      Scattering {
        PML {
        FieldComponents = Electric
        Accuracy {
          FiniteElementDegree = %(finiteElementDegree)e

The parameters 'pml' and 'finiteElementDegree' will be set by RPExpand. The former as the perfectly matched layers (PML) should usually be the same for each integration point. The latter, as this way it can be conveniently set in the script and as some post processes require this information in order to determine the integration order. If there exists a file 'pml.log' in the project directory, the PML will be based on its content, otherwise, it will be created with the first scattering problem. As the scattering problem has to be solved for many different frequencies, the source file must contain the definition Omega = %(omega)e.

Custom parameters to be substituted in an input file can be passed to the FEM solver via the interface 'Scattering'. We refer to the examples, the parameter reference and the documentation of the matlab interface of JCMsuite.

Expansion of optical quantities

Currently the following quantities are available for expansion:

  • the total field
  • dipole emission and normalized decay rate of a point source
  • electromagnetic field energy flux
  • electric field energy
  • electromagnetic field absorption
  • power radiated to the far field
  • radiation pattern

If you want to expand a quantity which is not yet readily available, you can either contact us or add a new post process to the existing ones. Doing the latter, the following paragraph will give you some guidance.

For an efficient integration along the contours, not the total field is integrated but the target quantities individually, which often are scalars. N dimensional arrays are flattened to the size nx1 for integration and can be reshaped to their original sizes afterwards. This is done for the expansion of the radiation pattern. The name of the quantity must match the name of the jcmpt-file located in the directrory 'postprocesses'. It is defined in JCMsuite with the primitive 'OutputQuantity' or, in cases where the integrand is a python expression, 'IntegralName'. If a scalar is returned, the creation of the new jcmpt-file is all you need to do. If the result needs further post processing this can be done in the method '@Scattering/collect.m'. Please be aware that you can only expand holomorphic expressions. For the visualization of your result you can add a new private method which is called by '@RieszProjection/plot.m'.

Notes on a usage independent of JCMsuite

If you want to use your own software to solve the linear systems at the abszissae, you can still use the class 'RieszProjection'. The integration is based on quadrature rules of the form:

\(I(\omega_0) = \sum_i c_i(\omega_0) f(\omega_i)\).

Using the method 'getContours' you generate the weights \(c_i\) and the frequencies \(\omega_i\). All you need to provide are the function values \(f(\omega_i)\), e.g., solutions of the second order Maxwell's equation:

\(\nabla \times \mu^{-1}\nabla\times\bf{E}(\bf{r},\omega)-\omega^2\epsilon(\omega)\bf{E}(\bf{r},\omega)=i\omega\bf{J}(\bf{r}).\)

Usually, it is not necessary to integrate the total field. Instead, RPExpand is designed to evaluate the target quantities at the integration points and, subsequently, integrate them separately. As quantities based on quadratic forms require solutions of the two independent function values \(f(\omega_i)\) and \(f^*(\omega_i^*)\) the partial differential equation (the scattering problem) and the postprocess to extract the quantity of interest must be done in two independent steps.

The interface to a custom solver must be defined as a callable object and will be used as follows:

The frequencies which discretize the contours are saved in a cell array 'contours' of size (1, n) where n is the number of contours. Each element is a complex double array of shape (k,m), which corresponds to a contour. Here, k is the number of nodes per subinterval. Unless using higher order quadrature methods, there is only one interval, i.e., m = 1 and k is the total number of integration points. In a first step, the cell array 'contours' is passed to your function, which is expected to return a cell array of size (1, n) or (n, 1) whose elements are objects of size (1, k*m). For instance, a cell array containing numeric arrays, each of them representing the solution of a partial differential equation. In a second step, each element of the returned cell vector is passed to your custom function again with some additional arguments: v = f(sc_results,quantity,keys). The additional arguments 'quantity' and 'keys' define the quantity which has to be extracted from the results returned previously and additional meta data (e.g., about the shape), respectively. The quantity is defined as a 'char' vector (e.g., 'ElectricFieldEnergy') and the last argument is a scalar struct. The returned value 'v' is expected to be a cell containing numeric arrays of size (d, k*m), where d is a custom number, which in many cases will be one, but can be larger as the examples for the expansion of the radiation pattern shows. In the case of a quantity quadratic in the solutions of the scattering problems, the first input is of size (n, 2) and the second column contains the solutions of the conjugated scattering problems. Otherwise it is of shape (n, 1). An example for a function meeting the described criteria is given in the file 'quantumExample.m' which implements the quantum transition problem in 1D.

In a last step, you have to edit the constant properties 'linearQ' and 'quadraticQ'. Insert all quantities available for expansion, depending on whether they are linear or quadratic in \(f(\omega)\).


[1] L. Zschiedrich, et al., Phys. Rev. A 98, 043806 (2018)

[2] F. Binkowski, et al., Phys. Rev. B 100, 155406 (2019)

[3] F. Binkowski, et al., Phys. Rev. B 102, 035432 (2020)

[4] F. Binkowski, et al., J. Comput. Phys. 419, 109678 (2020)

[5] P. Lalanne, et al., Laser Photonics Rev. 12, 1700113 (2018)

[6] T. Wu, et al., J. Opt. Soc. Am. A 38, 1224 (2021)

[7] F. Betz, et al., Phys. Status Solidi A 220, 2200892 (2023)


We acknowledge funding by the German Federal Ministry of Education and Research (BMBF Forschungscampus MODAL, project 05M20ZBM) and by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under Germany's Excellence Strategy - The Berlin Mathematics Research Center MATH+ (EXC-2046/1, project ID: 390685689). This project has received funding from the EMPIR programme co-financed by the Participating States and from the European Union's Horizon 2020 research and innovation programme (project 20FUN02 POLIGHT).


Files (128.7 kB)

Name Size Download all
128.7 kB Preview Download

Additional details