1 Introduction
Although the existence of a globally optimal multilayer thinfilm was proven [20, 21, 6], the optimization of multilayer thinfilms concerning reflectivity and transmittivity over wavelength and incidence angle remains a challenging task for the scientific community [2, 3, 12]. As such it gained attention in many publications [5, 17, 24, 7, 13] including various methods based on reinforcement learning [10, 22] and machine learning [9, 18]. Naturally, the need for a standardized, encompassing environment arises to make future research more trustful, comparable, and consistent. Here, we propose a comprehensive Python package that provides functionality to researchers to design and optimize multilayer thinfilms. At its core, the proposed TMMFast package (Transfer Matrix Method  Fast) is a parallelized and reimplemented revision of the existing TMM code that was initially published by S. Byrnes [4]. It implements Abèles transfer matrix method [1]
in Python to calculate transmission and reflection of an incident plane wave through a slab of layered thinfilm materials with different thicknesses. Given a broad, discretized spectrum of light irradiating a thinfilm under particular incident angles, the transmission and reflection must be calculated for each contributing wavelength at each angle of incident. The coefficients of transmission and reflection for a specific wavelength and angle of incidence depend on the thicknesses and dispersive and dissipative refractive index of the layers. Those coefficients can deviate significantly at wavelengths which differ only slightly for the same multilayer thinfilm. An intuitive approach to increase the computational speed of the response of a multilayer thinfilm is to parallelize the pairwise independent computations regarding wavelengths and angles. The parallelization is implemented via matrix operations based on the opensource numeric algebra package Numpy
[8] and the threadmanagement package Dask. The parallelization showed to reduce the computational time by with an additional factor on the order of available CPU cores by using Dask. Additionally, by using PyTorch Autograd [16], analytical gradients of the multilayer thinfilm layer thicknesses and refractive indices can be computed for local optimization. Since the implementation is done via PyTorch methods, it can be integrated into advanced Neural Networks. Moreover, the parallelized TMMFast code is used to implement an OpenAI Gym environment
^{1}^{1}1https://github.com/MLResearchAtOSRAM/gymmultilayerthinfilmthat can be used by physicists and reinforcement learning researchers. In the environment, the optimization of multilayer thinfilms is introduced as a sequence generation process and thus considered as a parameterized Markov decision process
[14], shown in section 6. An overview of the scope of the package is shown in Figure 1. All code is available on GitHub^{2}^{2}2https://github.com/MLResearchAtOSRAM/tmm_fast under MIT license.2 Physical Background of the Transfer Matrix Method (TMM)
Light of a particular wavelength passing from one into another material experiences a sudden change of the refractive index to that results in reflection and transmission of the incoming wave. A light wave with a wavevector k which is not parallel to the surface normal experiences refraction. The ratio of angle of incidence and angle of refraction is given by Snell’s law . Given the Fresnel equations, reflection and transmission coefficients and can be computed, where a distinction is made between s and p polarization by using a subscript
(1)  
(2) 
and being the respective magnetic permeability. Note that for = 0, i.e. vertical incidence, and . Based on the reflection and transmission coefficients, the reflectivity and transmittivity R and T can be computed as shown in Equation 3.
(3) 
where the subscript indicates the polarization. Note that, for absorptionless materials holds.
Now, consider a multilayer thinfilm with layers and the individual layers are denoted by . The light enters from an injection layer of semiinfinite thickness with a relative amplitude of 1 and exits the multilayer thinfilm in the outcoupling layer of semiinfinite thickness with . From the outcoupling layer, no light enters the thinfilm. The transmitted part of the light in layer , that travels in "forward" direction, ie. towards the layer with is denoted by and the reflected part that travels in "backward" direction is given by . By using the reflection and transmission coefficients r and t, the response of any layer can be written as
(4) 
where is the accumulated phase of the light wave when travelling through a layer with a specific thickness
and with wave vector
.Eventually, the total characteristic matrix of the multilayer thinfilm is given by
(5) 
Finally, to compute the reflection and transmission of the entire multilayer thinfilm, one needs to evaluate
(6) 
The transmission and reflection coefficients are separated
(7) 
and allow to compute the reflectivity and transmittivity via Equation 3. Note that one can easily calculate the partial transmission and reflection coefficients by only multiplying Equation 5 up to . The initial angle of incidence in the first layer is given by .
3 Implementation in Numpy
The key contribution of TMMFast is the parallelized handling of the characteristic matrix that reduces computational time. The matrix consists of three separate matrices, matrix which encompasses the accumulated phase, and the two matrices holding the coefficients of reflection and transmission, respectively. They are of shape , where and are the number of wavelengths and incident angles, respectively. To get the characteristic matrix , Numpy’s einsum^{3}^{3}3Documentation: https://numpy.org/doc/stable/reference/generated/numpy.einsum.html method allows to specify multiplication and contractions of different dimensions easily:
Finally is computed by multiplying out the thickness dimensions.
The entire function is called coh_tmm_fast or by coh_tmm_fast_disp. Both methods differ since the former assumes dispersionless materials whereas the latter accepts dispersive materials. An Example is give in Appendix A.1.
4 Core functionality speedup and dataset generation
To verify the speedup that the TMMFast package provides compared to the native implementation of Byrnes, ten multilayer thinfilms are generated with 21 layers are generated and evaluated in a spectral range from 400 nm to 700 nm at 100 equally spaced points. The angles of incidence range from 0° to 90° based on 20 equally spaced supporting points. The original TMM method requires a computation time of s for the evaluation while our proposed method computes the coefficients of the multilayer thinfilm in s which corresponds to an acceleration of 100.
Finally, the Python package Dask^{4}^{4}4https://dask.org/ manages parallel threads of the CPU to distribute the computation of different independent computations on all available CPU cores. The application is straightforward: by calling the coh_tmm_fast function implicitly inside the delayed() function of Dask, a list of all necessary computations is created. Then, the entire list is executed implicitly to create a computational graph for Dask, which is required to orchestrate the parallel threads efficiently. Lastly, the compute() method triggers the actual computation and the result is returned. By running coh_tmm_fast on all available threads, an additional speedup of the dataset creation on the order of the number of available CPU cores is possible. However, the benefit might decrease for very large computational clusters with many cores since the management of the parallel threads creates computational overhang. By calling multithread_coh_tmm from the TMMFast package, the computation is easily started. A Dataset with one million samples of 9layer multilayer thinfilms at 100 wavelengths and 10 angles of incidence can be created in approximately 40 min on an eightcore machine. Computing more layers, wavelengths points, and angles can increase the computational time significantly. An example is shown in the Appendix subsection A.2.
5 TMMTorch: Multilayer thinfilm gradients via Autograd
For optimization, computing gradients enable gradientbased optimization algorithms such as gradient descent [11]. These gradient based optimization methods generally converge to local minima with fewer iterations than other nongradient based algorithms such as the NelderMead downhill simplex algorithm [15]. The Python package PyTorch^{5}^{5}5https://PyTorch.org/ implements matrix multiplication methods which allow parameters to be automatically differentiated (Autograd
) via the chain rule
[16]. Autograd enables the user to compute the gradients of the input parameters without the necessity of deriving the gradient analytically and can be dynamically adapted to the problem. By using the TMMTorch subpackage of TMMFast, the gradients for a multilayer thinfilm can be readily computed. An example is shown in the Appendix A.3. The application of the transfer matrix takes slightly longer by using the PyTorch routines. Therefore, using the regular TMMFast algorithms is still advised to reduce computational time if gradients are not necessary.6 Environment for reinforcement learning
Reinforcement learning [19] is an area of machine learning concerned with how intelligent agents ought to take actions in an environment in order to maximize a notion of reward. The proposed code implements such an environment, where agents can stack, characterize and optimize multilayer thinfilms. Therefore, the generation of multilayer thinfilms is considered as a parameterized Markov decision processes [14] and thereby implemented as a sequence generation process: Beginning from , an agent subsequently executes parameterized actions that specify the thickness and material index of the th layer. These actions determine which material of which thickness to stack next, thereby consecutively forming a multilayer thinfilm as illustrated in Figure 1. The stacked layers and  optional  the optical characteristics of these intermediate multilayer thin films are provided to the agent as the environmental state . Given this state, the agent takes the next action until the predefined maximum number of layers is reached or the agent decides to terminate stacking. The optical characteristic of the final proposed multilayer thinfilm of layers, e.g. regarding reflectivity over wavelength and angle of incidence, is computed via the proposed transfermatrix method (TMMFast). Here, refers to the matrix of (dispersive and dissipative) refractive indices of the material of each layer. Each material is identified by the material identifier index . denotes the vector of layer thicknesses. When the agent starts to stack a thinfilm, and are initialized with zeros and get filled according to the taken actions with the layer thicknesses and (dispersive) refractive indices, respectively. The observed reflectivity is compared to a userdefined, desired reflectivity , in order to derive a notion of numeric reward, e.g. an inverted reconstruction error
Based on this reward, the agent learns — for example, based on Qlearning [10, 23] — to distinguish between good and bad actions and thus derive an optimal thinfilm design.
Whereas the contained physical methods are wellstudied and known for decades, the contribution of the code available at https://github.com/MLResearchAtOSRAM/gymmultilayerthinfilm lies in the implementation of an OpenAI gymrelated environment. Here, the intention is to enable AI researchers without optical expertise to solve the corresponding parameterized Markov decision processes based on common code in the future.
7 Conclusion
In this technical note, the comprehensive Python package TMMFast for multilayer thinfilm computation, optimization and reinforcement learning is presented. At its core, the package is comprised of revised and speedup transfer matrix method code from the original TMM package[4]. TMMFast enables the user to compute multilayer thinfilms with Numpy and PyTorch methods and gives full control over the data to the user. Since the TMMFast package evaluates multilayer thinfilms especially fast, it can be used to generate datasets for machine learning. The reduced computational time also enables evolutionary optimization to be executed in a reasonable amount of time. The TMMTorch implementation allows the user to compute analytical gradients via automatic differentiation. QuasiNewton or gradient based optimization algorithms such as gradient descent can be used and converge faster to local minima by using an analytical gradient. Finally, an OpenAi Gym environment is proposed which allows researchers to easily test and experiment with new reinforcement agents on solving the multilayer thinfilm problem. All code proposed in this paper is opensource and available at https://github.com/MLResearchAtOSRAM/ under the MITlicence.
References
 [1] (1950) La théorie générale des couches minces. Journal De Physique Et Le Radium 11, pp. 307–309. Cited by: §1.
 [2] (2012) Numerical methods for the design of gradientindex optical coatings. Appl. Opt. 51 (34), pp. 8277–8295. External Links: Link, Document Cited by: §1.
 [3] (2014) Design and realization of advanced multiindex systems. Appl. Opt. 53 (4), pp. A88–A95. External Links: Link, Document Cited by: §1.
 [4] (2019) Multilayer optical calculations. External Links: 1603.02720 Cited by: §1, §7.
 [5] (1990) Optimization of a thinfilm multilayer design by use of the generalized simulatedannealing method. Optics Letters. Cited by: §1.
 [6] (2018) Design and optimization of thin film polarizer at the wavelength of 1540 nm using differential evolution algorithm. Optical and Quantum Electronics 50 (192). Cited by: §1.
 [7] (2014) Design of broadband omnidirectional antireflection coatings using ant colony algorithm. Optics Express 22. Cited by: §1.
 [8] (202009) Array programming with NumPy. Nature 585 (7825), pp. 357–362. External Links: Document, Link Cited by: §1.

