1. INTRODUCTION

In optics, multilayered structures have been studied extensively for more than a century. The very first effort at calculating the reflectance of multilayers dates back to the early work of Lord Rayleigh [1,2]. The subject has gained increasing attention in the middle of the 20th century when quarter-wave stacks (also called Bragg mirrors or 1D photonic crystals) were fabricated and fully understood for the first time [3]. An important contribution has been the formalism introduced by Abelès [4], which provided a deeper insight into the physical mechanisms behind the optical properties of multilayer systems, as well as a convenient mathematical framework for their calculation. In the years that followed, this formal approach was extended and successfully applied to periodic structures [5–7] and anisotropic media [8]. The resulting formalisms have been widely incorporated into optical textbooks [9–11].

The design of diverse optical filters, guided by physical intuition, emerged as early as 1958 [12], and thrived in the subsequent years, notably due to the works of Thelen [13,14] and Apfel [15]. The first computer program enabling the automation of optical filter design was published by Dobrowolski and Lowe in 1978 [16]. Given the versatility of multilayered structures, very efficient design methods have been conceived since then [17], requiring the reliable computing of the optical properties of increasingly complex filters. A direct application has been the wavelength division multiplexing technology, which allowed a substantial increase of the transcontinental optical connection bandwidth in the last decades. Today, optical filter design and fabrication contests [18–20] are periodically organized, showing how powerful state-of-the-art techniques have become. In the late 20th century, multilayered structures gained renewed interest and became widespread due to advancements in fabrication techniques. Notably, two types of problems have received particular attention due to their practical importance: sensors based on surface plasmon resonances for the detection of biologically relevant molecules [21], and anti-reflective coatings used in photovoltaic devices [22]. In parallel, the increasingly complex problems associated with plasmonics [23], metamaterials [24], or nonspecular phenomena [25] have also made multilayer simulations more challenging, particularly due to the frequent need to include a large number of metallic layers in the structures.

Over the years, an increasing number of numerical tools have thus been proposed to compute the optical response of complex multilayers [26–28], with a significant portion of them tailored for specific applications such as optical filters [29] or surface plasmon resonance (SPR) [30]. The demand for accurate and efficient computation of properties has grown significantly, while well-studied and understood structures serve as a testing ground for optimization techniques [31,32]. In the recent past, thanks to the fast computation methods, multilayer stacks have also often been used as large and cost-effective datasets for rapid testing of deep learning algorithms [33–39]. Due to its relative simplicity, differentiable transfer matrix implementations have furthermore been used as actual part of deep learning models for accelerated training [40,41]. Because of this multiplication of objectives, the range of computed quantities of physical interest has also expanded.

Here, we introduce PyMoosh (Python-based multilayer optics optimization and simulation hub), a user-friendly and comprehensive set of numerical tools to compute as efficiently and as reliably as possible many optical properties of multilayered structures [42,43]. This includes reflectance and transmittance in any condition and with different state-of-the-art numerical methods, but also absorption and Poynting flux in any layer, 2D Green functions, and the computation of structure modes and their dispersion curves. Finally, PyMoosh is specifically designed to be used together with optimization tools. Developed in Python, PyMoosh is an open-source software aimed at maximizing usability. It is intended to be simple enough to be used for educational purposes, while also offering advanced features for research, with a particular focus on tackling demanding problems such as optimization and deep learning. It is the successor of Moosh, which was implemented in Octave/MATLAB [44].

This paper has a total of five sections. In the first section, we detail the different techniques included in PyMoosh to solve Maxwell’s equations in multilayered structures. PyMoosh specifically includes transfer and scattering matrices, the Abelès formalism, Dirichlet-to-Neumann (DtN) maps and the admittance formalism. The methods, all gathered and described in a single place for educational and convenience purposes, are compared in terms of computational costs and tested for reliability on common photonics problems. In a second section, we detail how the Poynting flux and absorption in each layer can be computed and explain the tools specific to photovoltaic problems that are implemented in PyMoosh. In a third section, we explain how the field inside the structure can be computed with high accuracy using scattering matrices. We demonstrate this for two example illuminations, namely an incoming beam and an oscillating line current source, placed at an arbitrary location inside the structure. The fourth section describes how resonant modes (leaky or guided) can be retrieved. The last section describes our general approach in terms of coding and implementation, along with a discussion of the motivations behind our choices. PyMoosh comes in fact with a very large number of examples, under the form of Jupyter Notebooks [45], which are used to illustrate its capabilities, and that may serve as starting points for applications of other users.

