## Abstract

In this paper, we present an approximation intended to find applications in numerical simulations of optical phenomena in layered structures. The method can be used to avoid approximating graded layers by using numerous homogeneous layers in the simulation. In our approach, a single layer with a graded refractive index profile or any number layers can be replaced with only two layers that for a selected wavelength and normal incidence imitate exactly the optical properties of the replaced layer or layers. The proposed approximation is valid for a wide range of wavelengths and incidence angles. It is especially useful in time-consuming simulations (especially in 3D), where it is of paramount importance to keep the number of layers in the simulated structure low.

© 2021 Optical Society of America under the terms of the OSA Open Access Publishing Agreement

## 1. Introduction

Graded-refractive-index layers are used for many purposes in a variety of optic and optoelectronic devices, including as wide-spectrum anti-reflecting coatings [1,2], in the active regions of semiconductor lasers [3,4], and in electrically conductive distributed Bragg reflectors (DBRs) for vertical-cavity surface-emitting lasers (VCSELs). In its basic form, a DBR consists of many alternating quarter-wavelength-thick layers of materials with different refractive indices. The real reflectivity of such mirrors can be close to $100\%$ (exceeding $99.9\%$), depending on the contrast between the indices in the materials and the number of pairs of layers. DBRs with such high reflectivities are used in VCSELs and femtosecond lasers, among other devices [5]. The application of graded layers deteriorates somewhat the optical properties of a DBR (reduces its reflectivity). In devices where electrical current is put through a DBR, graded interfaces are used to reduce electrical resistance. Since the 1990s, arsenide VCSELs, which are by far the most common kind of VCSEL, have usually been equipped with $\mathrm {Al}_x\mathrm {Ga}_{1-x}\mathrm {As/Al}_y\mathrm {Ga}_{1-y}\mathrm {As}$ or $\mathrm {Al}_x\mathrm {Ga}_{1-x}\mathrm {As/GaAs}$ conductive DBRs with graded interfaces [6–11]. In an arsenide VCSEL, the output-coupling DBR usually contains around 20 pairs, whereas bottom DBR contains 30 or more pairs. The whole structure contains around 100 thin graded index layers (two per DBR period). Because the optical length of the whole period is very close to $\lambda /2$ (where $\lambda$ is the wavelength of the radiation emitted by the laser), the optical thickness of the graded layers is well below $\lambda /4$.

Numerical modelling plays a very significant role in the design of VCSELs and other optoelectronic devices. The optical properties of the VCSEL structure can be modelled to determine the reflectivity of individual DBRs (this is usually the first step) and the optical modes in the whole laser cavity, which contains two DBRs as well as the active region. When modelling reflectivity, a 1D approximation can often be used. However, this is not possible when calculating the optical modes of the cavity. For such calculations, the distribution of the refractive index in the directions perpendicular to the propagation of light is fundamentally important, because it determines the optical losses of the modes. Often, the main mechanism for mode selection in the laser is the slight lateral variation in the refractive index that results from the non-homogeneous temperature distribution in the device. For this reason, full 3D (preferably vectorial) optical models should be used. In some cases, when the symmetry of the problem allows, a 2D approximation can be used.

Because numerical models based on discretization in 3 and even 2 dimensions are time- and memory consuming, many optical models treat one direction (usually perpendicular to the physical layers, and, in many cases, such as for VCSELs, at once the direction of light propagation) in a special way, in order to reduce the complexity of the equations that have to be solved numerically. However, many numerical approaches require the layers to be homogeneous in the distinguished direction. An important example are modal methods, such as RCWA [12–14] and PWAM [15], offering significant advantages [16,17] in certain important cases. When modelling layers with a refractive index graded in this distinguished direction, they can be approximated by an arbitrary number of steps (homogeneous layers). In [18], the authors analyzed in 1D the convergence of the calculated reflectivity of a DBR with graded layers (in 1D the exact solution can be obtained) with the reflectivities of DBRs where the graded refractive index profile in the layers was approximated by a step function of $N$ steps. More precisely, the graded layer was divided into $N$ parts of equal thickness, and with refractive indices equal either to the midpoint or to the boundary value. The results presented in [18] show that at least five steps are necessary to obtain a reasonably good approximation. In graded DBR VCSELs, where there can be 100 graded layers or more constituting half of all the layers, such a dense division is in most cases unacceptable for 2- or 3D optical models. The amounts of computational memory and time necessary become impractically high, and numerical errors become increasingly significant.