[9]
(2019)
Accelerating optics design optimizations with deep learning
. Optical Engineering 58. Cited by: §1.  [10] (2020) Multilayer optical thin film design with deep q learning. Sci Rep 10. Cited by: §1, §6.
 [11] (1998) Gradientbased learning applied to document recognition. In Proceedings of the IEEE, pp. 2278–2324. Cited by: Figure 1, §5.
 [12] (1981) Computeraided techniques for the design of multilayer filters. A. Hilger. Cited by: §1.

[13]
(1995)
Synthesis of optical multilayer systems using genetic algorithms
. Applied Optics 34. Cited by: §1. 
[14]
(2016)
Reinforcement learning with parameterized actions.
In Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence
(), pp. 1934–1940. Cited by: §1, §6.  [15] (1965) A simplex method for function minimization. Comput. J. 7, pp. 308–313. Cited by: §5.
 [16] (2017) Automatic differentiation in pytorch. Cited by: §1, §5.
 [17] (2013) Genetic algorithms, a natureinspired tool: a survey of applications in materials science and related fields: part ii. Materials and Manufacturing Processes 28. Cited by: §1.

[18]
(2018)
Modeling and optimization of thinfilm optical devices using a variational autoencoder
. Technical report Stanford University. Cited by: §1.  [19] (1998) Introduction to reinforcement learning. Cambridge: MIT Press. Cited by: §6.
 [20] (1993) Quasioptimal synthesis for antireflection coatings: a new method. Appl. Opt. 32 (22), pp. 4265–4275. External Links: Link, Document Cited by: §1.
 [21] (1993) Some theoretical aspects of thinfilm optics and their applications. Appl. Opt. 32 (28), pp. 5417–5426. External Links: Link, Document Cited by: §1.
 [22] (2021) Parameterized reinforcement learning for optical system optimization. J. Phys. D: Appl. Phys 54 (30), pp. . Cited by: Figure 1, §1.
 [23] (198905) Learning from delayed rewards. Ph.D. Thesis, King’s College, Cambridge, UK. External Links: Link Cited by: §6.