2. COMPUTING REFLECTANCE AND TRANSMITTANCE

Historically, and for understandable practical reasons, the focus of multilayer structure studies has been put on computing the reflectance and eventually the transmittance of a given structure for a given wavelength $\lambda$. This is equivalent to solving Maxwell’s equations in the structure for a time dependency ${e^{- i\omega t}}$, where $\omega = \frac{{2\pi c}}{\lambda}$ is the angular frequency of the electromagnetic wave. Since the structure is invariant by translation along the $x$ and $y$ axis, we choose without any loss of generality, to consider a ${e^{i{k_x}x}}$ spatial dependency with respect to $x$ and invariance along $y$ for all fields. In transverse electric (TE) [resp. transverse magnetic (TM)] polarization, ${E_y}$ (resp. ${H_y}$) satisfies the Helmholtz equation

Because the ${k_x}$ spatial dependency is constant through the stack, in any layer $i$, the vertical ($z$) propagation constant of the field can be defined by

The solution to Eq. (1) can be written in layer $i$ as

 

Fig. 1. Multilayer structure and definition of important variables. We consider the light coming from the top of the image. ${\epsilon _i}$ and ${\mu _i}$ are, respectively, the permittivity and the permeability of layer $i$. ${A_i}$ (resp. ${B_i}$) are defined as the amplitude of the upward propagating (resp. downward propagating) field in layer $i$, and the ${+}$ (resp. ${-}$) superscript indicates that the field value is taken at the top (resp. bottom) of layer $i$.

Download Full Size | PDF

The continuity conditions for the tangential electric and magnetic fields at each interface between layer $i$ and layer $i + 1$ lead to

When the amplitude of the incident plane wave is taken to be equal to 1 ($B_0^ - = 1$) and there is no incoming light from below ($A_{N + 1}^ - = 0$), then the reflection and transmission coefficient are $A_0^ - = r$ and $B_{N + 1}^ + = t$, respectively. This means we have a system of $2(N + 1)$ continuity equations and $2N$ propagation equations for just as many unknowns (the four $A_i^ \pm$ and $B_i^ \pm$ coefficients of each layer $i \in [1;N]$, $A_0^ -$ and $B_{N + 1}^ +$). In theory, such a system could be solved numerically using classical algorithms. However, the corresponding matrix is sparse and contains figures with very different magnitudes so that classical approaches are numerically unstable, to the point that it often cannot be used for more than typically a few layers. Specific formalisms to solve this system of equations have therefore been developed and are presented below.

The TM polarization can be treated simply by changing ${\mu _i}$ into ${\epsilon _i}$ in the system above. We underline, however, that there is a $\pi$ phase in normal incidence for the reflection coefficient defined in TE and in TM polarization.

Now, because of Eq. (2), we need to discuss how the square root of a complex number is defined and how this may have an impact on the accuracy of the calculation. By default, on any computer, the phase of a complex number $\zeta$ is defined in ]${-}\pi ,\pi$], meaning that taking its square root thus produces a complex number with an argument in ]${-}\pi /{2},\pi /{2}$] (i.e., with a positive real part and arbitrary imaginary part). This determination of the square root can be changed manually, to obtain a positive imaginary part and arbitrary real part for instance, by taking the opposite of the result of the square root if the argument of $\zeta$ is negative. In any case, this introduces a cut in the complex plane, which cannot always be ignored.

Here, Eq. (3) indicates that changing the determination of the square root for $\gamma$ is equivalent to exchanging $A$ and $B$ for all the layers in the structure, so that this is not expected to change the result of the computation, as Eqs. (4) and (5) are unchanged by this permutation. However, it is clear that if ${\gamma _i}$ has a negative imaginary part, then these equations introduce positive real exponentials, which would amplify any numerical error made on their amplitude. We have made the choice to use ${\gamma _i}$ with positive imaginary parts in PyMoosh to maximize the stability of the code by using only decreasing exponential functions in all our calculations (especially the computation of the whole field; see below), inspired by what is done in RCWA [47,48]. The change we have brought to the determination of the square root for the surrounding medium has more importance and will be discussed later.

