## Abstract

In this work, black silicon (BSi) structures including nanocones and nanowires are modeled using effective medium theory (EMT), where each structure is assumed to be a multilayer structure of varying effective index, and its optical scattering in the infrared range is studied in terms of its total reflectance, transmittance and absorptance using the transfer matrix method (TMM). The different mechanisms of the intrinsic absorption of silicon are taken into account, which translates into proper modeling of its complex refraction index, depending on several parameters including the doping level. The model validity is studied by comparing the results with the rigorous coupled wave analysis and is found to be in good agreement. The effect of the aspect ratio, the spacing between the structure features and the structure disordered nature are all considered. Moreover, the results of the proposed model are compared with reflectance measurements of a fabricated BSi sample, in addition to other measurements reported in the literature for Silicon Nanowires (SiNWs). The TMM along with EMT proves to be a convenient method for modeling BSi due to its simple implementation and computational speed.

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

## 1. Introduction

Silicon is one of the most widely used semiconductors in a variety of applications; however the rather high reflectivity of polished silicon substrates (ranging from 30 – 40% in the visible range) reduces its efficiency for applications such as photovoltaic devices and sensitive photodetectors. To suppress reflection, anti-reflection coatings can be used, but are only effective over a small spectral range. On the other hand, surface micro/nano texturing of silicon proves to be an alternative technique that is efficient over a wide spectral range and more cost effective [1-3]. Structured silicon consists of a thin layer of nano-/micro-scale features on its top surface. These surface features suppress reflection and enhance the absorption of the incident optical radiation, which is of interest for a variety of applications such as solar cells [2], photodetectors [4], and black body radiation sources [5,6]. Black body radiation sources are key components for spectroscopic and gas sensing applications, and they require high emissivity materials for adequate performance. This can be achieved using structured silicon, which also has the advantage of being easily integrated with other components on the same chip [5,7]. Other applications include surface-enhanced Raman scattering [8] and even anti-bacterial surfaces [9]. Structured silicon can have the shape of pseudo-periodic conical spikes [3,10] or nano-wires [11] depending on the fabrication technique, and is called in the literature with different names depending on its surface features such as black silicon (BSi) [2,3], silicon nano-wires (SiNWs) [11], or silicon nano-rods [12]. Given in Fig. 1 illustrations showing some of these structures fabricated using different techniques [3,11,13].

Among different fabrication techniques, BSi can be produced in a maskless and convenient and fast manner at wafer scale using cryogenic plasma processing under Sulphur hexafluoride (SF_{6}) and Oxygen (O_{2}) with specific conditions of temperature, pressure, bias voltage and the plasma flow rate [3]. Within a few minutes of plasma processing BSi is formed, and it has relatively small surface features, where the height of the formed cones can be in the range of 0.5-3 μm with periodicity in the range 0.25-1 μm [3]. Another technique uses femto-second pulsed laser in the presence of SF_{6} or Chlorine (Cl_{2}), which can yield structures of relatively larger periodicities exceeding 8 μm and heights in the range of 10-40 μm [13]. Other structures include SiNWs, which can be fabricated by wet chemical etching of Si substrates in the presence of a thin layer of silver for metal-induced etching leading to highly uniform arrays of SiNWs [11].

Highly doped silicon can be used to fabricate BSi for further enhancements in the absorptance especially at longer wavelengths beyond the silicon absorption edge, where the BSi layer reduces the reflection of incident radiation, and the highly doped substrate residing below absorbs all the power leading to zero transmission [14]. Also, the high doping of silicon has a significant impact on its complex refractive index, where increasing the silicon doping causes the real part of the refractive index to decrease. This leads to lower reflection at the air/silicon interface and also causes the extinction coefficient to increase due to higher free carrier absorption [15].

The increasing number of applications and wide use of BSi along with highly doped silicon requires efficient modeling of these structures to determine the required dimensions that can give the desired optical properties. There are rigorous techniques that can be used to model BSi by considering it to be a periodic structure. These techniques include rigorous coupled wave analysis (RCWA), finite-difference time domain (FDTD) and finite element method (FEM) [16], but they require high computational power and are not suitable for design optimization. Hence a computationally fast and easy to implement method is required for modeling these structures with some approximations and limitations that still suit many applications.

In this work, a method for modeling BSi based on the Transfer Matrix Method (TMM) and the Effective Medium Theory (EMT) is proposed, compared to rigorous methods and experimental data. In addition, a thorough parametric study of the structures is carried out.

