Multem 3: An updated and revised version of the program for transmission and band calculations of photonic crystals

We present here Multem 3, an updated and revised version of Multem 2, which syntax has been upgraded to Fortran 2018, with the source code being divided into modules. Multem 3 is equipped with LAPACK, the state-of-the art Faddeeva complex error function routine, and the Bessel function package AMOS. The amendments significantly improve both the speed, convergence, and precision of Multem 2. Increased stability allows to freely increase the cut-off value LMAX on the number of spherical vector wave functions and the cut-off value RMAX controlling the maximal length of reciprocal vectors taken into consideration. An immediate bonus is that Multem 3 can be reliably used to describe bound states in the continuum (BICs). To ensure convergence of the layer coupling scheme, it appears that appreciably larger values of convergence paramaters LMAX and RMAX are required than those reported in numerous published work in the past using Multem 2. We hope that Multem 3 will become a reliable and fast alternative to generic commercial software, such as COMSOL Multiphysics, CST Microwave Studio, or Ansys HFSS, and that it will become the code of choice for various optimization tasks for a large number of research groups. The improvements concern the core part of Multem 2, which is common to the extensions of Multem 2 for acoustic and elastic multiple scattering and to the original layer-Kohn-Korringa-Rostocker (LKKR) code. Therefore, the enhancements presented here can be readily applied to the above codes as well. New version program summary Program Title: Multem 3CPC Library link to program files: https://doi.org/10.17632/3hnz5wdvzr.1Developer's repository link: Multem-3Licensing provisions: MITProgramming language: Fortran 2018Journal reference of previous version: N. Stefanou, V. Yannopapas, A. Modinos Computer Physics Communications 132 (2000) 189Does the new version supersede the previous version?: Yes.Reasons for the new version: The linear algebra part and the standard of the programming language are outdated. More precise algorithms for calculating special functions (Faddeeva complex error function and Bessel functions) have been introduced.Summary of revisions:(i)Code refactoring. Language standard was updated to Fortran 2018. The source code was divided into modules with CMake building system as highlighted in the flowchart of Fig. 1.(ii)Transition to LAPACK. Original routines ZGE and ZSU were changed to the LAPACK routines zgetrf and zgetrs, resulting in a significant reduction of more than three orders of magnitude in the numerical errors that occur when solving systems of linear equations.(iii)Augmenting the implementation of the Faddeeva complex error function. Pendry's routine CERF [1] calculating the Faddeeva complex error function was substituted by the C++ SciPy code by Johnson [2,3] The latter shows significantly improved results, with the relative error limited to the range between 1.69×10−14 and 2.67×10−13 for the real and imaginary parts of the function values, respectively, for the arguments from 10−8 to 108.(iv)Using the Bessel function package AMOS. Original routines for Bessel functions calculation, used to calculate the T-matrix of a single sphere in routine TMTRX, were substituted by the routines originally written by Donald E. Amos [4], which can be found at [5]. Interfaces [6] were used to link the code with the original Amos routines. The benefit of this change was not directly evaluated. Nevertheless, the Amos library provides in our experience the most accurate results among all the libraries used to compute Bessel functions.(v)Increased stability allowing to freely increase the cut-off value LMAX on the number of spherical vector wave functions. LMAX is one of the most important convergence parameters of Multem. Nevertheless, any change of LMAX was very inconvenient as it required changing its value in almost any original routine. Furthermore, an upper cut-off value LMAXD was imposed on allowable LMAX values. Any attempt to increase LMAX above a preset value of LMAXD=7 required changing LMAX and LMAXD, together with additional dependent parameters NELMD and NDEND (cf. Table 1 of supplementary material). Contrary to that, Multem 3 allows calculations to be performed with any value of LMAX set at a single code location.(vi)Increased stability allowing to freely increase the cut-off value RMAX controlling the maximal length of reciprocal vectors taken into consideration. RMAX is the second of the main convergence parameters of Multem. Similarly to LMAX, varying RMAX required adjusting parameters in almost all original routines. Furthermore, an upper cut-off value RMAXD was imposed on allowable RMAX values. Any attempt to increase RMAX above a preset value of RMAXD required changing both RMAX and RMAXD, together with additional dependent parameter IGD (cf. Table 2 below). Contrary to Multem 2, Multem 3 allows calculations to be performed with much larger values of RMAX which is now set at a single code location.(vii)Selective multipole expansion. Multem 3 has been extended by a selective multipole expansion option in order to be able to isolate the contribution of any given multipole. This allows one to perform a selective multipole analysis of a resonant 2D structure. We have adapted the routines SETUP and PCSLAB, which are used to construct the coefficient matrices for linear systems of equations and to calculate the transmission/reflection matrices for a plane of spheres, accordingly.Nature of problem: Calculation of the transmission, reflection and absorption coefficients, and complex band structure of a photonic crystal slab.Solution method: Solution of Maxwell's equations in combination with ab-initio multiple-scattering theory for 2D periodic systems and a layer-by-layer coupling scheme.Additional comments including restrictions and unusual features: As its predecessor, Multem 3 can work only with layers of non-overlapping spheres arranged into simple Bravais lattice. By default, the underlying Bravais lattice has to be same in each plane of spheres (i.e. it is not possible to stack a triangular lattice on top of a square lattice, but it is possible to arrange different spheres (e.g. of different radius or of material) on the same underlying lattice in different planes). Additionally, an interface between planes may not cut any of the spheres. References [1]J. Pendry, Low Energy Electron Diffraction: The Theory and its Application to Determination of Surface Structure, Techniques of Physics, Academic Press, 1974.[2]Faddeeva Function — SciPy v1.11.3 Manual, [Online; accessed 26. Oct. 2023]. URL https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.wofz.html[3]Faddeeva Package - ab-initio, [Online; accessed 28. Oct. 2022]. URL http://ab-initio.mit.edu/wiki/index.php/Faddeeva_Package[4]D. E. Amos, Algorithm 644: A portable package for bessel functions of a complex argument and nonnegative order 12 (3) (1986). doi: https://doi.org/10.1145/7921.214331.[5]A Portable Package for Bessel Functions of a Complex Argument and Nonnegative Order, [Online; accessed 23. Oct. 2023]. URL https://netlib.org/amos[6]Ondřej Čertík, Fortran Utilities, [Online; accessed 14. Jun. 2023] (Oct. 2019). URL https://github.com/certik/fortran-utils