From there, there are different ways to eliminate all the unknowns to obtain the reflection or the transmission coefficients. The most common formalism is based on transfer matrices ($T$-matrices). However, a more physical and stable formalism consists in writing reflection and transmission coefficients for each layer or interface. This is the scattering matrix ($S$-matrix) method, widely used in more advanced simulation techniques where evanescent waves are ubiquitous. Since the work of Abelès [4], we also know that considering the impedance or admittance of each layer is a relevant approach. This approach can be translated both into a matrix formalism and into a recursive formula to compute how the admittance of a structure is modified by adding a supplementary layer [11]. Because these methods are not perfectly stable, DtN maps have been proposed more recently [49]. The DtN technique is related to the Abelès matrices formalism exactly in the same way that the scattering matrix formalism is related to the transfer matrix approach. In the following, we derive the main formalisms of these different approaches.

A. $T$-Matrices

Let us define ${\psi _i} = \frac{{{\gamma _i}}}{{{\mu _i}}}$ for all layers. For the interface between layers $i$ and $i + 1$, the continuity conditions in Eqs. (4) and (5) can be rewritten as

From the transfer matrix,

B. $S$-Matrices

Scattering matrices are simply another way to eliminate unknowns. The $S$-matrix formalism uses the same field definitions but rewrites the continuity conditions obtained from Eqs. (4) and (5) to link incoming fields at the interface $(A_i^ + ,B_{i + 1}^ -)$ with the outgoing fields $(A_{i + 1}^ - ,B_i^ +)$, giving the new interface matrices

The propagation matrices are also modified to

As we will see later, these propagation matrices, where both exponentials are identical, lead to a greater numerical stability. However, computing the matrix linking boundary conditions is no longer a straightforward matrix product, since the coefficients of a given layer are not on the same side of Eq. (10). Therefore, the relationship between the different $I$ and ${C^{(S)}}$ matrix equations can be written as

We can define a “cascading” operation, computing $S$ as a function of $U$ and $V$ coefficients

Applying this computation iteratively on every interface and layer leads to the scattering matrix $S$ for the whole structure, written as

The reflection and transmission coefficient can be read directly, knowing $B_0^ - = 1 ,A_{N + 1}^ + = 0$: $r = {S_{00}} ,t = {S_{10}}$.

C. Abelès Formalism

Contrary to ${T}$-matrices and ${S}$-matrices, the Abelès formalism (as well as the following DtN formalism) uses only one ${E_y}$ field position in each layer, thus removing the need for propagation matrices. We will also therefore directly speak in terms of ${E_{y,i}}$ and $\partial z{E_{y,i}}$. Taking, for instance, the field on the upper side of the layers and using the field continuity, we can write

These relations and the continuity conditions at the interface at ${z_{i + 1}}$ lead to the matrix equations

Once the matrix product $M = \prod\nolimits_i {M_i}$ is computed, the $r$ and $t$ coefficients are slightly more complicated to retrieve than for $T$-matrices and are written as

D. DtN Formalism

Finally, just like the ${S}$-matrix formalism was the equivalent of ${T}$-matrices but using the incoming field on one side of the equation and the outgoing fields on the other, the DtN formalism is the equivalent of the Abelès formalism but keeping ${E_y}$ on one side of the equation and $\partial z{E_y}$ on the other, written as

As with ${S}$-matrices, the complete structure matrix computation is not straightforward, and needs a cascading operation similar to the one used for ${S}$-matrices. To this end, if we call $A = {N_{i,i + 1}}$ and $B = {N_{i + 1,i + 2}}$, then the ${D_{i,i + 2}}$ matrix linking fields in layer $i$ and in layer $i + 2$ is given by

For all other formalisms, we typically consider a superstrate of thickness 0, to have a phase reference exactly at the interface between the superstrate and the first layer of the structure. However, for the DtN formalism, as we can see in Eq. (22) this would lead to division by 0 in the first matrix. So, we choose an arbitrary thickness for the superstrate ($\lambda /100$ in the PyMoosh code) and compensate for the additional phase at the end of the computation.

E. Admittance Formalism

The admittance formalism is based on the Abelès formalism and its main advantage is to make computations of the structure’s reflectivity much faster. We define

With the Abelès formalism [inverting Eq. (17)], we get

We then define $X = \frac{{{\partial _z}{E_y}}}{{{E_y}}}$, and get

Computing the matrix-vector product and dividing both rows, we get

Let us define ${Y_i} = \frac{i}{{{k_0}}}{X_i}$, leading to

Then, we use ${Y_N} = {n_{{\rm eff},N}}$ as the initial condition and compute all ${Y_i}$, with the final reflection coefficient being for TE polarization

In TM, the difference is that $n_{{\rm eff},i}^{\rm (TM)} = n_i^2/{n_{{\rm eff},i}}$ and

