Anisotropic Complex Refractive Indices of Atomically Thin Materials: Determination of the Optical Constants of Few-Layer Black Phosphorus

Abstract

In this work, studies of the optical constants of monolayer transition metal dichalcogenides and few-layer black phosphorus are briefly reviewed, with particular emphasis on the complex dielectric function and refractive index. Specifically, an estimate of the complex index of refraction of phosphorene and few-layer black phosphorus is given. The complex index of refraction of this material was extracted from differential reflectance data reported in the literature by employing a constrained Kramers–Kronig analysis combined with the transfer matrix method. The reflectance contrast of 1–3 layers of black phosphorus on a silicon dioxide/silicon substrate was then calculated using the extracted complex indices of refraction.

1. Introduction

Two-dimensional materials emerged as a promising option for the development of nanoelectronic and optoelectronic devices, due to their peculiar electronic and optical properties. In this respect, transition metal dichalcogenides (TMDs), when thinned to a single atomic layer, show an indirect–direct band gap transition, inequivalent valleys in the Brillouin zone and non-trivial topological order [1,2,3]. In the vast family of two-dimensional materials, elemental ones represent a niche, due to their high reactivity at ambient conditions, that makes their characterization and exploitation challenging [4,5]. Among them, phosphorene, i.e., monolayer black phosphorous (BP), is the most promising material for optoelectronic and photonic applications, as its electronic band gap lies in a spectral region between the ones of TMDs and the one of graphene [6]. For an effective employment of these materials in photonic and optoelectronic devices, it is important to determine their dielectric response.

Among the monolayers of transition metal dichalcogenides, MoS₂ was the first material on which optical measurements were performed. In fact, in 2013 and 2014, the complex dielectric function of MoS₂ was measured with spectroscopic ellipsometry by Shen et al. [7], Yim et al. [8] and Li et al. [9]. The complex dielectric functions of monolayers MoS₂, WS₂, MoSe₂ and WSe₂ were measured by Li et al. by employing reflectance contrast spectroscopy [10]. In 2019, the monolayers, together with the bilayers and the trilayers, of the four above mentioned chalcogenides were studied by Hsu et al. [11]. In 2020, Ermolaev et al. [12] measured the complex dielectric function of MoS₂ in a broad range of wavelengths, i.e., 290–3300 nm. Ermolaev et al. [13] also measured the complex index of refraction of WS₂ in the range of 375–1700 nm.

Additionally, a number of studies investigated the anisotropic components of TMD and phosphorous few-layer systems via the incorporation of dielectric and plasmonic nanostructures [14,15,16,17,18,19,20,21]. In one study, the coupling of electromagnetic fields with the out-of-plane component of the MoS₂ dielectric function was possible due to the peculiar morphology of chemical vapor deposition (CVD)-grown few-layer MoS₂ on rippled substrate, hence promoting 1D nanostructures as a promising path for the full exploitation of the few-layer MoS₂ dielectric response [19]. Similarly, a top-down approach was revealed to be effective for the strain-induced modification of the electronic band structure of atomically thin MoS₂ and black phosphorous, paving the way for the on-demand tuning of their optical properties [20,21]. One study investigated exfoliated MoS₂ on a lithographically defined SiO2 nanocone substrate, demonstrating a deterministic red-shift of the A exciton photoluminescence (PL) due to elastic strain [14]. Another study utilized spatially and time-resolved PL diffusion measurements of WSe₂ exfoliated onto a 1.5 µm diameter SiO₂ pillar, resulting in a similar red-shift and demonstration of exciton diffusion towards high tensile strain regions; that work was a proof of concept for engineerable excitonic diffusion that may lead to a new generation of opto-excitonic devices [15].

Strain engineering may also be used to enable new radiative pathways for the optical excitation of selection rule forbidden TMD dark excitons under normal incidence conditions; these dark and “gray” excitons possess much longer radiative lifetimes (>100 ps) [16,22,23]. In the absence of intentional strain, one study utilized a high numerical aperture objective lens (NA = 0.82) to characterize the dark exciton in WSe₂: even for normal incidence, a Gaussian beam with a tight focus has a non-zero electric field component along the beam propagation direction [16]. However, optical excitation of the dark exciton in WSe₂ was more easily enabled by coupling to the surface plasmon polariton of silver (Ag) [17]. Another study incorporated a hexagonal boron nitride (hBN)-encapsulated WSe2 monolayer into a plasmonic modulator device, requiring precise knowledge of the TMD, gold (Au) and hBN complex indices of refraction [18]. It thus emerges how the development of simple and robust tools for the careful analysis of the dispersion of the complex refractive index is beneficial for an appropriate engineering of devices that include two-dimensional materials.

