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

The vacuum wavelength $\lambda _0$ and refractive index $n(z)$ are related to $\omega$, $\varepsilon (z)$ and $\mu$ by

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

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

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

Substituting 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

The general solution of Eq. (33) can be given in terms of parabolic cylinder functions $D_\nu$ as a linear combination of

The 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.

 

Fig. 1. Simulated DBR with graded layers.

Download Full Size | PDF

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:

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.

 

Fig. 2. Error in power reflectance for the simulated DBR using different approximations. In the top figure, signed error is marked on a linear scale. On the bottom, its absolute value is shown on a logarithmic scale. In the top figure, the curves for SC$_2$ and SC$_5$ are shown only near 980 nm for better readability.

Download Full Size | PDF

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.

 

Fig. 3. Error in the determination of the phase of the reflected wave for normal incidence from the substrate.

Download Full Size | PDF

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:

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.

 

Fig. 4. Error of reflection amplitude calculations for non-perpendicular incidence and TE (or s) polarization. Top: spectra for approximations GP and SC$_5$, for $k_{\mathrm {t}}$ equal to $0$, $4\sigma _k$ and $k_{\textrm {TIR}}$. Bottom: Error at 980 nm as a function of $k_{\mathrm {t}}$. In the figure on the top, the continuous lines represent the same case as the corresponding lines in Figs. 2 and 3. Number $k_{\textrm {TIR}}$ is calculated for $\lambda _0= 980\,$nm.

Download Full Size | PDF

 

Fig. 5. Error of reflection amplitude calculations at 980 nm for non-perpendicular incidence and TM (or p) polarization as a function of $k_{\mathrm {t}}$.

Download Full Size | PDF

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.

 

Fig. 6. Error in the transmission amplitude for three different thicknesses of the grating layer. Solid lines denote the whole layer approximated by two layers (i.e. GP approximation). The dashed lines show the situation when the gradient is first cut into two gradient sub-layers, which are then GP-approximated (each one separately). The dashed-dot line is similar to the latter, but the 100-nm gradient is first divided into three sub-gradients.

Download Full Size | PDF

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.

 

Fig. 7. Error in transmission amplitude for different values of absorption on the low-refractive-index edge of the analyzed 20-nm thick gradient. At the other edge $\mathrm {Im}(n_r) = 0$. The values for absorption refer to a vacuum wavelength of 980 nm.

Download Full Size | PDF

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.

 

Fig. 8. The difference in the results obtained using SC approximation and our approximation as a function of the number of SC steps ($N$) in both 2D and 3D (the curves for both cases are indistinguishable in the scale of the graph). The slopes for the straight lines in the graph are: $-0.60$ for the Q factor and $N\leq 20$; and $-0.52$ for the wavelength and $N\leq 10$.

Download Full Size | PDF

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.