This formalism allows to compute only the reflectance of a structure, and is thus much less general than the methods described above. It can be, however, of particular interest for pedagogical reasons or in time-critical applications, given its simplicity and, as a direct consequence, its computational efficiency.

F. Comparing the Different Methods

Choosing the right method to compute the reflectance or transmittance of a structure depends on the context and whether stability or computational speed are particularly sought after. No method is both the fastest and the most stable. We have thoroughly tested the stability and estimated the efficiency of each of the methods as implemented in PyMoosh, to assist users in determining the most suitable approach for their specific needs.

1. Computational Speed

As can be guessed from the mathematical expressions above and as is shown Fig. 2, the computational time required for computing the reflectance of a given structure increases linearly with the number of considered layers. Also, as can be expected, the scattering matrix method is the slowest. It is typically twice as slow as the ${T}$-matrix or DtN methods. On the other hand, the Abelès formalism is twice as fast, but the most computationally efficient technique of all is the admittance formula. It is thus ideally suited for optimizing the reflectance of a multilayer. Alas, this technique is also the most prone to stability issues.

 

Fig. 2. Computation time of the $r$ and $t$ coefficients for all methods as a function of the number of layers in a stack, averaged over 500 runs. As exact computation times are heavily influenced by hardware capabilities, it remains important to prioritize the relative comparison between formalisms over these precise times.

Download Full Size | PDF

2. Numerical Stability

For a random dielectric structure, all formalisms can be considered stable, as the transmission can be expected to be relatively high. Numerical accuracy issues occur when evanescent waves are present, whether they are due to total internal reflection, the presence of metals, or if evanescent Bloch waves are excited in the bandgaps of Bragg mirrors. such Although such problems are unlikely for optical filters, evanescent waves are ubiquitous in prism couplers or plasmonic structures, not to mention photonic crystal-based structures.

The scattering matrix method has proven to be the most stable and accurate of all, regardless of the situation [28,47,48]. For this reason, we consider its results as our reference in all the following. Figures 3–5 show the absolute difference between the values given by the scattering matrix method and the other methods for the reflection and transmission coefficient in several cases. We explicitly discuss here the computation errors for TE illumination, but TM polarization faces roughly the same instabilities.

 

Fig. 3. (a) Absolute error between S-matrix (considered as our reference) and other formalisms. Transmission coefficient for a quarter-wave stack of materials $n = 1.2$ and $n = 1.5$ as a function of the total number of layers of the stack. (b) Schematic representation of the structure.

Download Full Size | PDF

 

Fig. 4. (a) Absolute error between the $S$-matrix (considered as our reference) and other formalisms. Transmission coefficient in a frustrated total internal reflection situation for a layer of air of width $D$ between materials of index $n = 1.5$. (b) Schematic representation of the structure.

Download Full Size | PDF

 

Fig. 5. (a) Absolute error between the $S$-matrix (considered as our reference) and other formalisms. Surface plasmon resonance coupler in the Kretschmann–Raether configuration with a prism of index $n = 1.5$ coupling to air through a silver layer of width $D$. (b) Schematic representation of the structure.

Download Full Size | PDF

The first case, depicted in Fig. 3, is a simple Bragg mirror, illuminated at an oblique incidence (15°) with a wavelength in the middle of the bandgap. As the number of layers increases, the transmission coefficient goes down to zero. As can be seen, the DtN and scattering matrix methods differ only by an error of the order of machine precision. Therefore, when the relative error stays the same, as the value of the transmission coefficient decreases, so does the absolute error. It can be seen however that the Abelès formalism and the $T$-matrix method are subject to numerical instability, which typically occurs for more than 150 layers. This means that these methods should be considered with caution when structures involve Bragg mirrors, for instance in microcavities [51] for which the transmission through the stack is crucial. It should also be noted that the TE case is slightly more stable than the TM case, when no magnetic materials are involved, and that the Abelès formalism is significantly more stable under normal incidence illumination.

The second case, depicted in Fig. 4, is a simple frustrated total internal reflection, i.e., a layer of air comprised between two layers of high index material $n = 1.5$, illuminated with a 42° incidence angle, slightly beyond the critical angle. When the air layer thickness is increased beyond a few wavelengths, the evanescent field within causes numerical instabilities for the ${T}$-matrices and the Abelès matrices, similar to what occurs for the Bragg mirror. This means these two algorithms should be avoided when simulating prism couplers, where evanescent waves occur. Above a thickness of 2 µm for the air layer, as the theoretical transmission tends toward zero, the DtN and ${S}$ matrix methods often agree down to machine precision, leading to a zero absolute numerical error. They differ only by a relative error of the order of machine precision in other cases.

