TMM-Fast: A Transfer Matrix Computation Package for Multilayer Thin-Film Optimization

11/24/2021
by   Alexander Luce, et al.
2

Achieving the desired optical response from a multilayer thin-film structure over a broad range of wavelengths and angles of incidence can be challenging. An advanced thin-film structure can consist of multiple materials with different thicknesses and numerous layers. Design and optimization of complex thin-film structures with multiple variables is a computationally heavy problem that is still under active research. To enable fast and easy experimentation with new optimization techniques, we propose the Python package TMM-Fast which enables parallelized computation of reflection and transmission of light at different angles of incidence and wavelengths through the multilayer thin-film. By decreasing computational time, generating datasets for machine learning becomes feasible and evolutionary optimization can be used effectively. Additionally, the sub-package TMM-Torch allows to directly compute analytical gradients for local optimization by using PyTorch Autograd functionality. Finally, an OpenAi Gym environment is presented which allows the user to train reinforcement learning agents on the problem of finding multilayer thin-film configurations.

READ FULL TEXT VIEW PDF

page 3

page 7

08/16/2019

Studying Wythoff and Zometool Constructions using Maple

We describe a Maple package that serves at least four purposes. First, o...
11/06/2019

DISROPT: a Python Framework for Distributed Optimization

In this paper we introduce DISROPT, a Python package for distributed opt...
06/24/2021

Graphlets in multilayer networks

Representing various networked data as multiplex networks, networks of n...
05/05/2022

Perron communicability and sensitivity of multilayer networks

Modeling complex systems that consist of different types of objects lead...
07/09/2021

Universal Multilayer Network Exploration by Random Walk with Restart

The amount and variety of data is increasing drastically for several yea...
08/23/2018

Multilayer Optimization for the Quantum Internet

We define a multilayer optimization method for the quantum Internet. Mul...
10/29/2020

Multilayer Clustered Graph Learning

Multilayer graphs are appealing mathematical tools for modeling multiple...

1 Introduction

Although the existence of a globally optimal multilayer thin-film was proven [20, 21, 6], the optimization of multilayer thin-films 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 thin-films. At its core, the proposed TMM-Fast package (Transfer Matrix Method - Fast) is a parallelized and re-implemented 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 thin-film materials with different thicknesses. Given a broad, discretized spectrum of light irradiating a thin-film 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 thin-film. An intuitive approach to increase the computational speed of the response of a multilayer thin-film is to parallelize the pair-wise independent computations regarding wavelengths and angles. The parallelization is implemented via matrix operations based on the open-source numeric algebra package Numpy

[8] and the thread-management 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 thin-film 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 TMM-Fast code is used to implement an OpenAI Gym environment

111https://github.com/MLResearchAtOSRAM/gym-multilayerthinfilm

that can be used by physicists and reinforcement learning researchers. In the environment, the optimization of multilayer thin-films 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 GitHub222https://github.com/MLResearchAtOSRAM/tmm_fast under MIT license.

Figure 1: Schematic overview of the contents of the TMM-Fast package. The package consists of three subdivisions. Part a) of Figure 1 contains the core functionality of the TMM-package, computing the optical response of a multilayer thin-films quickly via the transfer matrix method over a broad range of wavelengths and incident angles. Another part of the core functionality of solving multilayer thin-films is the possibility to generate huge datasets with >1e5 data samples for machine learning models.
Part b) encompasses an OpenAI Gym Environment which allows easy comparisons of different reinforcement learning agents. Here, the environment state is given by the total layer number, layer thicknesses, and material choice. The agent takes an action that specifies the material and thickness of the layer to stack next. The environment implements the multilayer thin-film generation as consecutive conduction of actions and assigns a reward to a proposed multi-layer thin film based on how close the actual (solid orange line) fulfils a desired (dashed orange line) characteristic - in this case, a one-dimensional optical characteristic, e.g. reflectivity over wavelength for perpendicular incidence. The experience accumulated by the taken actions of the agent is used to adapt subsequent actions in order to increase the reward and thus generate more and more sophisticated multilayer thin-films. An investigation of a reinforcement learning agent which was trained with the proposed Gym Environment is given by Wankerl et. al. [22].