There are some modeling approaches for BSi based on EMT that are presented in the literature [17,18]. These models however treat BSi as one layer of a certain effective index, meaning there’s an abrupt change in the refractive index at the air/BSi interface and at the BSi/Si-substrate interface. This doesn’t resemble the actual gradual variation of the BSi effective index and can increase the simulated reflectance. This is taken into account in our proposed model by treating BSi as a stack of a large number of layers, where each layer has an effective complex refractive index. Also, the proposed model provides a thorough approach in modeling the complex refractive index of Si, achieved by obtaining its real and imaginary parts as a function of the wavelength, temperature and doping concentration, rendering it suitable for a wide range of applications.

This article is organized as follows: model formulation is discussed in details in section 2, where special attention is paid to taking into account different mechanisms of the intrinsic absorption in Silicon, which translates into proper modeling of the complex refraction index of Silicon depending on several parameters, especially the doping level. Then, the proposed method is verified by comparing its results with the RCWA method in section 3 identifying the model limitations. In section 4 a parametric analysis is performed for the effect of the dimensions of BSi on its optical properties. In section 5 the model is used to simulate a fabricated sample of BSi (having nanocones surface features) and the results are compared with practical measurements. The model results for some SiNWs structures are compared with those reported in the literature as well. Finally, the work is concluded in section 6.

## 2. Proposed Model Formulation

This section discusses the approach used to model the different types of micro/nano structured silicon – including BSi – using the TMM and the EMT. First the geometric topography of the modeled structures is discussed in section 2.1 and the TMM is briefly described. Then the details of obtaining the complex refractive index of silicon are discussed in section 2.2. Then the method for obtaining the complex effective permittivity is discussed in section 2.3.

#### 2.1 Modeled Surface Topography

Micro/nano structured silicon has a thin top layer of nano/micro-scale features which can be considered a layer of a certain effective index profile that depends on the shape of these features. For example, in the case of highly-ordered SiNWs the refractive index profile is expected to consist of three values (refractive index of air, effective index of SiNWs layer and the refractive index of silicon), while for the case of nanocones the effective index is expected to increase with the depth into the BSi layer as shown in Fig. 2, which is discussed in further details in the next sub-section.

To model a BSi structure its surface is assumed to consist of cones and these cones are divided into many fine stacked layers. Each layer is assumed to be a homogenous medium of a certain effective index in an approach analogous to the one implemented in references [19] and [20]. In each layer, the circular truncated cones are instead assumed to be cylinders of average radius for simplicity, which is a valid assumption when the number of layers is large and hence each layer is thin enough. An example for the assumed cones is demonstrated in Fig. 3. To compute the total reflectance and transmittance for these multi-layer structures, and hence absorptance or emissivity, the TMM was chosen [21]. Considering the assembly in Fig. 4, the electric and magnetic fields can be related at any two interfaces by multiplying by the appropriate matrices for the layers in between these interfaces.

For the entire multi-layer structure with M layers, the electric and magnetic fields in the incident and emergent mediums can be related as follows [21]:

_{r}is the medium admittance:

_{o}is the free space permittivity and μ

_{o}is the free space permeability, ${\delta}_{r}={n}_{r}k{d}_{r}\mathrm{cos}({\theta}_{r})$,

*n*is the layer complex refractive index,

_{r}*k*is the magnitude of the wave vector,

*d*is the layer thickness, and

*η*is the admittance in the emergent medium. The coefficients B and C can then be used to determine the reflectance and transmittance of the entire multilayer structure [21]:

_{N}#### 2.2 Obtaining the Complex Refractive Index of Silicon

### 2.2.1 Real Refractive Index of silicon

First the real part of the refractive index is obtained as a function of the operating wavelength *λ* and temperature *T* as follows [22]:

The above relations for the real part of the refractive index were found valid for the wavelength range of (0.5 – 10 μm) and for temperatures up to 800 K according to reference [23]. The real part of the refractive index of silicon is dependent on the doping concentration [15], and the resulting change in the refractive index is taken into consideration as follows [24]:

where ${n}_{d}$ is the doped silicon refractive index, ${n}_{u}$ is the undoped silicon refractive index,*N*is the free-carrier concentration,

*e*is the electron charge,

*ω*is the angular frequency of incident radiation, ${\u03f5}_{o}$ is the free space permittivity and

*m*is the electron mass.

### 2.2.2 Extinction Coefficient of silicon

The imaginary part of the refractive index (i.e. the extinction coefficient, denoted κ) is obtained from the total absorption coefficient α as follows [23]:

The total absorption coefficient α comprises three types of absorption: band-gap absorption α_{BG}, free carrier absorption α_{FC} and lattice vibrations α_{LV}. First, the band-gap absorption α_{BG} is found as a function of the photon energy and temperature as follows [23,25]:

*i*refers to the i

^{th}phonon with

*l = 1*refers to phonon emission and

*l = 2*refers to phonon absorption, with:

*E*is in electron volts, K

_{B}is Boltzmann constant and K

_{B}θ

_{i}is referring to the phonon energy of the i

^{th}phonon and given by 212 K

_{B}, 670 K

_{B}, 1050 K

_{B}and 1420 K

_{B}[26]. In Eq. (18) the energy gap is a function of temperature as follows [25]:where T is in Kelvin and ${E}_{g}^{0}=1.155eV,\text{A}=4.73\text{*}{10}^{-4}\text{and}\text{B}=635\text{K}$. Second, the free-carrier absorption α

_{FC}is considered using the following equations; where an electron in the conduction band can absorb an incident photon and move to a higher level in the same band [27]:where

*n*and

*p*are the total electron and hole concentrations respectively, σ

_{n}and σ

_{p}are the electron and hole absorption cross-sections, respectively, and are given by [27]: with λ in μm and T in Kelvins.

Another approach to calculate free-carrier absorption was introduced in reference [28], which yields better results for lightly doped silicon:

*λ*in μm and

*T*in Kelvins. Third, the absorption due to lattice vibrations is calculated as follows [29]:with λ in μm and T in Kelvins.

According to reference [29], Eq. (24) was derived to estimate lattice absorption in the wavelength range of (6-15 μm) because it is negligible at shorter wavelengths.

In general, for wavelengths below the absorption edge of silicon the strongest absorption scheme is band-gap absorption. Beyond the absorption edge, band-gap absorption diminishes while free-carrier and lattice vibration absorption can have different contribution depending on the temperature and doping concentration. The different absorption coefficients are compared in Fig. 5, with the free-carrier absorption calculated using Eqs. (20)–(22).

#### 2.3 Calculating Complex Effective Permittivity

To simulate the irregular surface of BSi, the cones are assumed to be asymmetric with slight variations in height and radius. A unit cell is then chosen and is assumed to be repeated laterally across the BSi surface, where the volume fill-factor of silicon in this unit cell is equivalent to that in the entire layer. Finding the fill-factor is, hence, a simple operation of dividing the volume of the silicon content in the unit cell by that of the cell, and is then used to find the effective complex permittivity of each layer using the EMT. Each of the stacked layers can be viewed as a composite of two materials: silicon and air. This composition can be approximated to a homogenous medium of certain effective permittivity, when the periodicity is much smaller than the wavelength of light. Many efforts have been exerted to find approximate formulas for the effective permittivity of composite mediums. Some of the most widely used formulas are:

- - General Lichtenecker mixing rule giving that [30,31]:
where

*ε*,_{e}*ε*and_{l}*ε*are the effective permittivity of the composite medium, the permittivity of the low-phase material (lower permittivity material) and the permittivity of the high-phase material (higher permittivity material), respectively, while V_{h}_{l}and V_{h}are their volume fractions and α is a parameter that determines the type of mixing (−1 ≤ α ≤ 1). - - Maxwell Garnett equation, which assumes spherical inclusions in a homogenous medium [32]:
where

*ε*and_{m}*ε*are the permittivity of the host medium and the permittivity of the inclusions, respectively, with_{i}*F*is the volume fill-factor of the inclusions. - - Bruggeman symmetric equation [29,30,33]:
where

*V*is the critical volume fraction. As discussed in the previous sub-section each layer of multilayer stack is assumed to contain truncated cones of silicon in air and are approximately treated as discs. In this case the silicon inclusions can be viewed as being distributed in a square lattice and hence Vc = 0.5 according to percolation theory [34]. The Bruggeman equation is the one chosen in our proposed model due to its superior agreement with other more rigorous methods for modeling BSi, which is discussed in detail in the next section._{c}

## 3. Proposed Model Validation

As mentioned in section 2, the proposed model relies on using the EMT in obtaining the effective index of the layers of the structure, and then the transfer matrix method is used to calculate the total reflectance, transmittance and absorptance of the multilayer structure. Since this is an approximate method the results of the proposed model are compared to that of RCWA to determine its range of validity.

#### 3.1 Implementing RCWA for Microstructured Si