As a final example, we consider a metallic (silver) layer comprised between an upper layer with a high index ($n = 1.5$) and air, which is a typical prism coupler for a surface plasmon in a Kretschmann–Raether configuration [52]. In the previous cases, we underline that the reflection coefficient is computed with no instability for ${ T},\;{S}$, Abelès, or DtN matrices, even when the transmission coefficient is unstable. However, as can be seen in Fig. 5, the admittance formula is not accurate for thin metallic layers, and while the instability is not exponential, the value computed is not at all accurate. This shows that the admittance formalism is most adapted to optical filters with no metallic layers. While this cannot be seen in the linear scale graph, shown in Fig. 5(a), the difference between the values computed using the three other methods is down to machine precision.

Despite their stability in these typical examples, ${S}$ and DtN matrices also have specific breaking points. For scattering matrices, when ${\gamma _i}/{\mu _i} \sim - {\gamma _{i + 1}}/{\mu _{i + 1}}$ (in TE polarization, but this also occurs in TM with ${\epsilon _i}$ instead of ${\mu _i}$), a division by zero occurs. This can only happen when lossless metamaterials are considered, which is not realistic, but is often studied for educational purposes (showing the ideal response of Pendry’s lens, for instance [53]). As for the DtN formalism, its better stability is due to the terms in $1/\!\sin$ and $1/\!\tan$, which are less prone to numerical errors than exponentials. However, there are still specific breaking points. For instance, for a layer of optical thickness $\lambda /2n$, where $n$ is the refractive index of the layer, in normal incidence, the argument of the $1/\!\sin$ is $\pi$, leading to a division by 0. As this may occur with any kind of transparent material, and it is typically used as a cavity in distributed-feedback lasers, for instance, the DtN formalism must be used with care.

3. Discussion

From the results above, it seems clear that scattering matrices should be the default choice when stability matters because it presents only very specific breaking points, and we underline that the supplementary computational effort is reasonable. Transfer matrices seem an obvious choice in the literature; however, our results show that the Abelès formalism is equivalently stable but twice as computationally efficient, while being no more complex to implement. When speed really matters, when no metal is present, and when the computation of reflectance is sufficient, the admittance formalism is a less explored but very interesting alternative, for instance, in the context of optimization [31] or for database generation for deep learning, being one order of magnitude more efficient than $S$-matrices. Finally, the DtN formalism is a good compromise for large multilayer stacks that may exhibit evanescent fields, as long as no $\lambda /2n$ layer is present. PyMoosh allows such choices to be made and to easily assess the accuracy of any method by comparing it to the $S$-matrix result.

3. COMPUTING ABSORPTION IN AN ARBITRARY LAYER

In many applications (e.g., sensing or photovoltaic stacks), knowing where the absorption takes place is of paramount importance. This requires computation of the field amplitudes $A_i^ \pm$ and $B_i^ \pm$ that are ignored when only reflectance is to be computed. This can be done in any of the four main formalisms, and the present section describes how these fields can be computed for the ${S}$-matrix and the Abelès formalisms, because they have been shown to be, respectively, the most stable and the fastest. The admittance formalism can obviously not be used because it is designed to only compute the external coefficients.

To compute how much of the incident power is absorbed in each layer, the PyMoosh library computes the (time-averaged) projection along the $z$ axis of the Poynting vector at the top of every layer $i$, ${N_i}$. This projected Poynting vector is proportional to the power flux through the top interface of the layer. Since we normalize all power by the incident power, then the difference in power given by the difference in Poynting vector ${{\cal A}_i} = {N_i}$ - ${N_{i + 1}}$ is the fraction of the incident power absorbed in layer $i$.

A. Abelès Matrices

With ${M_i}$ defined in Eq. (17), we can define $M_{\rm{tot}}^{(k)} = \prod\nolimits_{i = 1}^k {M_{i - 1}}$, which allows us to compute the fields in the structure by

Knowing this, the equation used to compute the Poynting vector at the top of one layer is

B. $S$-matrices

Using the cascading operation defined in Eq. (14), we can introduce the $\odot$ operator, where $A \odot B$ is the result of the cascading operation. This allows the partial scattering matrices ${U^{(k)}}$ and ${D^{(k)}}$ to be written more simply, corresponding to the scattering matrices of the incomplete structures