In contrast, phosphorene (single-layer BP) and few-layer BP have not been studied with the same fervor as their TMD counterparts, likely due to their relative sensitivity to degradation. One significant distinction between TMDs and BP is that the band structure of BP exhibits considerable anisotropy: that is, the effective electron mass along the armchair and zig-zag directions is a factor of five times larger in the former direction [24]. This anisotropy is manifest in polarization-sensitive absorption and PL measurements [24,25,26,27], as well as spatially resolved excitonic diffusion measurements [24]. The relatively fragile phosphorene system holds considerable promise due to the ease of PL/absorption tunability by layer number tuning, which can shift the resonance energies between 0.32 and 1.7 eV [24].

We focus on estimating the complex refractive index of phosphorene, which shows a morphological in-plane anisotropy resulting in an anisotropic optical response [25]. We retrieve the refractive index along the armchair (AC) and zig-zag (ZZ) direction by exploiting the constrained Kramers–Kronig analysis (CKKA) and fit experimental contrast reflectivity with a model obtained by means of the transfer matrix method (TMM).

2. Materials and Methods

Methodology for extraction of the complex index of refraction of few-layer BP from reflectance contrast data: The linear optical properties of thin films are commonly measured via static absorption or reflectivity contrast measurements. Experimental techniques include normal incidence reflectivity contrast (dR/R) over a wide plane of view, spatially resolved confocal dR/R in the case of exfoliated small-area TMDs, few-layer BP, graphene and graphene oxide and simultaneous absolute reflectivity and transmission measurements for the determination of static absorption [10,11,25,27,28,29,30]. Although these methods reveal the static absorption A and reflectivity dR/R and transmission dT/T contrasts of the thin films relative to the linear substrate response, the determination of the complex index of refraction $\tilde{n} = n+ik$, or equivalently the complex dielectric function $\epsilon = {\epsilon }_r+i{\epsilon }_i$, where $n = \sqrt{\epsilon }$, typically requires more sophisticated optical characterization methods such as ellipsometry [8,31,32,33].

In this section, we extract the complex index of refraction of hBN-encapsulated few-layer BP on a sapphire substrate from confocal reflectivity contrast data taken from [25]. Specifically, the layered sample structure is as follows (see Figure 1): air → 15 nm hBN → few-layer BP → $\alpha$-Al₂O₃ (sapphire); BP thicknesses range nominally from 0.5 nm (1L) to 2.5 nm (5L), and bulk BP is 100 nm thick. We assume that the exfoliated hBN layer is crystalline, rather than isotropic as for the case of cubic BN or polycrystalline hBN [33,34,35] with an in-plane dielectric constant given by

$${\epsilon }_{hBN,\parallel } = n_{hBN,\parallel }^2 = 1+ \frac{3.336 {\lambda }^2}{{\lambda }^2-26322}$$

(1)

where the wavelength $\lambda$ is given in nanometers [33]. Additionally, we also use a Sellmeier equation to describe the $\alpha$-Al₂O₃ substrate index of refraction [36].

One commonly used method for the determination of the complex linear material response functions is the Kramers–Kronig analysis (KKA): it is well known that the real and imaginary parts of the dielectric function are related to one another via the Kramers–Kronig (KK) relations, which impose causality conditions on the response functions [37,38,39]. However, the integrals involved in these relations are taken over infinite limits; practical implementation requires the acquisition of data over a wide energy range, and may be limited to materials which exhibit an optical response only over a narrow wavelength range.

We utilize an alternative method called the constrained Kramers–Kronig analysis (CKKA), which was introduced by Kuzmenko to address some of the problems associated with KKA [37]. This method takes advantage of the fact that basis functions such as the complex Lorentzian can be chosen that individually satisfy the KK relations. Such a function is given by

$${\epsilon }_k = \frac{{\omega }_{p,k}^2}{{\omega }_{0,k}^2-{\omega }^2+i\omega {\gamma }_k}$$

(2)