The RCWA method was first proposed by M. G. Moharam and T. K. Gaylord for modeling light diffraction on planar gratings [35], and has been widely used and developed ever since in the literature [36-38]. The RCWA method is based on the exact solution of Maxwell’s equations for light diffraction on periodic gratings in both cases of planar and general conical diffraction. The solutions are obtained in the different simulated regions (the light input medium, grating region and light exit medium) and then the tangential electric and magnetic field components are matched at the boundaries [36]. For example, assume a 1D binary grating that is periodic along the x-axis (and having a normal along the z-axis), with *h* is the grating height, *Δ* is the grating period and *d* is the duty cycle. In the grating layer, where 0 < z < h, the permittivity *ε* can be expanded in a Fourier series in the following form [36]:

*ε*is the

_{m}*m*Fourier component of the permittivity in the grating layer. Generally, the RCWA is a stable method where the energy is always conserved and the accuracy of the obtained solution depends on the number of terms in the field-space expansion [36].The BSi and SiNWs structures considered in this work are simulated using RCWA under the same assumptions discussed in section 2 regarding the profile of the refractive index (shown in Fig. 3) and using the same formulas for the real and imaginary parts of the Si refractive index. The RCWA calculations are made using Reticolo software for grating analysis [40]. Once the total reflectance and transmittance are calculated, the absorptance and, hence the emissivity, is found as follows according to Kirchhoff’s law [39]:

^{th}#### 3.2 Comparing Proposed Model Results with RCWA Results

In this sub-section the simulation results obtained for different microstructured Si using the proposed model are compared with those obtained using the RCWA method. It should be noted that the proposed model is much simpler to implement and orders of magnitude computationally faster than the RCWA method for all the simulated structures.

### 3.2.1. Model Validity at Different Wavelengths

First, two BSi structures are simulated for the wavelength range of [0.5 – 10 μm]: the first structure has a maximum cone height of 1 μm, periodicity of 100 nm and base-spacing of 10 nm, while the second structure has a maximum cone height of 1.5 μm, periodicity of 150 nm and base-spacing of 15 nm. Both structures are assumed to be on the top of a Si substrate with a thickness of 500 μm and with donor doping concentration of 10^{19} cm^{−3}. The light is unpolarized and at normal incidence. The temperature is 300 K. The results of the proposed model and the RCWA method in terms of the total reflectance of the structure are compared in Fig. 6. The model yields results of high agreement with the RCWA, especially at longer wavelengths when the EMT is more accurate. The model is also able to capture the dependence of the reflectance on the geometry at small wavelengths, but the deviation between the two methods is larger due to the surface features being comparable to the wavelength.

### 3.2.2. Model Validity at Different Base-Spacing

Next, the reflectance of a BSi structure with height of 1 μm and periodicity of 0.15 μm is calculated but with the base-spacing varied from 0 to 0.1 μm. The results are shown in Fig. 7. It is noticed that as the base-spacing increases the results of the proposed model conform more to the RCWA method. This can be related to a fundamental assumption in the EMT formulas in general including the Bruggeman symmetric formula used in the proposed model. It treats each particle of the inclusions as an isolated individual in a homogenous medium, ignoring the interaction between the inclusion particles [41]. This assumption can lead to inaccurate results especially at high volume-fill factors of the inclusions in the composite medium, where these interactions are not negligible. This is further explained in Fig. 8, where the normalized electric field for another BSi structure is plotted using the RCWA method and the two cases of zero base-spacing and 0.08 μm spacing are compared. This structure is assumed to have a maximum cone height of 2 μm and periodicity of 0.2 μm. The donor doping concentration is 10^{15} cm^{−3}. The contour plots shown in Fig. 8 depict the magnitude of the electric field in a plane parallel to the surface normal of the simulated BSi structure. The plots demonstrate that in case of zero or small spacing (left figure) the resulting electric field distribution shows dependence on the surrounding BSi cones, while the in case of relatively larger spacing (right figure) the electric field distribution is confined around each cone and isn’t affected by the neighboring cones.

### 3.2.3. Model Validity at Oblique Incidence

Next, the proposed model is tested at different angles of incidence, where the incidence angle is swept from 0 to 88 degrees. The results are shown in Fig. 9 for a BSi structure of maximum cone height of 1 μm, periodicity of 0.15 μm, base-spacing of 0.05 μm and donor doping concentration of 10^{19} cm^{−3}. The proposed model shows good agreement at smaller incidence angle at the different wavelengths. The deviation of the proposed model from the RCWA method at larger angles of incidence can be related to the increased optical length in the simulated stacked layers, which is a function of the effective index calculated in each layer leading to larger error.