The same can be written for the lower coefficients using the odd index $2i + 1$. The idea here is that if we assume an incoming plane wave from above with a unity amplitude (i.e., $B_0^ - = 1$ and $A_N^ + = 0$), we can use the intermediate matrix that usually appear when calculating the cascading formula to retrieve the unknowns

Again, using the odd index $2i + 1$ instead of $2i$ allows the retrieval of the $A_i^ -$ and $B_i^ -$ coefficients.

The equation used to compute the Poynting vector at a given interface, using the corresponding coefficients, is

C. Applications to Solar Cells

Solar cells can most often be modelized as a multilayered structure. Once the absorption in the active part ${{\cal A}_{\rm{active}}}$ is computed for the whole spectrum, it can be converted into an estimate of the short-circuit current ${j_{\rm{sc}}}$ delivered by the cell under standard illumination (an AM1.5 solar spectrum) assuming a unity quantum yield and using [54]

Because of how important this type of application is, there are functions in PyMoosh specifically designed to simulate photovoltaic multilayers, considering the photon density of the solar spectrum and the conversion efficiency of solar cells.

Figure 6(a) shows the photon density of the visible solar spectrum, received on the Earth, as included in the PyMoosh library. Figure 6(b) shows the absorptance of two structures: one that is a single layer of amorphous silicon (aSi) of 1 µm thickness, and a second structure with an additional dielectric anti-reflecting (AR) layer on top ($n = 1.5$, $h = \lambda /4n$, $\lambda = 600\,{\rm nm}$). These absorptance curves show that bulk aSi already has an absorptance of about 50%, and that a single AR coating can improve this absorptance to more than 85% for half of the visible spectrum.

 

Fig. 6. (a) Photon spectral density (photons per wavelength unit) of the AM1.5 visible solar spectrum. (b) Absorptance in the active layers of two different basic solar cells ($S$-matrix). 1 µm of SiA bulk with or without adding a layer of optical index $n = 1.5$ and thickness $\lambda /4n$ with $\lambda = 600\,{\rm nm}$.

Download Full Size | PDF

4. FIELD COMPUTATION

While reflectance or the transmittance spectra are among the most represented quantities, it is highly informative to be able to compute the electric or magnetic fields inside a structure in response to an illumination, either with an incoming beam, or when a source is placed in the structure. The magnetic field is especially relevant for all situations involving metals. The computational cost is less of a concern in such a situation, while the numerical stability is crucial. We thus use the ${ S}$-matrix method in that context, and the technique explained above to compute all the intermediate coefficients $A_i^ \pm$ and $B_i^ \pm$. The electric (resp. magnetic) field in TE (resp. TM) polarization can be written, assuming $x = 0$ and for a single incident plane wave, as

Using this expression and the determination of the square root mentioned above (positive imaginary part), presents the advantage that, since all the coefficients are computed with the best possible accuracy, all the exponentials are decreasing away from the interfaces. No numerical error can thus be amplified.

A. Illumination with a Gaussian Beam

A Gaussian beam can be written as a superposition of plane waves. In TE polarization, the ${E_y}$ field in layer $i$ can be written as

In Pymoosh, the field is computed for each given ${k_x}$ in all the layers using Eq. (39) and for all the points defined in the vertical discretization. Then, the full field at different horizontal positions is computed by multiplying this vertical vector by ${e^{i{k_x}x}}$ for all the values of $x$ included in the horizontal discretization. The computed fields for each plane wave ${k_x}$ are subsequently multiplied by the corresponding complex amplitude ${{\cal E}_0}({{k_x}})$ [see Eq. (41)] and all the plane wave contributions are summed, leading to the calculation given by Eq. (40). Importantly, as the number of incidence angles (and thus ${k_x}$ values) is finite and thus discretized, the resulting field is consequently periodic. The period is what is called the window size $d$. If $\Delta {k_x}$ is the difference between two successive values of ${k_x}$, we have the relation $d = \frac{{2\pi}}{{\Delta {k_x}}}$. The number of plane waves is chosen such that the amplitude for the maximum and minimum wavevectors given by Eq. (41) is smaller than 1% and has thus no impact on the final field image.

The exact same procedure can be applied to obtain images corresponding to ${H_x}$ and ${H_z}$ in TE polarization, calculated as

 

Fig. 7. Schematic representation of the field values surrounding the source.

Download Full Size | PDF

B. Current Line Source