where ${\omega }_{p,k}^2$ is the plasma frequency, ${\omega }_{0,k}^2$ is the oscillator frequency and ${\gamma }_k$ is the linewidth. Since linear combinations of such basis functions also satisfy the KK relations, the method introduces a mesh in energy space of a large number of oscillators $N_{osc}$~$N_{data}$; the oscillator frequencies are fixed, and the linewidth is chosen as the spacing between oscillator energies. Thus, the only remaining parameters left to minimize are the oscillator strengths (or plasma frequencies): hence, the constrained KK analysis. To further optimize the CKKA fitting procedure, we chose functions that are similar to the complex Lorentzian of Equation (2) and also satisfy the KK relations, but are more locally weighted, avoiding the contributions of the tails of the Lorentzian distributions far away from the central oscillator energy. One such function is a triangular function for ${\epsilon }_i$, with a corresponding analytically integrable ${\epsilon }_r$ [37,38].

We note that although this method can easily reproduce reflectivity contrast data when a sufficiently large oscillator mesh is used, it is sensitive to a lack of information about the optical response outside of the experimental acquisition range [10,37,38]. This problem is not entirely unexpected, considering that it is fundamentally a KK-type analysis method. Often, high-frequency information about a given optically active system is known via X-ray photoelectron emission spectroscopy or UV absorption experiments; additionally, low-frequency information is often revealed via Fourier transform IR (FTIR) and Raman scattering experiments [5,40,41,42]. However, such information is currently lacking for few-layer BP, although bulk BP is well characterized. Although similar calculations using the CKKA method have extracted the complex index of refraction for TMD monolayer and few-layer systems assuming a bulk-like response at high energies [10], no such assumptions are made here for the few-layer BP fitting procedure.

For a given complex index of refraction, the optical response of a multi-layered system is then determined by the transfer matrix method (TMM) [43]. It is noted that perturbative “linearized” methods have also been developed for the treatment of TMD monolayers, few-layer BP and graphene [10,25,30], which treat the optical response of the thin film using the sheet conductivity

$${\sigma }^d = \sigma d = -id\omega \left(\epsilon -1\right).$$

(3)

For basic layered systems such as air → monolayer TMD → SiO₂, relatively simple forms for the reflectivity contrast dR/R can be derived [29]. However, for the system under study in this section (air → hBN → BP → $\alpha$-Al₂O₃), with no assumptions being made about reflection coefficients between interfaces, no tractable analytical form for dR/R is derived; note that this stands in contrast to previous studies [25] which assume that the BP → $\alpha$-Al₂O₃ interface does not contribute to the overall optical response of the layered system.

We use the standard TMM [43] to relate the input and reflected fields at the air → hBN interface (1) to the transmitted fields at the BP → $\alpha$-Al₂O₃ interface (3). Specifically, for normal incidence fields, the electric and magnetic fields, which are implicit sums of all fields (incident, reflected, transmitted) are given by

$$\left[\begin{array}{c}{E}_1\\ H_1\end{array}\right] = {\widehat{M}}_{1-2}{\widehat{M}}_{2-3}\left[\begin{array}{c}{E}_3\\ H_3\end{array}\right], {\widehat{M}}_{i-j} = \left[\begin{array}{cc}\cos {\varphi }_{ij}& \frac{i}{n_{ij}}\sin {\varphi }_{ij}\\ i{n}_{ij}\sin {\varphi }_{ij}& \cos {\varphi }_{ij}\end{array}\right], {\varphi }_{ij} = \frac{2\pi n_0{n}_{ij}{d}_{ij}E}{hc}$$

(4)

which can then be related to the field reflection coefficient by the following relation

$$r = \frac{E_r}{E_i} = \frac{n_0{m}_{11}+n_0{n}_s{m}_{12}-m_{21}-n_s{m}_{22}}{n_0{m}_{11}+n_0{n}_s{m}_{12}+m_{21}+n_s{m}_{22}}$$

(5)

where $m_{ij}$ is the ij-th element of the product of the two matrices in Equation (4), $n_0,n_s$ are the indices of refraction of air and the back substrate ($\alpha$-Al₂O₃), $hc$ is the product of the speed of light and Planck’s constant, $E = \raisebox{1ex}{$hc$}\!\left/ \!\raisebox{-1ex}{$\lambda $}\right.$ is energy, $n_{ij}$ is the index of refraction for the material layer between the i and j interfaces and $d_{ij}$ is the thickness between the i and j interfaces (Figure 1). Finally, the reflectance contrast is given by

$$\raisebox{1ex}{$dR$}\!\left/ \!\raisebox{-1ex}{$R$}\right. = \left({\left|r_{ref}\right|}^2-{\left|r_{sample}\right|}^2\right)/{\left|r_{ref}\right|}^2$$