In this paper, we show that the optical properties of any graded layer of a sub-wavelength thickness can indeed be approximated with very good accuracy, by replacing this layer with only two homogeneous layers (with a collective thickness equal to that of the graded layer), provided the thicknesses and refractive indices of the two steps are calculated according to the procedure we have elaborated. Although our approximation is exact for one arbitrarily chosen wavelength, it is valid for a very wide range of wavelengths, which crucially allows for its utilization in practical applications where the frequency of the electromagnetic wave is not *a priori* known, such as for the determination of cavity modes.

## 2. Theory and derivation

Consider a layered optical medium with a scalar $z$-dependent permittivity $\varepsilon$ and a constant scalar permeability $\mu$. The plane wave solutions of Maxwell equations can be obtained by postulating the variable-separated form of the physical fields $\mathbf {E}$ and $\mathbf {H}$ (see Appendix). Denote by $E_{\parallel }$ and $H_{\parallel }$ any of the components of $\mathbf {E}$, and respectively $\mathbf {H}$, parallel to the layers. For TE mode ($s$ polarization, $E_z=0$) the assumption

leads to the following Helmholtz equation: for $k(z) = \omega \sqrt {\varepsilon (z) \mu }$ and constant $k_{\mathrm {t}}$. Similarly, for TM mode ($p$ polarization, $H_z=0$) from we obtainThe vacuum wavelength $\lambda _0$ and refractive index $n(z)$ are related to $\omega$, $\varepsilon (z)$ and $\mu$ by