It is relatively simple to compute the field corresponding to a source of oscillating current $I$, invariant in the $y$ direction, in TE polarization (the electric field, hence the oscillation direction along the infinite axis), even in a complex multilayered structure. The equation that must be solved inside the homogeneous layer where the source is placed can be written as

Now, it is always possible to write that $B_{{i_s}}^ + = {r_u}A_{{i_s}}^ +$ and $A_{{i_s}}^ - = {r_d}B_{{i_s}}^ -$, whatever the surrounding of the source. The coefficients ${r_d}$ and ${r_u}$ are simply found in the scattering matrix corresponding to the part of the structure below (respectively above) the source, as in Eq. (36). A simple representation of these fields and the contribution of the line source is presented in Fig. 7. A straightforward calculation yields an expression for the coefficients representing the wave amplitudes traveling into directions, respectively, upward ($A_{{i_s}}^ +$) and downward ($B_{{i_s}}^ -$) from the source, written as

These two coefficients represent the amplitudes of the plane waves excited by the source, and, as can be clearly seen, both depend highly on the surroundings through ${r_d}$ and ${r_u}$, which means the surroundings do influence the source itself. This allows the computation of the field above and under the source when no other incoming light is considered. Practically, we use Eq. (40) with amplitude ${{\cal E}_0}({k_x})$ unity for all the values of ${k_x}$.

5. GUIDED MODES TRACKING

It is extremely useful to be able to find the electromagnetic modes of a multilayered structure, whether these modes are strictly guided or leaky. Both types of modes correspond to a solution of the system of equations evoked above [see Eqs. (4) and (5)] without any illumination ($B_0^ - = A_N^ + = 0$), meaning the determinant of the system must be zero [57] (i.e., a zero of the dispersion relation).

The dispersion relation also corresponds to a zero of the determinant of ${S^{- 1}}$ or a pole of $S$, where $S$ is the scattering matrix of the whole structure. We underline, however, that this does not work for a single interface and the surface plasmon modes. Note that the determinant of an interface scattering matrix is always equal to ${-}{1}$.

A solution of the dispersion relation is also a pole of the reflection coefficient in the ${k_x}$ complex plane because it corresponds to the existence of a solution or an outgoing wave without any incoming plane wave, and thus a diverging reflection coefficient. Obtaining the guided modes in this way works in the general case, including the case of only a single interface. It is therefore the approach implemented in PyMoosh.

The solutions of the dispersion relation, and thus the poles of the reflection coefficient, can generally be found either for a real value of the wavevector ${k_x}$ and a complex value of the frequency $\omega$ (the imaginary part representing the typical lifetime of the mode) or for a real value of the frequency $\omega$ and a complex value of the wavevector ${k_x}$ (the imaginary part representing the typical propagation distance of the mode before leaking). In PyMoosh we have, so far, made the latter choice. However, the default choice for the computation of the complex square root is not adequate for finding the solution of the dispersion relation. Given the physical meaning of the imaginary part of the poles, it is, most of the time [56], positive. The determination of the square root has to be modified to make the poles appear [58] and the default choice is to consider that the phase of ${k_z}$ should be comprised between ${-}\pi /5$ and $4\pi /5$. This allows, in most cases, the poles to be seen correctly and the complex plane close to them to be explored.

 

Fig. 8. (a) Dispersion of the guided modes of a 2000 nm thick slab of dielectric with a refractive index $n = 1.5$ in air, plotted as the effective wavelength of the mode as a function of the wavelength in vacuum in a manner similar to [59]. (b) Corresponding field profiles of the six modes with lowest effective index.

Download Full Size | PDF

To find the guided modes, we define a minimum and a maximum effective index, defining bounds on the wavevectors. Then we start from a set of points regularly placed on the real axis between these values for ${k_x}$ and run a steepest descent over $\frac{1}{{|r|}}$ using a finite difference scheme to estimate the gradient of the function. Steepest descent algorithms are typically employed to locate any minima, but here we leverage the fact that the minima correspond to zeroes, to determine when the algorithm should stop. Several of the initial points will converge to the same pole, and only the unique poles are kept. While this technique is probably not the most efficient nor the most advanced, it has proven relatively reliable so far.

Once the modes have been found, their profiles can be computed. An evanescent wave with a unity amplitude at the first interface is assumed. The rest of the mode is reconstructed using scattering matrices. Any inaccuracy in the determination of ${k_x}$ will result in a discontinuity at the first interface. Figure 8(b) presents the modes found in a 2000 nm thick slab of dielectric with a refractive index $n = 1.5$ in air.