(6)

where $r_{ref}$ is the field reflection coefficient for the reference air → hBN → $\alpha$-Al₂O₃, which is determined by setting ${\widehat{M}}_{2-3}= \widehat{1}$. This procedure makes no perturbative assumptions about the index of refraction of BP, and allows for generalization to thicker layers such as the bulk BP case.

3. Results and Discussion

Complex index of refraction for few-layer and bulk BP: The complex index of refraction for few-layer BP is then determined by fitting the reflectance contrast data [25] using non-linear least squares fitting of the weight (plasma frequencies) parameters in the CKKA method with triangular basis functions; this method is referred to as CKKA+TMM [37,38]. The dR/R data were extracted from the literature using a plot digitizer with high accuracy [44]. The resulting dR/R fits, real indices of refraction n and extinction coefficients k, for both x- and y-polarized excitation for one to five layers and bulk BP, are displayed in Figure 2. We note the high quality of the dR/R fits (Figure 2a–e), even small variations in dR/R are easily reproduced with CKKA+TMM; a high-resolution mesh of 300 triangular oscillators was chosen to fit these data. We argue that the constrained KK method utilized here provides credibility to the resulting complex indices of refraction, since the dispersive component (n, or ${\epsilon }_r$) is connected directly via the KK relations to the absorptive component (k, or ${\epsilon }_i$), which for atomically thin films is the major contribution to the dR/R measurement [10,29].

As it has been reported previously [25], the reflectance contrast data for 1–5L BP (taken at 77 K) shows either one or two peaks within the experimental acquisition range of 0.75–2.4 eV, corresponding to both the intralayer exciton at lower energies and a blue-shifted exciton that arises due to interlayer effects similar to quantum wells [24]. We note that x- and y-polarized excitation here corresponds to polarization aligned along the AC and ZZ directions of the puckered few-layer BP, respectively. This puckering effect leads to significant in-plane anisotropy, resulting in a band structure for which the effective electron mass along the AC direction is five times larger than the ZZ direction [24].

As revealed both in the dR/R data, as well as in our extracted complex indices of refraction, sharp excitonic features are observed only in the AC (x) case in dR/R and k, with strongly dispersive lineshapes in n imposed on a relatively flat background which varies from 2–4 (Figure 2k–o). In the case of ZZ excitation, both the index n and extinction coefficient k are relatively featureless: the index is relatively flat between 2.5–3 from 1–2.2 eV in the cases of 2–5 L (Figure 2g–j, pink dashed curve), while the extinction coefficient k increases from 0 to 4 between these energy ranges (Figure 2l–o, pink dashed curve). In fact, for both n and k, the results for ZZ excitation seem to constitute the background features in the AC case, indicating the presence of shared higher energy resonances which are isotropic. Thus, by examining the extinction coefficient k directly, which corresponds to absorption, we confirm previous statements in the literature that there exists a broad and increasing background absorption in few-layer BP towards higher energy, on which the strong and narrow excitonic resonances are super-imposed.

The CKKA+TMM analysis used here is more general than the linearized models often used for the treatment of static absorption by atomically thin materials such as TMDs and few-layer BP [10,25]. Indeed, this model is used here to extract quantitatively meaningful results for the complex index of refraction, without resorting to approximations such as small optical path length for the encapsulation layer, minimal absorption in the BP layers or the neglect of the substrate reflection at the BP → $\alpha$-Al₂O₃ interface.

This generality allows us to also fit the 100 nm bulk BP sample (Figure 2, fourth row). The dR/R data are no longer easily directly correlated with the absorption coefficient k, due to the non-negligible thickness of bulk BP, as well as the interference effects with the hBN encapsulating layer. High quality fits of the dR/R data reveal again large anisotropy resulting in birefringence of bulk BP, with a modulation of the index of refraction around 0.32 eV in the AC case, which is not observed in the ZZ excitation case. For the extinction coefficient k, in the AC case, the absorption shows an onset around the band edge at 0.32 eV, but also shows absorption at lower energies likely due to the Drude response from free carrier excitation [24]. Surprisingly, for the ZZ excitation case, the absorption also increases steadily towards higher energies, without any dispersion modulation around the band edge.