and where $c_0$ is the vacuum speed of light. Moreover, $k_{\mathrm {t}}/k$ can be interpreted as the sine of the angle between the normal to the layers and the direction of propagation.Our goal is to use Eq. (2) for fixed $\omega$ and $k_{\mathrm {t}} =0$ (i.e. for a fixed vacuum wavelength and normal propagation) to exactly replace a single layer of variable $k$ by two homogeneous layers of some constant $k_A$ and $k_B$. In the next section, we provide numerical evidence that this replacement is of high accuracy (although it is not exact) not only for the selected $\omega$ and $k_{\mathrm {t}} =0$, but also for a wide range of wavelengths around $\lambda _0$ and various directions of propagation. Here, the term *exact replacement* can be explained in the following way. For a layer within an interval $[z_0,z_1]$ and an initial condition for Eq. (2) organized into a column $[f(z_0), f'(z_0)]^{T}$, one can express $f(z_1)$ and $f'(z_1)$ as

Notice that the above equality assumes continuity of $f$ and $f'$ at the interface between layers $A$ and $B$. This is consistent with the physical behavior of the corresponding electric field. Since $z_0 = 0$ can always be chosen and constant $k_B$ makes $\mathcal {T}_B$ translation-invariant, one can simply write

for $d=z_1-z_0$, $d_A=z_p -z_0$, $d_B=z_1-z_p$ and treat all matrices $\mathcal {T}$ as acting on initial conditions at $z=0$.Superficially, the system (9) seems to involve four equations for three unknowns $k_A$, $k_B$ and $z_p$. However, these equations are not independent, as Abel’s identity shows. Indeed, a transition matrix for the second order ODE

can be calculated according to the following formula where the fundamental matrix $\mathcal {M}(z)$ is arranged from values and derivatives of two linearly independent special solutions $g_1, g_2$ of Eq. (10)The Abel’s identity expresses the $z$-dependence of Wronskian $W(z) = \det \mathcal {M}(z)$ solely by $p(z)$

This formula, when applied to Eq. (2), shows that all corresponding Wronskians must be constant. Consequently, the determinants of all transition matrices $\mathcal {T}_L$, $\mathcal {T}_A$, $\mathcal {T}_B$ must also be constant and equal to $1$. As a result, one of the equations in system (9) can be obtained from the remaining three, at least locally. Thus, the actual number of equations and the number of unknowns match.

Let us examine system (9) in some detail. Since, for a moment, no particular form of $\varepsilon (z)$ is assumed for $\mathcal {T}_L$, let

In layers $A$ and $B$, $\varepsilon$ is constant and the following form of matrix $\mathcal {T}_A$ can be obtained analytically:

From Eqs. (16) and (19) it immediately follows that

andSubstituting relation (20) into Eqs. (18) and (17), one can solve them for cotangents and obtain

The product of the above equalities gives a formula for $\cot ( k_A d_A) \cot ( k_B d_B)$, which can be compared with formula (21). The resulting condition reads:

Thus, our system of equations for $d_A$, $d_B$, $k_A$ and $k_B$ can be rewritten as

Solving Eqs. (25), (26) and (27) for $d_B$, $k_B$ and $d_A$ in terms of $k_A$ one obtains

The above dependencies can be substituted into formula (28), producing a relation solely for $k_A$

Eventually, Eq. (32) can (in fortunate cases) be solved numerically. The resulting $k_A$ inserted into formulae (29), (30), and (31) yields the remaining unknown quantities.

So far, we have not assumed any special form of the matrix $\mathcal {T}_L (d)$ in relation (14). In the next section, we consider the case of a linear refractive index profile within graded layers. Linearity of $n(z)$ implies linearity of $k(z)$, so Eq. (2) can be rewritten for $k_{\mathrm {t}} = 0$ as

with $\beta$ and $\gamma$ determined by input and output refractive indices $n_{\mathrm {in}}$, $n_{\mathrm {out}}$, together with vacuum wavelength $\lambda _0$ and layer thickness $d$The general solution of Eq. (33) can be given in terms of parabolic cylinder functions $D_\nu$ as a linear combination of

andThe explicit form of the resulting fundamental matrix, and especially the $\mathcal {T}_L (d)$ matrix, seems rather cumbersome to present here. It can, however, be derived easily from formulae (11) and (12) using a computer algebra system. The code in which it was implemented (also for nonzero $k_{\mathrm {t}}$) is provided in the Code 1 [19].

Finally, let us remark that for normal propagation ($k_{\mathrm {t}} = 0$) physically equivalent results can be obtained from Eq. (4), which describes TM modes. This approach is nevertheless more complicated, as it involves discontinuities of $h'(z)$ for abrupt $\varepsilon (z)$ profiles (see Appendix).

## 3. Applications and analysis

In this section, we show how the approximation described above can be applied in simulations of photonic structures. We present the results of 1D, transfer-matrix-method (TMM) simulations, comparing our approximations with the exact solutions and the staircase midpoint approximation. Unless stated otherwise, we assume a linear refractive index profile in the graded layer. First, let us analyse a DBR with graded-index layers, designed for 980 nm. The simulated DBR has 15 periods, each consisting of an AlAs and a GaAs layer separated by 20-nm-thick gradient layers. The DBR is deposited on an infinitely thick GaAs substrate. The incident (and reflected) light propagates in the GaAs substrate. A scheme is presented in Fig. 1. The refractive index change in the graded layers is the highest for the AlGaAs structures, so this is the most challenging case. We assume 3.516 and 2.939 as the refractive indices for GaAs [20] and AlAs [21], respectively (corresponding to a wavelength of 980 nm). We assume that the refractive index is wavelength independent, in order to show that, although exact only at one wavelength (980 nm), our approximation is very accurate for a very wide range of wavelengths. As a consequence, especially when strong absorption appears in GaAs at $\lambda \lesssim 870\,$nm, the results do not describe a real mirror.

The reflectivity of a DBR is usually very close to 100% (for the wavelength it is designed for), and differences in the order of 0.01% can have a significant impact on the performance of devices such as VCSELs. For this reason, optical models of DBRs need to be very accurate. We compare the following approximations:

- GP where each gradient is replaced with a pair of homogeneous layers using the parameters obtained from our model;
- HP where each DBR half period is replaced with a pair of homogeneous layers using the parameters obtained from our model;
- SC
_{N}where each gradient is replaced by a staircase approximation with $N$ steps of equal widths.

Each DBR period then contains 6, 4, $2N+2$ homogeneous layers, for the SG, HP, SC$_N$ approximations, respectively. The wavelength for which the GP and HP structures were calculated is 980 nm. In the GP structure, the gradients are divided into layers of thicknesses $d_1\approx 10.18\,$nm and $d_2\approx 9.82\,$nm, with the corresponding refractive indices $n_1=3.038$ and $n_2=3.421$. For the approximation SC$_2$, which also uses two layers, the refractive indices are $3.083$ and $3.372$, so the two approximations are significantly different. In the case of the HP structure, the physical interface between the gradient and the adjacent homogeneous layer does not overlap with the interface of the approximating layers, but the total thickness of the two layers is the same. Figure 2 compares the accuracy of the power reflectance calculations for the DBR using the considered approximations. The incident light falls perpendicularly on the mirror surface. We calculate the error relative to a TMM solution, which is exact (for a plane wave) within the accuracy of floating-point arithmetic. The GP approximation provides excellent accuracy across the whole, very wide, spectral range considered. To achieve similar accuracy requires as many as around 50 steps in each gradient layer using SC approximation, which is often unacceptable in real simulations. A more realistic SC$_5$ gives a generally non-negligible error of around 0.01% at 980 nm (although it still generates many additional layers) and much higher error for the wavelengths in the DBR stop-band.

In addition to power reflectance, another equally important parameter is the phase of the reflected wave. Figure 3 shows the error of the calculated phase. As a benchmark to assess which phase changes matter, the following example can be used. A 980 nm vacuum-wavelength wave in GaAs traveling a distance of 1 nm changes its phase by roughly $1.8\cdot 10^{-3}$. A ten times lower error in the calculations of the phase of a wave reflected by a DBR would modify the cavity vacuum wavelength by about 0.7 nm, for a $\lambda /2$-thick GaAs cavity. This can be a significant difference in some applications, so an error higher than $10^{-4}$ may be considered significant.

So far, we have considered the normal incidence of a plane wave, for which our approximation gives an exact result for the selected wavelength. However, when a monochromatic non-plane wave (for instance, with a Gaussian lateral profile), or different incidence angles are considered, the non-zero lateral components $k_{\mathrm {t}}$ of the wavevectors become important. To determine the most important range of values for $k_{\mathrm {t}}$, let us analyse the case of a Gaussian mode in the laser cavity. Its diameter is determined by various parameters, which can be controlled to a certain degree. However, the diameter is limited from below by the wavelength. We can assume that in a 980 nm, arsenide VCSEL a Gaussian mode with a standard deviation of $\sigma = 1\,$µm is an extremely narrow mode, requiring the widest range of $k_{\mathrm {t}}$. The Fourier transform of such a profile is also a Gaussian distribution, with the standard deviation $\sigma _k = 1/\sigma$. Almost the entire energy of the wave (over 99%) is related to the components of their $k_{\mathrm {t}}$ in the range $[0, 3\sigma _k]$.

Another characteristic value for $k_{\mathrm {t}}$ that can be used is the total-internal-reflection limit for GaAs–vacuum transitions. This value is given by the following formula:

where $\lambda _0$ is the vacuum wavelength (980 nm in the cases considered here). In Fig. 4, we compare the accuracy of the considered approximations for different values of $k_{\mathrm {t}}$. Instead of analysing power reflectance and phase separately, as in the previous plots, we calculated the error of determination for the complex reflection amplitude $r$ for the electric field, which represents both the power reflectance $R$ and the phase $\varphi$ as follows:We present only one SC structure, SC$_5$, because the approximation SC$_2$ is too inaccurate and SC$_{50}$ is impractical for many applications. As may be expected, the wavelength for which approximations GP and HP are exact shifts towards shorter wavelengths with increasing $k_{\mathrm {t}}$, so that the vertical component remains close to that of a 980 nm wave at normal incidence. We consider both possible polarizations: TE (or s) where the electric field of the wave is parallel to the layers, and TM (or p) where the magnetic field of the wave is parallel to the layers. When the incidence angles are smaller than the TIR angle ($k_{\mathrm {t}}<k_{\textrm {TIR}}$), both GP and HP approximations have small errors across the whole considered wide spectral range (HP approximations are not shown on the spectrum in Fig. 4 for the sake of readability). However, as was seen in the previous graphs, GP is significantly more accurate than HP. In the graph on the bottom in Fig. 4 and in Fig. 5, the TIR limit is clearly visible. The error is on average the highest in the region above TIR, with the exception of SC$_5$ for TM polarization. For all the approximations in the TIR regime, the power reflectance $|r|^{2}=1$, so the error is in the calculation of the phase of the reflected wave. Our approximations are more accurate for TE polarization; however, for both polarizations and angles of incidence below the TIR limit the error is small.

So far, we have proven that our approximations work for a thin (20 nm) layer. In fact, the thinner the layer that is replaced with two substitute layers, the more accurate the approximation. This is why the GP approximation works better than the HP approximation. To investigate this effect, we calculated the transmission amplitude $t$ through a single gradient layer, for three different thicknesses: 20 nm, 50 nm, and 130 nm. The light is transmitted from the high-refractive-index side (GaAs) to the low-refractive-index side. The gradient connects the homogeneous layers (GaAs and AlAs), so the refractive index is continuous. This time, instead of $r$ we calculate $t$, because for such a structure reflectance is close to 0, while $|t|\approx 1$. When the thickness of the approximated gradient (or more generally non-homogeneous) layer increases, so does the inaccuracy of the approximations. If the basic division into two layers does not provide sufficient accuracy, the non-homogeneous layer should first be divided into two (or if necessary more) sub-layers, and then each of the resulting sub-layers should be approximated using the method described in this paper. In Fig. 6, we show the accuracy of the approximations for gradient layers with thicknesses of 20 nm, 50 nm, and 130 nm. The red line represents the approximation used previously for analysis of the DBR, so it can be considered as a reference for a very accurate approximation. To achieve similar (and in fact improved) accuracy, the 50-nm layer should first be divided into two sub-layers (so the whole layer will be approximated finally by four homogeneous layers). In this case, each of the approximated sub-layers is 25-nm thick, so the whole layer is thicker than the layer represented by the red line. However, the refractive-index contrast in the sub-layers is lower than the contrast in the single 20-nm layer, and as a result the overall inaccuracy is higher in the case of the thinner layer. In the case of the 130-nm thick layer, it is enough to have 6 layers in total (3 gradient sub-layers, each then GP-divided) to achieve similar (in fact better) accuracy than in the case of the 20-nm layer.

Because the number of layers necessary to obtain good accuracy seems to have a logarithmic dependence on the thickness of the gradient layer, even very thick gradients can be approximated by a limited number of homogeneous layers.

#### 3.1 Absorption

So far, we have considered non-absorbing (or, more generally, non-homogeneous) grating layers, which means that the refractive indices are real-valued. Formally, there is no problem with having one of each of the peripheral refractive indices in a grating layer complex-valued. The equations for finding the parameters of the two approximating homogeneous layers remain valid and can be solved. However, in general the resulting values for the two thicknesses and refractive indices will be complex. A complex value for thickness has no physical sense, although it can be perfectly acceptable in some computational methods, such as TMM. In general, however, we have to use its real part as the thickness of the approximating layer. Unless metallic layers are considered, the imaginary part of the refractive index of a layer is a few orders of magnitude smaller than the real part. The neglected imaginary part of the calculated thickness is then very low and this modification introduces only a small error. Complex values for the refractive index are perfectly valid and can be used directly when approximating the layers.

In the worst-case scenario, the absorption is very different at both ends of the gradient, so we can assume that it is 0 at one end with high absorption at the low-refractive-index end. In order to show the impact of absorption and the projection of the calculated thicknesses onto the real axis, we calculate the error of the GP approximation for a 20-nm gradient. Except for the absorption, this case is identical to the gradient denoted by the red line in Fig. 6. Figure 7 shows the approximation error for different values of absorption on the AlAs-edge, when the absorption on the other edge is kept 0. The values for absorption refer to the vacuum-wavelength $\lambda _0 = 980$ nm. The imaginary part of the refractive index $\mathrm {Im}(n_r)$ is kept constant, so the actual absorption $\alpha$ depends on the wavelength according to the formula

The considered values do not refer to a real case; they simply create the worst-case absorption distribution. Absorption of $10^{4}\,$1/cm is extremely high and is roughly equal to the inter-band absorption in GaAs, so it can be considered the upper limit for the class of materials considered here. In this case, the error is relatively high, close to $10^{-4}$. Whether it is too high to be acceptable depends on the structure containing a similar gradient that is to be simulated. Generally, layers with such high absorption are avoided in optical devices, or are used in regions where the optical field is very weak, reducing its usually negative impact. However, for some applications, such as in photodetectors, it is necessary to approximate layers with high absorption with good accuracy. Accuracy can be improved by sub-dividing the gradient as described earlier. Each case should be considered separately, because we present an artificial, worst-case situation. It is important to note that our approximation also works with layers with high absorption across a very wide spectral range.

#### 3.2 Multidimensional analysis

So far, we analyzed our approximation and the staircase approximation in 1D case, where the exact solution can be found. In a structure, which is not homogeneous in lateral directions, exact solutions are generally not available and the optical field is not a plane wave. Above, we showed the validity of our aproximation for different lateral components of the plane wave vector, which implies its validity in the case of anis arbitrary lateral distribution of the optical field. Here, we directly compare the results of simulations of a laterally non-homogeneous structure where our approximation and the staircase approximation was used. The simulated structure is a simplified version of a typical VCSEL cavity designed for emission at $980$ nm. It consists of a 30-period bottom DBR, a one-$\lambda$ GaAs cavity and a 15-period top DBR. The DBR structure is the same as in the previous simulations (see Fig. 1) and 0 absorption in the whole structure is assumed. The width of the cavity is 6 µm and the cavity is surrounded by air, except from an infinitely thick and wide GaAs substrate. The finite width of the cavity introduces a lateral non-homogeneity and results in appearance of transverse modes. We perform the simulations in 2D and 3D, using PWAM model, designed for such simulations [15]. In 2D the cavity has a rectangular shape, while in 3D is in form of a cuboid. We calculated the fundamental (approximately Gaussian) mode wavelength and the quality factor of the cavity in several different approximations of the gradient layer. As a reference point we use the wavelength ($\lambda ^{\mathrm {2D}}_0\approx 980.63$ nm, $\lambda ^{\mathrm {3D}}_0 \approx 980.32$) and the Q factor ($Q^{\mathrm {2D}}_0=13861$, $Q^{\mathrm {3D}}_0 = 13934$) of the cavity where the gradient layers were approximated using our approximation. Because the results obtained using the staircase approximation seem to converge towards those values, we refer to the respective differences as errors of the staircase approximation. In Fig. 8 we show the relative difference in the cavity Q factor and the wavelength of the mode:

Surprisingly, they are practically identical (except a small difference in wavelength in the case of SC$_1$) in 2D and 3D, so the graph refers to both cases simultaneously (with the 2D values at $N=1$). In all the analysed cases, the staircase approximation resulted in a shorter wavelength and a lower Q factor. The difference in wavelengths is very small, rather negligible in many applications. The cavity mode wavelength is mainly defined by the optical path length in the cavity material between the mirrors. In our example, for the sake of simplicity (in order to avoid having gradients with different thicknesses), we resigned from thick gradient layers usually used in VCSELs between the active region (at the center) and the mirrors. For this reason, in our example, the type of approximation of the gradient layers has only little impact on the wavelength, which might not be the case in more complex (and realistic) structures. However, a convergence of the staircase approximation towards our approximation is clearly visible at least up to the 20-step division. In the cavity Q factor, however, the difference is much more significant. The cavity Q factor determines the laser threshold gain and hence threshold current. The error in the Q factor is significant up to at least $N=5$, where it drops to around $1\%$, while for $N=1$ it exceeds $20\%$, a huge error which is unacceptable in the cases where quantitatively reliable results are sought after. What is, however, most important is the fact that the staircase approximation converges slowly (roughly as $N^{-0.6}$ in the presented case) to our 2-layer approximation. This proves that our approximation outdoes the staircase approximation also in multidimensional calculations.

## 4. Dispersion and similar phenomena

So far, we have only considered nondispersive media. In real situations, the refractive indices of the materials used in a device do depend on the wavelength and temperature. Both phenomena often are crucially important, although the induced variations in the indices are very small. Of course, the parameters of the approximating layers can be calculated separately for each wavelength or temperature, but if the thicknesses of the layers depend on these parameters it may be difficult to apply the approximation in a simulated structure. Fortunately, when real-life values for dispersion or thermal dependence were assumed, the thicknesses of the two layers were almost identical (and consequently approximately equal to $d/2$ where $d$ is the thickness of the approximated layer), in all the cases we analysed. Moreover, if the derivative $dn_r/d\xi$, where $\xi$ is the considered parameter (for instance $\xi =\lambda _0$ or $\xi = T$), is the same at both ends of the gradient, it seems that the same value of the derivative applies to the approximating layers. In any case, the parameters of the layers can be plotted using the Python module provided in the Code 1 [19] (in the 1D case), as functions of $\xi$. The appropriate derivatives for the approximating layers can be fitted to the results. In the simulations of a structure, a linear approximation of dispersion or thermal dependence should be perfectly sufficient, because the range of wavelengths which need to be taken into account is not very wide. Thermal changes usually remain linear for temperatures changing by several hundreds of K. This may be not true if the considered photon energies are close to, for instance, the band gap of a material in the simulated structure.

## 5. Summary

We have presented a method for approximating the optical properties of gradient layers that is far more accurate and effective than the widely used staircase approximation. We used 1D simulations, for which exact solutions were available, to analyse the accuracy of our method in various situations that appear in real situations. The presented validity of our approach in a wide range of directions of light propagation shows that the approximation remains applicable in multidimensional simulations. In 2D and 3D we compared our approximation with the staircase approximation, and we showed a slow convergence of the latter, which further confirms superiority of our approach. Even when absorption was taken into account, the approximations proved to be very accurate across a wide spectrum and incidence angles, if the thicknesses of the approximated layer or layers were not too high. The real power of our approximations applies to 2D and especially 3D simulations, where the number of layers used to describe the simulated structure has to be kept as low as possible, in order to keep calculation time and memory consumption within available limits.

The proposed procedure enables the replacement not only of a gradient layer but also of a group of different layers with only two layers. This can be helpful when the number of layers in an optical simulation has to be reduced for some reason. However, one should bear in mind that the thicker the approximated region, the less reliable the approximation becomes when the wavelength varies.

In this paper, we analysed materials with parameters similar to those of arsenide semiconductors. However, the software we provide in the Code 1 [19] can be used in principle with any sensible parameters to calculate the parameters of the approximated layers, according to the needs of the analysis.

## Appendix: Helmholtz equations in 1D

To keep our presentation self-contained, let us derive the Eqs. (2) and (4) from basic principles. We start with Faraday’s and Ampére’s laws:

where $\mathbf {D}(\mathbf {x},t)=\varepsilon (\mathbf {x}) \mathbf {E}(\mathbf {x},t)$ and $\mathbf {H}(\mathbf {x},t)=\frac {1}{\mu (\mathbf {x})} \mathbf {B}(\mathbf {x},t)$, for scalar but position-dependent $\varepsilon$ and $\mu$. Applying $\nabla \times$ to both equations and mutually eliminating the variables one obtainsGauss’s laws

together with the identity $\nabla \times (\nabla \times \mathbf {A}) = \nabla (\nabla \cdot \mathbf {A}) - \nabla ^{2} \mathbf {A}$, where $\nabla ^{2}$ is the vector Laplacian produces the following, fairly general, inhomogeneous form of the wave equationsEven for non-vanishing $\rho$ and $\mathbf {j}$, it is important to study the homogeneous variant of these relations. By postulating a variable-separated form of the solutions

the homogeneous wave equations become Helmholtz equations for $\mathbf {E}(\mathbf {x})$ and $\mathbf {H}(\mathbf {x})$Let us specialize now to a medium layered along the $z$ axis (hence $\varepsilon (\mathbf {x})=\varepsilon (z)$) with constant $\mu$. For the TE case ($E_z=0$) the system of equations given by formula (52) decouples, producing

with $E_i$ standing for $E_x$ or $E_y$. Let $E_{\parallel }$ be either of these components. By postulating the separation of variables again the standard form of Helmholtz equation for the $xy$-dependence of $E_{\parallel } (\mathbf {x})$ can be obtainedSimilarly, for TM modes ($H_z=0$) the system (53) becomes

Finally, let us remark on continuity conditions for $f(z)$ and $h(z)$. Suppose there is an abrupt change in $\varepsilon$ at $z_0$. Denote by $\varepsilon (z_{0+})$ and $\varepsilon (z_{0-})$ the right- and left-hand limits of $\varepsilon (z)$ at $z_0$, and similarly for limits of $f$, $f'$, $h$, and $h'$. From the integral representation of Maxwell equations it can be derived that tangential ($x$ and $y$) components of $\mathbf {E}$ and $\mathbf {H}$ must be continuous across the interface. For the TE mode, this means that

Then, substitution of relations (50) and (51) into the system (42) yields

with $\eta '=x$ for $\eta =y$ and vice versa. From the relation (63) it follows that for globally constant $\mu$ the derivative $\frac {\partial E_\eta }{\partial z}$ must be continuous, henceSimilarly, for TM mode we have

and the system (43) produces showing that $\frac {1}{\varepsilon }\frac {\partial H_\eta }{\partial z}$ must be continuous. In turn, and this relation must be taken into account when gluing together solutions of Eq. (4) for two constant but distinct values of $\varepsilon$.## Funding

Narodowe Centrum Badań i Rozwoju (MAZOWSZE/0032/19-00).

## Acknowledgements

This research was partially supported by the National Centre for Research and Development (grant MAZOWSZE/0032/19-00).

## Disclosures

The authors declare no conflicts of interest.

## Data availability

Data underlying the results presented in this paper are not publicly available at this time but may be obtained from the authors upon reasonable request.

## References

**1. **M.-L. Kuo, D. J. Poxson, Y. S. Kim, F. W. Mont, J. K. Kim, E. F. Schubert, and S.-Y. Lin, “Realization of a near-perfect antireflection coating for silicon solar energy utilization,” Opt. Lett. **33**(21), 2527–2529 (2008). [CrossRef]

**2. **W. Lowdermilk and D. Milam, “Graded-index antireflection surfaces for high-power laser applications,” Appl. Phys. Lett. **36**(11), 891–893 (1980). [CrossRef]

**3. **D. Fekete, M. Yasin, A. Rudra, and E. Kapon, “Very low transparency currents in double quantum well ingaas semiconductor lasers with *δ*-doped resonant tunneling,” Appl. Phys. Lett. **92**(2), 021109 (2008). [CrossRef]

**4. **S. Stańczyk, T. Czyszanowski, A. Kafar, J. Goss, S. Grzanka, E. Grzanka, R. Czernecki, A. Bojarska, G. Targowski, M. Leszczyński, T. Suski, R. Kucharski, and P. Perlin, “Graded-index separate confinement heterostructure ingan laser diodes,” Appl. Phys. Lett. **103**(26), 261107 (2013). [CrossRef]

**5. **A. Jasik, P. Wasylczyk, P. Wnuk, M. Dems, A. Wójcik-Jedlińska, K. Regiński, L. Zinkiewicz, and K. Hejduk, “Tunable semiconductor double-chirped mirror with high negative dispersion,” IEEE Photonics Technol. Lett. **26**(1), 14–17 (2014). [CrossRef]

**6. **P. Zhou, J. Cheng, C. F. Schaus, S. Z. Sun, K. Zheng, E. Armour, C. Hains, W. Hsin, D. R. Myers, and G. A. Vawter, “Low series resistance high-efficiency GaAs/AlGaAs vertical-cavity surface-emitting lasers with continuously graded mirrors grown by MOCVD,” IEEE Photonics Technol. Lett. **3**(7), 591–593 (1991). [CrossRef]

**7. **M. Sugimoto, H. Kosaka, K. Kurihara, I. Ogura, T. Numai, and K. Kasahara, “Very low threshold current density in vertical-cavity surface-emitting laser diodes with periodically doped distributed Bragg reflectors,” Electron. Lett. **28**(4), 385–387 (1992). [CrossRef]

**8. **Y. Xiong, Y. Zhou, Z. Zhu, Y. Lo, C. Ji, S. Bashar, A. Allerman, T. Hargett, R. Sieg, and K. Choquette, “Oxide-defined GaAs vertical-cavity surface-emitting lasers on Si substrates,” IEEE Photonics Technol. Lett. **12**(2), 110–112 (2000). [CrossRef]

**9. **M. Geębski, J. A. Lott, and T. Czyszanowski, “Electrically injected VCSEL with a composite DBR and MHCG reflector,” Opt. Express **27**(5), 7139–7146 (2019). [CrossRef]

**10. **V. P. Kalosha, V. A. Shchukin, N. Ledentsov, and N. N. Ledentsov, “Comprehensive analysis of electric properties of oxide-confined vertical-cavity surface-emitting lasers,” IEEE J. Sel. Top. Quantum Electron. **25**(6), 1–9 (2019). [CrossRef]

**11. **R. Michalzik, * VCSELs* (Springer, 2013), vol. 166 of

*Springer Series in Optical Sciences*, chap. VCSEL Fundamentals, pp. 19–75.

**12. **M. Moharam and T. K. Gaylord, “Rigorous coupled-wave analysis of planar-grating diffraction,” J. Opt. Soc. Am. **71**(7), 811–818 (1981). [CrossRef]

**13. **M. Neviere and E. Popov, * Light propagation in periodic media: differential theory and design* (CRC Press, 2002).

**14. **P. Lalanne and E. Silberstein, “Fourier-modal methods applied to waveguide computational problems,” Opt. Lett. **25**(15), 1092–1094 (2000). [CrossRef]

**15. **M. Dems, R. Kotynski, and K. Panajotov, “Planewave admittance method–a novel approach for determining the electromagnetic modes in photonic structures,” Opt. Express **13**(9), 3196–3207 (2005). [CrossRef]

**16. **M. E. Solano, M. Faryad, A. Lakhtakia, and P. B. Monk, “Comparison of rigorous coupled-wave approach and finite element method for photovoltaic devices with periodically corrugated metallic backreflector,” J. Opt. Soc. Am. A **31**(10), 2275–2284 (2014). [CrossRef]

**17. **G. Yoon and J. Rho, “MAXIM: Metasurfaces-oriented electromagnetic wave simulation software with intuitive graphical user interfaces,” Comput. Phys. Commun. **264**, 107846 (2021). [CrossRef]

**18. **P. Debernardi and R. Orta, “Analytical electromagnetic solution for Bragg mirrors with graded interfaces and guidelines for enhanced reflectivity,” IEEE J. Quantum Electron. **43**(3), 269–274 (2007). [CrossRef]

**19. **M. Dobrski and M. Wasiak, “layers.py: Python module for approximating gradient layers and systems of layers by a pair of homogeneous layers,” figshare (2021) https://doi.org/10.6084/m9.figshare.14977086.

**20. **T. Skauli, P. S. Kuo, K. Vodopyanov, T. Pinguet, O. Levi, L. Eyres, J. Harris, M. Fejer, B. Gerard, L. Becouarn, and E. Lallier, “Improved dispersion relations for GaAs and applications to nonlinear optics,” J. Appl. Phys. **94**(10), 6447–6455 (2003). [CrossRef]

**21. **S. Gehrsitz, F. K. Reinhart, C. Gourgon, N. Herres, A. Vonlanthen, and H. Sigg, “The refractive index of Al_{x}Ga_{1} − *x*As below the band gap: Accurate determination and empirical modeling,” J. Appl. Phys. **87**(11), 7825–7837 (2000). [CrossRef]