Part c) encompasses an implementation of the transfer matrix method via PyTorch functions which allow for backpropagation through the entire computation. This enables easy differentiation of the parameters of interest and allows gradient based and quasi-newton optimization methods to be used for optimization of multilayer thin-films. In the example shown, a random thin-film is optimized with respect to the given optimization target, weighted by the mean squared error. Gradients via autograd are used with the BFGS algorithm

[11]. The iterative evolution of the layer thicknesses and the correspondingly decreasing figure of merit are shown below. For the given initial values the optimization converged in about 25 iteration steps.

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 thin-film with layers and the individual layers are denoted by . The light enters from an injection layer of semi-infinite thickness with a relative amplitude of 1 and exits the multilayer thin-film in the outcoupling layer of semi-infinite thickness with . From the outcoupling layer, no light enters the thin-film. 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 thin-film is given by

(5)

Finally, to compute the reflection and transmission of the entire multilayer thin-film, 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 TMM-Fast 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 einsum333Documentation: https://numpy.org/doc/stable/reference/generated/numpy.einsum.html method allows to specify multiplication and contractions of different dimensions easily:

1M_l = np.zeros((num_lambda, num_angles, num_layers, 2, 2),
2                   dtype=complex)
3F = r_list[:, 1:]
4M_l[:,:,1:-1,0,0] = np.einsum(’hji,ji->jhi’, 1/A, 1/t_list[:, 1:])
5M_l[:,:,1:-1,0,1] = np.einsum(’hji,ji->jhi’, 1/A, F/t_list[:, 1:])
6M_l[:,:,1:-1,1,0] = np.einsum(’hji,ji->jhi’, A, F/t_list[:, 1:])
7M_l[:,:,1:-1,1,1] = np.einsum(’hji,ji->jhi’, A, 1/t_list[:, 1:])
8
9Mtilde = np.empty((num_angles, num_lambda, 2, 2), dtype=complex)
10Mtilde[:, :] = make_2x2_array(1, 0, 0, 1, dtype=complex)
11for i in range(1, num_layers-1):
12        Mtilde = np.einsum(’ijkl,ijlm->ijkm’, Mtilde, M_l[:,:,i])

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 TMM-Fast package provides compared to the native implementation of Byrnes, ten multilayer thin-films 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 thin-film in s which corresponds to an acceleration of 100.

Finally, the Python package Dask444https://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 TMM-Fast package, the computation is easily started. A Dataset with one million samples of 9-layer multilayer thin-films at 100 wavelengths and 10 angles of incidence can be created in approximately 40 min on an eight-core 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 TMM-Torch: Multilayer thin-film gradients via Autograd

For optimization, computing gradients enable gradient-based optimization algorithms such as gradient descent [11]. These gradient based optimization methods generally converge to local minima with fewer iterations than other non-gradient based algorithms such as the Nelder-Mead downhill simplex algorithm [15]. The Python package PyTorch555https://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 TMM-Torch subpackage of TMM-Fast, the gradients for a multilayer thin-film 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 TMM-Fast 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 thin-films. Therefore, the generation of multilayer thin-films 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 thin-film 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 pre-defined maximum number of layers is reached or the agent decides to terminate stacking. The optical characteristic of the final proposed multilayer thin-film of layers, e.g. regarding reflectivity over wavelength and angle of incidence, is computed via the proposed transfer-matrix method (TMM-Fast). 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 thin-film, 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 user-defined, 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 Q-learning [10, 23] — to distinguish between good and bad actions and thus derive an optimal thin-film design.

Whereas the contained physical methods are well-studied and known for decades, the contribution of the code available at https://github.com/MLResearchAtOSRAM/gym-multilayerthinfilm lies in the implementation of an OpenAI gym-related 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 TMM-Fast for multilayer thin-film 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]. TMM-Fast enables the user to compute multilayer thin-films with Numpy and PyTorch methods and gives full control over the data to the user. Since the TMM-Fast package evaluates multilayer thin-films 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 TMM-Torch implementation allows the user to compute analytical gradients via automatic differentiation. Quasi-Newton 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 thin-film problem. All code proposed in this paper is open-source and available at https://github.com/MLResearchAtOSRAM/ under the MIT-licence.