We also note that in this energy range, the extracted index of refraction for bulk BP is below 1. Sub-unity indices of refraction are commonly reported in noble metal materials such as Ag and Au, where there is an interplay between interband absorption and the plasmonic response [45,46,47] as well as in Al thin films at wavelengths shorter than 660 nm [48]. Indeed, a zero-crossing of the real part of the dielectric function ${\epsilon }_r$ is required to satisfy the conditions for a localized surface plasmon resonance (LSPR) [49]; the plasmonic nature of BP is currently under investigation, due its novel features such as hyperbolic plasmonics arising from the large ratio between the AC and ZZ direction effective electron mass [24].

Systematics: In this section, we remark on the sensitivity of the CKKA+TMM analysis for the extraction of the complex indices of refraction to variations in parameters such as substrate and encapsulation layer refractive indices, as well as fixed fitting parameters. Previous reports involving the utilization of the CKKA method have commented on the need for incorporating information about low and high energy extremes of the dielectric functions outside of the experimental acquisition range. For instance, Li et al. [10] assumed during their studies of TMD monolayers that at high energies (>3 eV), the bulk optical properties also contributed to the monolayer optical response; absorption resonances were included up to 30 eV. No such assumptions are made in our determination of n,k for few-layer BP or the bulk BP.

However, in Figure 3, a convergence of the fitting results for 1L-x BP is shown as a function of the buffer energy that is included in the fits. More specifically, the experimental data were collected from 0.75–2.4 eV; additionally, we examined the convergence of the fit solution as a function of how far the oscillator energy mesh extended beyond these extremal values, i.e., $E_{mesh}\in \left[E_{data,i} - E_{buf},E_{data,f}+E_{buf}\right]$. Although the imaginary part of the CKKA basis triangular dielectric function is well localized in energy, the real part exhibits long tails away from the oscillator center; thus, by extending our oscillator energy range, we can simulate how these low and high energy contributions increase the quality of fit. This buffer energy was tested from 0 to 0.3 eV; it is shown (Figure 3, center, inset) that increasing the buffer energy from 0 to 0.1 eV reduces the sum of the residuals by nearly two orders of magnitude, with generally no enhancement beyond 0.1 eV. This convergence can be seen in the results for n and k. However, it is noted that at energies greater than 2.4 eV that the real index of refraction n varies considerably for different buffer energies; random variations are also observed in this energy range for slight variations in the number of energy mesh points (not shown). We believe that this lack of convergence arises because of the lack of constraints on the dielectric function we have imposed on energies higher than 2.5 eV; this aspect may be potentially improved by using the higher energy absorption data for bulk BP, in a similar manner as done with TMD monolayers.

The dependence of the results for 1L-x BP were also studied as a function of the hBN refractive index. It is known that exfoliated hBN is birefringent, but estimates of the for in-plane and out-of-plane indices of refraction have varied significantly between reports [33,34,35]. For the results displayed in Figure 2, as well as Figure 3, a Sellmeier equation is used that describes the index of refraction of the in-plane index, which ranges from 2.25 to 2.1 from 400–1600 nm [33]. The fitting procedure was also performed for a constant index n = 2.108, or the average of the in-plane index, and n = 1.86, a commonly used value for isotropic BN [33,34,35] (Figure 3a,b). Using the constant in-plane value yields a generally lower value than the Sellmeier case, and using the isotropic value leads to much larger estimates of n; the variations in k are not as significant (Figure 3c). We believe that the use of the functional Sellmeier equation form is appropriate; the measurement resulting in that form utilized confocal oblique incidence ellipsometry to capture both the in-plane and out-of-plane indices from 400–1600 nm [33].

Calculation of reflectance contrast of BP on SiO₂/Si substrates: One immediately apparent use for the complex indices of refraction of BP determined in this paper using CKKA+TMM is the simulation of the optical properties of complex multi-layered devices incorporating BP. During the characterization and development of novel optical devices that incorporate exfoliated materials such as TMDs and BP, one step in the standard procedure for exfoliation involves the isolation of small flakes by thickness. This identification often involves measuring the reflectance contrast over the visible spectrum, i.e., examining a microscope image by color contrast (by eye or CCD) [50]. Thus, maximizing the color contrast between the active layer (TMD or BP) and the substrate is essential for this exfoliation procedure. If the reflectance contrast curves can be predicted ahead of time given knowledge of the complex indices of refraction of the TMD/BP and the substrate, the thickness of the active layer can be determined with certainty.