[24]
(2013)
Design of reflective color filters with high angular tolerance by particle swarm optimization method
. Opt. Express 21. Cited by: §1.
Appendix A Appendix
a.1 TMMFast Example
To demonstrate the functionality, a minimal example is given for a multilayer thinfilm with random thicknesses and refractive indices.
The result of the computation can then be further evaluated. By using the plot_stacks function from the TMMFast package, the multilayer thinfilm can be visualised. An example is shown in Figure 2.
a.2 Dataset Generation
A dataset containing tens of thousands of datasamples can be created by using the multithread_coh_tmm function. The external package tqdm is used to display a progress bar in order to keep track of the generation progress. In this example, a dataset with 1e6 thinfilms at 100 wavelengths should be created for vertical incidence.
Now the dataset can be saved for example by using .hdf or Numpy’s .npz format.
a.3 TMMTorch Example
Here, an example how to compute gradients via automatic differentiation with the TMMTorch package is shown.
Since the error is computed with respect to zero reflectivity for all wavelengths at all incidences, the gradient points towards the steepest descent for a broadband antireflection coating. By using the minimize function from the Python package Scipy and a gradient based optimization algorithm, for example "LBFGSB"^{6}^{6}6https://docs.scipy.org/doc/scipy/reference/optimize.minimizelbfgsb.html, a local optimum is easily found.