For the same BSi structure modeled in Fig. 9 the reflectance is calculated using the proposed model for the wavelength range of 1-10 μm and for angles of incidence of 0-90 degrees, the results are shown in Fig. 10 for S-polarized and P-polarized light. For S-polarized light it can be seen that reflectance increases as the angle of incidence increases for almost all wavelengths. While for P-polarized light the reflectance remains almost constant in the range of 0-60 degrees, then drops in the range 60-80 degrees (especially for longer wavelengths), then increases sharply in the range of 80-90 degrees.

## 4. Analysis of the Surface Topology Effect

In this section, the effect of some geometrical parameters of the BSi structure is studied using the proposed model. Also, the effects of the random surface topology on the optical response are investigated.

#### 4.1. Studying BSi Height and Base-spacing Variations

First, the effect of increasing the BSi height is studied by simulating a BSi structure with the periodicity of 0.25 μm and base-spacing of 0.05 μm with the height varied from 0.5 to 2 μm at wavelengths of 4, 5 and 6 μm. This simulation is run at 300 K, with a donor doping of 10^{18} cm^{−3}, at normal incidence for unpolarized light, and for a substrate thickness of 500 μm. Figure 11(a) shows that as the height increases the reflectance of the simulated samples decreases, which is attributed to the smoother variation of the effective index of the BSi structures that have larger aspect ratio between the height and period. This prediction agrees well with the behaviors experimentally reported in the literature [3,12].

Next the effect of increasing the BSi cones base-spacing is studied in Fig. 11(b) by simulating a BSi structure with the height of 1 μm, cone base-diameter of 0.4 μm with the base-spacing varied from 0 to 0.4 μm at wavelengths of 4, 5 and 6 μm. This simulation is run at 300 K, with a donor doping of 10^{18} cm^{−3}, at normal incidence for unpolarized light, and for a substrate thickness of 500 μm. Figure 11(b) shows that as the base-spacing between the BSi cones increases the reflectance of the simulated samples increases, which is attributed to the more abrupt change in the effective refractive index at the bottom of the BSi cones.

#### 4.2. Surface Features Disorder

The surface topology of BSi can vary according to the different fabrication processes as discussed earlier, and some structures can be more ordered than others [11]. For better modeling of these structures the random variations in the dimensions of the surface features of the fabricated samples should be considered. These variations can affect the effective refractive index of the assumed stacked layers and hence the reflectance, transmittance and emissivity of the structure.

An example of a BSi structure is shown in Fig. 12, where the structure is assumed to consist of cones of height 1.5 μm, periodicity of 0.3 μm, base-spacing of 0.05 μm, on a substrate of thickness of 500 μm, at normal incidence of unpolarized light, at a temperature of 300 K, and with donor doping concentration of 5x10^{19} cm^{−3}. Three cases are studied for this structure:

- - Case 1 for perfectly ordered structure with no disorder, where identical cones are assumed for the entire structure
- - Case 2 for low disorder, where the cone heights and radii can vary across the structure by 15%
- - Case 3 for higher disorder, where the cone heights and radii can vary across the structure by 30%

This disorder in the simulated structure requires the use of a larger unit cell of cones, where there’s a variation across this unit cell in the cones’ heights and radii within the percentages specified for the aforementioned cases. This results in different Si/air content in the stacked layer of each case, and hence can influence the simulation results.

It is noticed that for case 1 (no disorder) the reflectance is higher at shorter wavelengths but lower at longer wavelengths, while for case 3 (high disorder) the opposite occurs, and case 2 (low disorder) lies in between. The variations in the optical response of the shown 3 cases can be attributed to the different changes in the volume fill factor for each case as light moves through the BSi structure till it reaches the substrate. Also, there is an abrupt change in the refractive index where the bottom of the BSi cones meets the Si substrate, and this change is a function of the disorder assumed. Therefore, the proposed model can be used to predict the optical response more accurately for the different samples by altering the model input dimensions of the BSi features with variations resembling those of the sample.

## 5. Experimental Results

In this section, the proposed model is used to simulate different types of BSi (nanocones and SiNWs) and the results are compared with experimental results. First, a BSi sample of nanocones fabricated using cryogenic DRIE is measured with the help of an integrating sphere and an FTIR spectrometer, and the results are compared with those obtained using the model. Also, highly ordered SiNWs are simulated and the results are compared with those reported in the literature. It should be noted that the approach used in this model allows the simulation of other structures as well, such as pyramidal shapes, simply by adjusting the Si/air content in the assumed stacked layers.

#### 5.1 Black Silicon - Nanocones Results