It is often practical to obtain the dispersion curve (i.e., the values of the wavevector of a given mode for different frequencies). In PyMoosh, once the modes have been found for the largest considered frequency, they can be tracked using steepest descents while the frequency changes. An example is shown in Fig. 8(a). We begin with the smallest wavelength so that no mode appears during the sweep; however, modes do disappear when they reach a cut-off.

6. CODE DESIGN PRINCIPLES

PyMoosh has been developed with the core principles of open science in mind. The choice of Python as the programming language plays a central role: Python packages are extremely easy to find and install through the use of the PyPI repository [60]. Python is often called good glue, thanks to its ability to interface other programming languages seamlessly. Moreover, modern repositories such as GitHub [61], being partly social networks, combined with the GNU General Public License, allow the accessibility and reusability of the code to be maximized. Our approach thus aligns with the principles of the FAIR data framework, which underscore the importance of making data and methods findable, accessible, interoperable, and reusable, which is the FAIR data approach [62].

The code architecture was chosen with simplicity of use in mind, especially for researchers and physicists in photonics. Although Python is an object-oriented language, we intentionally minimized the use of objects to facilitate code comprehension for nonprogramming experts. For convenience, some parts such as material handling are implemented as classes, but we emphasize that the part concerning materials can be used and expanded upon completely independently from the rest of the code. Finally, due to the inherent ease of implementation in Python, some interesting numerical optimizations that could be deemed fundamental (e.g., parallel computing) are left entirely in the hands of the user.

Our code is highly versatile, which may initially disorient some users, especially those who are not familiar with Python. To address this, in addition to this tutorial, we have provided extensive documentation and illustrated various use cases using Jupyter Notebooks (a beginner’s guide to simulation and another for optimization) and then specific use cases: surface plasmon resonance, total internal reflection, photovoltaics devices, ellipsometry problems, Green function computation, use of various formalisms, guided mode representation and tracking, Prism couplers, photonic crystals, and managing materials. A sample of the results these notebooks provide is shown in Fig. 9. Jupyter Notebooks have gained significant popularity, particularly in the context of open science [45], making them an accessible and user-friendly tool for exploring and utilizing our code. Finally, we illustrate possible applications for deep learning training dataset generation by a set of open source, extensive Python notebook examples, which are exhaustively discussed in a recent tutorial article [63,64].

 

Fig. 9. PyMoosh use case examples. (a) and (b) Ellipsometry results $\frac{{{r_p}}}{{{r_s}}} = \tan (\phi){e^{{i\Delta}}}$ for a layer of $n = 1.33$ and width 400 nm on top of a gold substrate. (c) Mode profiles for two coupled waveguides of index $n = 1.5$, width 1.2 µm and separated by 300 nm of air. (d) Bulk 1 µm aSi cell absorbance in two cases. Simple anti-reflective coating (ARC): 100 nm layer of $n = 1.5$ as shown in Fig. 6(b). Complex ARC: optimized six-layer Bragg mirror and impedance matching layers above and below (see notebook). (e) Surface plasmon resonance excitation with a prism in Kretschmann configuration: Prism is made out of BK7 on top of a 55 nm layer of gold illuminated at 45.5° and $\lambda = 600\,{\rm nm}$. (f) Electric field intensity of a current source encased in a layer of width 1 µm and $\epsilon = 1.4 + 0.1i$ at $\lambda = 600\,{\rm nm}$.

Download Full Size | PDF

7. CONCLUSION

In conclusion, we have presented PyMoosh, a user-friendly Python tool for the analysis and design of photonic multilayers. PyMoosh has been developed with open science principles in mind to make it genuinely findable, accessible, interoperable and reusable. With its comprehensive set of features, numerical stability, readiness for optimization and deep learning problems, and its extensive set of examples, we hope our platform will be adopted by the community at large to build upon. Already used in teaching and research, PyMoosh offers a versatile platform for studying the rich physics of multilayered structures, particularly those incorporating metallic or metallic-like layers. In a subsequent publication, we aim to demonstrate how PyMoosh can be seamlessly integrated with advanced optimization algorithms and optimization libraries and that its features allow for an efficient exploration of inverse design problems [63,65]. Future plans include incorporating more advanced material descriptions, such as spatial dispersion in metals [66,67], to further expand the capabilities of PyMoosh. We are dedicated to continuously improving PyMoosh and welcome any suggestions and contributions from the community.