References

  • [1] F. Abelès (1950) La théorie générale des couches minces. Journal De Physique Et Le Radium 11, pp. 307–309. Cited by: §1.
  • [2] S. W. Anzengruber, E. Klann, R. Ramlau, and D. Tonova (2012) Numerical methods for the design of gradient-index optical coatings. Appl. Opt. 51 (34), pp. 8277–8295. External Links: Link, Document Cited by: §1.
  • [3] H. Becker, D. Tonova, M. Sundermann, H. Ehlers, S. Günster, and D. Ristau (2014) Design and realization of advanced multi-index systems. Appl. Opt. 53 (4), pp. A88–A95. External Links: Link, Document Cited by: §1.
  • [4] S. J. Byrnes (2019) Multilayer optical calculations. External Links: 1603.02720 Cited by: §1, §7.
  • [5] C. P. Chang, Y. H. Lee, and S. Y. Wu (1990) Optimization of a thin-film multilayer design by use of the generalized simulated-annealing method. Optics Letters. Cited by: §1.
  • [6] M. Ebrahimi and M. Ghasemi (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] X. Guo, H. Y. Zhou, S. Guo, X. X. Luan, W. K. Cui, Y. F. Ma, and L. Shi (2014) Design of broadband omnidirectional antireflection coatings using ant colony algorithm. Optics Express 22. Cited by: §1.
  • [8] C. R. Harris, K. J. Millman, S. J. van der Walt, R. Gommers, P. Virtanen, D. Cournapeau, E. Wieser, J. Taylor, S. Berg, N. J. Smith, R. Kern, M. Picus, S. Hoyer, M. H. van Kerkwijk, M. Brett, A. Haldane, J. F. del Río, M. Wiebe, P. Peterson, P. Gérard-Marchant, K. Sheppard, T. Reddy, W. Weckesser, H. Abbasi, C. Gohlke, and T. E. Oliphant (2020-09) Array programming with NumPy. Nature 585 (7825), pp. 357–362. External Links: Document, Link Cited by: §1.
  • [9] R. S. Hedge (2019)

    Accelerating optics design optimizations with deep learning

    .
    Optical Engineering 58. Cited by: §1.
  • [10] A. Jiang, Y. Osamu, and L. Chen (2020) Multilayer optical thin film design with deep q learning. Sci Rep 10. Cited by: §1, §6.
  • [11] Y. Lecun, L. Bottou, Y. Bengio, and P. Haffner (1998) Gradient-based learning applied to document recognition. In Proceedings of the IEEE, pp. 2278–2324. Cited by: Figure 1, §5.
  • [12] H. M. Liddell and H. G. Jerrard (1981) Computer-aided techniques for the design of multilayer filters. A. Hilger. Cited by: §1.
  • [13] S. Martin, J. Rivory, and M. Schoenauer (1995)

    Synthesis of optical multilayer systems using genetic algorithms

    .
    Applied Optics 34. Cited by: §1.
  • [14] W. Masson, P. Ranchod, and G. Konidaris (2016) Reinforcement learning with parameterized actions.

    In Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence

    (), pp. 1934–1940.
    Cited by: §1, §6.
  • [15] J. Nelder and R. Mead (1965) A simplex method for function minimization. Comput. J. 7, pp. 308–313. Cited by: §5.
  • [16] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer (2017) Automatic differentiation in pytorch. Cited by: §1, §5.
  • [17] W. Paszkowicz (2013) Genetic algorithms, a nature-inspired tool: a survey of applications in materials science and related fields: part ii. Materials and Manufacturing Processes 28. Cited by: §1.
  • [18] J. Roberts and E. W. Wang (2018)

    Modeling and optimization of thin-film optical devices using a variational autoencoder

    .
    Technical report Stanford University. Cited by: §1.
  • [19] Sutton and Barto (1998) Introduction to reinforcement learning. Cambridge: MIT Press. Cited by: §6.
  • [20] A. V. Tikhonravov and J. A. Dobrowolski (1993) Quasi-optimal synthesis for antireflection coatings: a new method. Appl. Opt. 32 (22), pp. 4265–4275. External Links: Link, Document Cited by: §1.
  • [21] A. V. Tikhonravov (1993) Some theoretical aspects of thin-film optics and their applications. Appl. Opt. 32 (28), pp. 5417–5426. External Links: Link, Document Cited by: §1.
  • [22] H. Wankerl, M. L. Stern, A. Mahdavi, C. Eichler, and E. W. Lang (2021) Parameterized reinforcement learning for optical system optimization. J. Phys. D: Appl. Phys 54 (30), pp. . Cited by: Figure 1, §1.
  • [23] C. J. C. H. Watkins (1989-05) Learning from delayed rewards. Ph.D. Thesis, King’s College, Cambridge, UK. External Links: Link Cited by: §6.
  • [24] C. Yang, L. Hong, W. Shen, Y. Zhang, X. Liu, and H. Zhen (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 TMM-Fast Example

To demonstrate the functionality, a minimal example is given for a multilayer thin-film with random thicknesses and refractive indices.

1import tmm_fast as tmmf
2L = 12 # number of layers
3d = np.random.uniform(20, 150, L)*1e-9 # thicknesses of the layers
4d[0] = d[-1] = np.inf # set first and last layer as injection layer
5n = np.random.uniform(1.2, 5, L) # random constant refractive index
6n[-1] = 1 # outcoupling into air
7
8wl = np.linspace(500, 900, 301)*1e-9
9theta = np.deg2rad(np.linspace(0, 90, 301))
10# here s and p polarization is computed and averaged to simulate incoherent light
11result = (tmmf.coh_tmm_fast(’s’, n, d, theta, wl)[’R’]
12                + tmmf.coh_tmm_fast(’p’, n, d, theta, wl)[’R’])/2

The result of the computation can then be further evaluated. By using the plot_stacks function from the TMM-Fast package, the multilayer thin-film can be visualised. An example is shown in Figure 2.

Figure 2: A multilayer thin-film with random layer thicknesses and refractive index under unpolarized illumination. In this example, the materials are dispersion- and dissipationless. The injection region, which is below the thin-film in this depiction, possesses a refractive index of . The outcoupling region above the thin-film possesses a refractive index of . The reflectivity is computed over a wavelength range of 500 nm to 900 nm and from 0 to 90 on a 300300 grid.

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 thin-films at 100 wavelengths should be created for vertical incidence.

1import numpy as np
2import tqdm
3
4n_samples = 1e6
5n_lambda = 100
6n_layers = 12
7wavelength = np.linspace(1000, 1700, n_lambda)*1e-9
8
9theta = np.array([0]) # vertical incidence
10stack_layers = np.random.uniform(5, 180, (n_samples, n_layers))*1e-9
11stack_layers[:,0] = stack_layers[:,-1] = np.inf # injection and outcoupling layer
12optical_index = np.array([2.5] +[2.0, 1.4]*5 +[1.])
13optical_index = np.tile(optical_index, (n_samples, 1))
14n = 10000 # the dataset is computed in steps of 1e5
15dataset = np.empty((n_samples, n_lambda))
16for i in tqdm.tqdm(np.array(range(n_samples))[::n]):
17    tmm_res = np.empty((n, n_lambda))
18    tmm_res = multithread_coh_tmm(’s’, optical_index[i:i+n], stack_layers[i:i+n], theta, wavelength, TorR=’R’).squeeze()
19    dataset[i:i+n] = tmm_res

Now the dataset can be saved for example by using .hdf or Numpy’s .npz format.

a.3 TMM-Torch Example

Here, an example how to compute gradients via automatic differentiation with the TMM-Torch package is shown.

1import torch
2import tmm_fast_torch as tmmt
3n_layers = 12 # number of layers
4stack_layers = np.random.uniform(20, 150, n_layers)*1e-9 # thicknesses of the layers
5stack_layers[0] = stack_layers[-1] = np.inf # set first and last layer as injection layer
6optical_index = np.random.uniform(1.2, 5, n_layers) # random constant refractive index
7optical_index[-1] = 1 # outcoupling into air
8
9stack_layers = torch.tensor(stack_layers, requires_grad=True)
10
11wl = np.linspace(500, 900, 301)*1e-9
12theta = np.deg2rad(np.linspace(0, 90, 301))
13
14result = tmmt.coh_tmm_fast(’s’, optical_index, stack_layers, theta, wavelength)[’R’]
15mse = torch.nn.MSELoss()
16error = mse(result, torch.zeros_like(result)
17error.backward()
18
19gradients = stack_layers.grad

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 anti-reflection coating. By using the minimize function from the Python package Scipy and a gradient based optimization algorithm, for example "L-BFGS-B"666https://docs.scipy.org/doc/scipy/reference/optimize.minimize-lbfgsb.html, a local optimum is easily found.