The setup used to measure the reflectance of the BSi sample is shown in Fig. 13. The setup consists of a white light source emitting a wide spectrum in the near infrared range, lenses to collimate and focus the optical power on the BSi sample and an integrating sphere. The integrating sphere is coated on the inside with a highly diffuse reflectance material and is used to equally scatter the power reflected from the BSi sample over the sphere area. A portion of this power is coupled into an FTIR spectrometer. The spectrometer is used to measure the spectrum reflected from the BSi sample and this spectrum is divided by the spectrum obtained when the BSi sample is replaced by a closing port, which has the same coating as the sphere. Hence, the reflectance of the sample can be found for the entire wavelength spectrum. It should be noted that the BSi sample faces the light source with a small tilt angle of about 10 degrees so that the reflected power does not leave the sphere from the entrance port.

The BSi sample used in this measurement has the following parameters: a height of about 0.75 μm, period of 0.3 μm, almost no base-spacing and is fabricated on a lightly doped substrate of thickness 500 μm. An SEM of this BSi sample is shown in Fig. 14.

The spectrum obtained is shown in Fig. 15 and is compared to the simulation results using the proposed model. The TMM model shows good agreement with the measured data within 2%.

#### 5.2 Silicon Nano-wires Results

In this sub-section the proposed model results are compared to reflectance measurements reported in the literature for SiNWs of different dimensions [42-44]. The results are shown in Fig. 16. The SiNWs structures have the following parameters: Fig. 16(a) diameter = 38 nm, periodicity = 100 nm, height = 0.383 μm obtained from reference [42] and Fig. 16(b) diameter = 45 nm, periodicity = 100 nm, height = 0.383 μm obtained from reference [42] as well Fig. 16(c) diameter = 90 nm, periodicity = 400 nm, height = 1 μm obtained from reference [43]; Fig. 16(d) diameter = 105 nm, periodicity = 400 nm, height = 1 μm obtained from reference [44]. The angle of incidence is 65 degrees in all cases.

As expected, the agreement between the proposed model and the reported measurements depend on the structure dimensions, where in cases of smaller periods and diameters (Fig. 16(a) period = 100 μm, diameter = 38 nm, Fig. 16(b) period = 100 μm, diameter = 45 nm) there’s a good agreement for the whole simulated wavelength range, while in the cases of larger periods and diameters (Fig. 16(c) period = 400 μm, diameter = 90 nm, Fig. 16(d) period = 400 μm, diameter = 105 nm) there’s deviation at relatively shorter wavelengths and good agreement at relatively longer wavelengths. This can be related to the limitations of the EMT in modeling structures of dimensions (in this case the diameters of the SiNWs) comparable with the wavelength.

## 6. Conclusion

A method of modeling micro/nano-structured Si based on the TMM and the EMT was successfully presented. This model assumes stacked layers of a certain refractive index profile and each layer is assumed to be a homogenous medium of a certain complex effective index. The complex effective index is calculated using the EMT by means of Bruggeman symmetric formula. The reflectance of several simulated structures was compared to that obtained using the RCWA and compared with practical measurements. The model yields results of very good agreement specifically at longer wavelengths compared to the periodicity and for structures of relatively larger spacing between the surface features within the period. The proposed model was used to study the impact of the geometrical parameters of BSi including the structure disorder. The reported model is very efficient in terms of computational speed and can be used to predict and optimize the performance of BSi in the different applications including solar cells, photodetectors and Raman-scattering.

## References and links

**1. **J. Zhao and M. A. Green, “Optimized antireflection coatings for high-efficiency silicon solar cells,” IEEE Trans. Electron Dev. **38**(8), 1925–1934 (1991). [CrossRef]

**2. **H. Savin, P. Repo, G. von Gastrow, P. Ortega, E. Calle, M. Garín, and R. Alcubilla, “Black silicon solar cells with interdigitated back-contacts achieve 22.1% efficiency,” Nat. Nanotechnol. **10**(7), 624–628 (2015). [CrossRef] [PubMed]

**3. **K. N. Nguyen, P. Basset, F. Marty, Y. Leprince-Wang, and T. Bourouina, “On the optical and morphological properties of microstructured Black Silicon obtained by cryogenic-enhanced plasma reactive ion etching,” J. Appl. Phys. **113**(19), 194903 (2013). [CrossRef]

**4. **Y. Su, S. Li, Z. Wu, Y. Yang, Y. Jiang, J. Jiang, Z. Xiao, P. Zhang, and T. Zhang, “High responsivity MSM black silicon photodetector,” Mater. Sci. Semicond. Process. **16**(3), 619–624 (2013). [CrossRef]

**5. **M. Anwar, Y. Sabry, P. Basset, F. Marty, T. Bourouina, and D. Khalil, “Black silicon-based infrared radiation source,” Proc. SPIE **9752**, 97520E (2016). [CrossRef]