In TMD systems, one common substrate employed for high color contrast is 90–120 nm of SiO₂ on a Si substrate [28]. Here, we employed the TMM to calculate the expected reflectance contrast of one to three layers of BP on that same substrate, using the complex dielectric constants determined in this article via CKKA+TMM. Thus, we have provided a useful reference for optical physicists attempting to incorporate BP in their devices. The results of the calculation are shown in Figure 4. It becomes immediately clear that this standard substrate does not provide high reflectance contrast, at least not as high as using sapphire as a substrate (Figure 2). This lack of reflectance contrast is due to the relative similarities of the real parts of the index of refraction for Si and 1–3 L BP, with n ~ 3–4 from 0.8 to 2.5 eV (Figure 2). However, the contrast between the x- and y-polarized excitation reflectance contrast allows one to clearly locate sharp features in the optical spectrum, which can then be used for layer thickness determination. For easy thickness determination, we suggest utilizing either sapphire as a substrate, or polarization-dependent reflectance contrast.

Discussion of study limitations, and comparisons to the literature: Some of the limitations of this study have been discussed above, and we re-iterate in more detail here. First, as discussed in detail by Li et al. [10] and Kuzmenko [37], both the standard Kramers–Kronig method and the CKKA method utilized in this article are sensitive to the optical response of the system outside of the experimental acquisition range. The standard KK approach requires that the absorptive optical response is localized over an energy range smaller than the experimental energy range; while the CKKA method does not explicitly require this, we showed in the Systematics section of this article that the extracted complex index of refraction has high uncertainty at the edges of the data acquisition range. These uncertainties would be greatly reduced by incorporating both higher and lower energy data that could be acquired via electron energy-loss measurements and UV-XUV absorption measurements [41], as well as long wavelength FTIR spectroscopy. We note that it is non-trivial to incorporate the results of those higher/lower energy experiments into CKKA+TMM, since those experiments span over highly disparate energy ranges and apply to the bulk crystal black phosphorous, but point the reader towards referenced literature for future studies [51,52,53,54].

Second, both the KK and CKKA methods require precise characterization of the thicknesses of the substrate layers with nanometer precision, and knowledge of the indices of refraction for those layers. In this article, the hBN layer thickness has been taken as approximately 15 nm, as given by Li et al. [25]. We suggest that a more consistent approach would involve multiple reflectance contrast measurements of BP with varying hBN encapsulation layer thicknesses, with thicknesses of all layers characterized by atomic force microscopy. The same concern applies to the extraction of the complex index of the bulk BP described in this article; this thickness was taken to be equal to approximately 100 nm, as given by Li et al. [25], but the uncertainty on this thickness is unknown to us. Additionally, we are uncertain of the crystallinity of the hBN encapsulation layer, as discussed in the Systematics section of this article, but have provided systematic variations over the possible cases.

In the context of the literature regarding the optical characterization of TMDs and few-layer BP, and more generally exfoliated materials, this article takes a similar approach to that of Li et al. [10], in which the complex dielectric functions of monolayer TMDs were determined via CKKA+TMM. The analysis provided here clearly separates the dispersive (real) and absorptive (imaginary) parts of the optical dielectric constants of few-layer and bulk BP, which is not obvious in reflectance contrast data.

4. Conclusions

For the estimate of the index of the refraction of black phosphorus, we have demonstrated that the constrained KK analysis combined with a general transfer matrix method (CKKA+TMM) that takes into account finite layer optical path lengths can be used to extract the entire complex index of the refraction of 1–5 L BP and bulk BP from reflectance contrast data. It is confirmed that the excitonic peaks correspond to peaks in the extinction coefficient k, and that a significant absorptive background that increases towards higher energies exists. In contrast, the real index of refraction is relatively flat with superimposed dispersive features located at the excitonic resonances. The determination of the full complex index of refraction $\tilde{n}$ of few-layer BP is significant: these results may be incorporated into the development of novel optical devices such as responsive photonic crystals and distributed Bragg reflectors (DBRs) [55,56,57,58], plasmonic nanostructures [59] and on-chip waveguides, all of which require knowledge of both the dispersive and absorptive optical properties. We calculated the reflectance contrast of 1–3 L BP on a SiO₂/Si substrate, providing a useful reference for physicists working to determine thicknesses of exfoliated BP. We suggest further improvements to this extractive method that would involve the experimental acquisition of both reflectance and transmission contrast simultaneously, as well as constraining the low and high energy ranges in the CKKA+TMM analysis with optical properties deduced from FTIR and UV absorption.