**6. **Y. M. Sabry, D. Khalil, T. E. Bourouina, and M. Anwar, “Structured silicon-based thermal emitter,” U.S. patent WO2017011269 A3 (2017).

**7. **W. Liu, A. Ming, Y. Ren, Q. Tan, W. Ou, X. Sun, W. Wang, D. Chen, and J. Xiong, “CMOS MEMS infrared source based on black silicon,” in Proceedings of IEEE 11th Annual International Conference on Nano/Micro Engineered and Molecular Systems (2016). [CrossRef]

**8. **Z. Xu, Y. Chen, M. R. Gartia, J. Jiang, and G. L. Liu, “Surface plasmon enhanced broadband spectrophotometry on black silver substrates,” Appl. Phys. Lett. **98**(24), 241904 (2011). [CrossRef]

**9. **E. P. Ivanova, J. Hasan, H. K. Webb, G. Gervinskas, S. Juodkazis, V. K. Truong, A. H. F. Wu, R. N. Lamb, V. A. Baulin, G. S. Watson, J. A. Watson, D. E. Mainwaring, and R. J. Crawford, “Bactericidal activity of black silicon,” Nat. Commun. **4**, 2838 (2013). [CrossRef] [PubMed]

**10. **T. Her, R. J. Finlay, C. Wu, S. Deliwala, and E. Mazur, “Microstructuring of silicon with femtosecond laser pulses,” Appl. Phys. Lett. **73**(12), 1673–1675 (1998). [CrossRef]

**11. **Y. J. Hung, S. L. Lee, K. C. Wu, Y. Tai, and Y. T. Pan, “Antireflective silicon surface with vertical-aligned silicon nanowires realized by simple wet chemical etching processes,” Opt. Express **19**(17), 15792–15802 (2011). [CrossRef] [PubMed]

**12. **X. Liu, P. R. Coxon, M. Peters, B. Hoex, J. M. Cole, and D. J. Fray, “Black silicon: fabrication methods properties and solar energy applications,” Energy Environ. Sci. **7**(10), 3223–3263 (2014). [CrossRef]

**13. **P. G. Maloney, P. Smith, V. King, C. Billman, M. Winkler, and E. Mazur, “Emissivity of microstructured silicon,” Appl. Opt. **49**(7), 1065–1068 (2010). [CrossRef] [PubMed]

**14. **M. Steglich, D. Lehr, S. Ratzsch, T. Kasebier, F. Schrempel, E.-B. Kley, and A. Tunnermann, “An ultra-black silicon absorber,” Laser Photonics Rev. **8**(2), L13–L17 (2014). [CrossRef]

**15. **S. Hava and M. Auslender, “Theoretical dependence of infrared absorption in bulk-doped silicon on carrier concentration,” Appl. Opt. **32**(7), 1122–1125 (1993). [CrossRef] [PubMed]

**16. **D. Abi Saab, S. Mostarshedi, P. Basset, S. Protat, D. Angelescu, and E. Richalot, “Effect of black silicon disordered structures distribution on its wideband reduced reflectance,” Mater. Res. Express **1**(4), 045045 (2014). [CrossRef]

**17. **N. M. Ravindra, S. R. Marthi, and S. Sekhri, “Modeling of optical properties of black silicon/crystalline silicon,” Silicon. J. Sci. Ind. Metrol. **1**(1), 100001 (2015).

**18. **T. Rahman and S. A. Boden, “Optical modeling of black silicon for solar cells using effective index techniques,” IEEE J. Photovolt. **7**(6), 1556–1562 (2017). [CrossRef]

**19. **D. Abi Saab, *Black Silicon Optical Properties, Growth Mechanisms and Applications (PhD thesis)*, Université Paris-Est, (2015).

**20. **S.-H. Hsu, E.-S. Liu, Y.-C. Chang, J. N. Hilfiker, Y. D. Kim, T. J. Kim, C.-J. Lin, and G.-R. Lin, “Characterization of Si nanorods by spectroscopic ellipsometry with efficient theoretical modeling,” Phys. Status Solidi., A Appl. Mater. Sci. **205**(4), 876–879 (2008). [CrossRef]

**21. **H. A. Macleod, *Thin Film Optical Filters* (CRC, 2010), Chap. 2.

**22. **H. H. Li, “Refractive index of silicon and germanium and its wavelength and temperature derivatives,” J. Phys. Chem. Ref. Data **9**(3), 561–658 (1980). [CrossRef]

**23. **B. L. Sopori, W. Chen, S. Abedrabbo, and N. M. Ravindra, “Modeling emissivity of rough and textured silicon wafers,” J. Electron. Mater. **27**(12), 1341–1346 (1998). [CrossRef]

**24. **A. Yariv and P. Yeh, *Photonics: Optical Electronics in Modern Communications* (Oxford University, 2006), Chap. 5.

**25. **G. E. Jellison Jr and D. H. Lowndes, “Optical absorption coefficient of silicon at 1.152 μ at elevated temperatures,” Appl. Phys. Lett. **41**(7), 594–596 (1982). [CrossRef]

**26. **G. G. Macfarlane, T. P. Mclean, J. E. Quarrington, and V. Roberts, “Fine structure in the absorption-edge spectrum of Si,” Phys. Rev. **111**(5), 1245–1254 (1958). [CrossRef]

**27. **J. C. Strum and C. M. Reaves, “Silicon temperature measurement by infrared absorption. Fundamental processes and doping effects,” IEEE Trans. Electron Dev. **39**(1), 81–88 (1992). [CrossRef]

**28. **P. Vandenabeele and K. Maex, “Influence of temperature and backside roughness on the emissivity of Si wafers during rapid thermal processing,” J. Appl. Phys. **72**(12), 5867–5875 (1992). [CrossRef]

**29. **A. R. Abramsona, P. Nieva, H. Tada, P. Zavracky, I. N. Miaoulis, and P. Y. Wong, “Effect of doping level during rapid thermal processing of multilayer structures,” J. Mater. Res. **14**(6), 2402–2410 (1999). [CrossRef]

**30. **Y. Wu, X. Zhao, F. Li, and Z. Fan, “Evaluation of mixing rules for dielectric constants of composite dielectrics by MC-FEM calculation on 3D cubic lattice,” J. Electroceram. **11**(3), 227–239 (2003). [CrossRef]

**31. **K. Lichtenecker, “The dielectric constants of natural and synthetic mixtures,” Phys. Z. **27**, 115 (1926).

**32. **J. C. M. Garnett, “Colours in metal glasses and in metallic films,” Phil. Trans. R. Soc. A **203**(359-371), 385–420 (1904). [CrossRef]

**33. **D. A. G. Bruggeman, “Berechnung verschiedener physikalischer konstanten von heterogenen substanzen. I. Dielektrizitätskonstanten und leitfähigkeiten der mischkörper aus isotropen substanzen,” Ann. Phys. **416**(7), 636–664 (1935). [CrossRef]

**34. **R. Zallen, *The Physics of Amorphous Solids* (John Wiley & Sons, 1983), Chap. 4.

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

**36. **M. G. Moharam, E. B. Grann, and D. A. Pommet, “Formulation for stable and efficient implementation of the rigorous coupled-wave analysis of binary gratings,” Opt. Soc. Am. **12**(5), 1068–1076 (1995). [CrossRef]

**37. **P. Lalanne and G. M. Morris, “Highly improved convergence of the coupled-wave method for TM polarization,” J. Opt. Soc. Am. **13**(4), 779–784 (1996). [CrossRef]

**38. **P. Lalanne and M. P. Jurek, “Computation of the near-field pattern with the coupled-wave method for transverse magnetic polarization,” J. Mod. Opt. **45**(7), 1357–1374 (1998). [CrossRef]

**39. **Y. A. Cengel, *Heat and Mass Transfer: A Practical Approach* (McGraw Hill, 2006), Chap. 11.

**40. **J. P. Hugonin and P. Lalanne, (2005), *Reticolo Software for Grating Analysis*, Institut d'Optique, Orsay, France.

**41. **V. Myroshnychenko and C. Brosseau, “Finite-element modeling method for the prediction of the complex effective permittivity of two-phase random statistically isotropic heterostructures,” J. Appl. Phys. **97**(4), 044101 (2005). [CrossRef]

**42. **M. Khorasaninejad, N. Abedzadeh, J. Sun, J. N. Hilfiker, and S. S. Saini, “Polarization resolved reflection from ordered vertical silicon nanowire arrays,” Opt. Lett. **37**(14), 2961–2963 (2012). [CrossRef] [PubMed]

**43. **S. Patchett, M. Khorasaninejad, N. O, and S. S. Saini, “Effective index approximation for ordered silicon nanowire arrays,” J. Opt. Soc. Am. B **30**(2), 306–313 (2013). [CrossRef]

**44. **M. Khorasaninejad, S. Patchett, J. Sun, N. O, and S. S. Saini, “Polarization-resolved reflections in ordered and bunched silicon nanowire arrays,” J. Opt. Soc. Am. B **9**(11), 3063–3068 (2012). [CrossRef]