Brought to you by:
Paper

A simple and robust machine learning assisted process flow for the layer number identification of TMDs using optical contrast spectroscopy

, and

Published 14 November 2022 © 2022 IOP Publishing Ltd
, , Citation Nikhil Joseph Joy et al 2023 J. Phys.: Condens. Matter 35 025901 DOI 10.1088/1361-648X/ac9f96

0953-8984/35/2/025901

Abstract

Layered transition metal dichalcogenides (TMDs) like tungsten disulphide (WS2) possess a large direct electronic band gap (∼2 eV) in the monolayer limit, making them ideal candidates for opto-electronic applications. The size and nature of the bandgap is strongly dependent on the number of layers. However, different TMDs require different experimental tools under specific conditions to accurately determine the number of layers. Here, we identify the number of layers of WS2 exfoliated on top of SiO2/Si wafer from optical images using the variation of optical contrast with thickness. Optical contrast is a universal feature that can be easily extracted from digital images. But fine variations in the optical images due to different capturing conditions often lead to inaccurate layer number determination. In this paper, we have implemented a simple Machine Learning assisted image processing workflow that uses image segmentation to eliminate this difficulty. The workflow developed for WS2 is also demonstrated on MoS2, graphene and h–BN, showing its applicability across various classes of 2D materials. A graphical user interface is provided to enhance the adoption of this technique in the 2D materials research community.

Export citation and abstract BibTeX RIS

1. Introduction

The isolation of the two-dimensional (2D) material, graphene by mechanical exfoliation and its subsequent characterization in 2004 has promoted a paradigm shift in materials research [1]. The unique and fascinating properties exhibited by graphene acted as a catalyst for exploring new 2D materials (2DMs), and at present, the class of 2DMs are so vast that it covers graphene, elemental compounds known as Xenes, h-BN, black phosphorus, the family of transition metal dichalcogenides, layered double hydroxides, perovskite type oxides and the newly synthesized class of materials without a known 3D form like MoSi2N4, WSi2N4 [2, 3]. The new classes of 2DMs are shown to have excellent physical, chemical and mechanical properties, which make them ideal candidates for many applications [49]. Among them, transition metal dichalcogenides (TMDs) are interesting due to a finite bandgap in the monolayer limit [10, 11]. In addition, most TMDs show an indirect to direct bandgap transition as the monolayer limit is approached [12, 13]. For instance, TMDs like WS2 and MoS2 have a direct bandgap of the order 1 eV in the monolayer [14] and are shown to exhibit superconductivity in the monolayer limit in WS2 [15]. These make them ideal candidates for (a) opto-electronic applications, and (b) future nanoelectronic devices with atomic-scale thickness. Since most properties of these TMDs are strongly dependent on the number of layers, layer number determination becomes a crucial step in the realization of these two-dimensional materials in various device applications.

The most common method for thickness identification of thin layer samples is atomic force microscopy (AFM) [16]. However, the method is time consuming and is not suitable for mechanically exfoliated flakes. AFM tends to wrongly estimate the number of layers since it does not account for the substrate–sample interaction [17]. Raman spectroscopy is a more reliable tool and is widely used for the layer number determination of various 2DMs like graphene, h-BN and TMDs. But for each material, we need to consider different parameters like intensity, intensity ratio, peak position, the difference in the position of nearby peaks, full width at half maximum (FWHM) and so on. The excitation wavelength also plays a crucial role in this method since the layer-dependent features are not clear in certain wavelength excitations. For example, the Raman spectra of WS2 under 488 nm excitation give clear layer-dependent signatures, whereas for 532 nm excitation it does not. Thus, when it comes to TMDs, photoluminescence (PL) spectra contain layer-dependent signatures, at-least for layer number $\lt$5 [18, 19]. But as the layer number goes higher, the PL intensity of most of the TMDs like MoS2 and WS2 decreases drastically and thus limits its use.

The need for different experimental tools under particular conditions to characterize various TMDs becomes a bottleneck for accelerated research in this field. Utilizing the features in the optical images can be a universal alternative for determining the number of layers. The conventional approach in optical identification relies on the optical contrast of the samples and the field is called optical contrast spectroscopy [20, 21]. The observed contrast can be correlated with the layer number using standard characterization techniques such as Raman, PL or AFM. This relation can then be utilized for all further layer number determination. However, random averaging of the pixel values and differences in laboratory conditions while capturing the image makes this approach volatile. The slight variations in the intensity of light coming from uniform regions of the sample can lead to false identification of the number of layers. There have been several attempts trying to take into account various parameters like the SiO2 thickness, illumination effects etc that makes the problem complex. Ouyang et al concludes that there can be no universal function that relates contrast with the layer number for all laboratory conditions suggesting the need for custom functions [22]. The modern approach involves the classification of optical images using Machine Learning (ML) techniques for automatic detection of the number of layers [2330]. This is often complex to implement and computationally intense. Reliable performance from these methods requires a large labeled dataset for training. Construction of a labeled training set is a formidable task since the yield of mechanical exfoliation is much lower for TMDs in comparison with graphene. This makes the method unsuitable for beginners without a large training set. This is also prone to laboratory-dependent variations. In this paper, we propose a new approach to overcome the above difficulties in optical contrast identification by complementing the conventional approach with ML techniques. Image segmentation based on the k-means clustering algorithm is utilized to obtain the optical contrast of the TMD layers on a SiO2/Si substrate. To our knowledge, this is the first work on the layer number determination of WS2 samples on top of the common SiO2/Si substrate using optical images. A recent study identifying the number of layers of WS2, however, was demonstrated only for transparent glass substrate and the layers were not reliably distinguished using contrast [17]. Here, the inclusion of image segmentation allows one to identify the layers without having to perform a full contrast scan on the image. In addition to making the contrast calculation more robust, this approach reduces the time required.

We demonstrate the workflow on two of the common TMDs of interest WS2 and MoS2 showing its universality in application to TMDs. In addition, graphene and h-BN samples are also analysed using the workflow hinting that the method can be adopted for most 2DMs.

2. Background

Fresnel's multilayer interference theory can be used to understand the origin of contrast for 2DMs on top of oxidised silicon substrates [31]. An incident beam on the substrate undergoes reflection from the air-substrate interface as well as the SiO2–Si interface. The incident beam on the sample undergoes reflection from the air–2DM interface, the 2DM–SiO2 interface and the SiO2–Si interface, see figure 1(a). The wavelength dependent contrast for the 2DM–substrate combination under normal incidence can then be calculated as [21]:

Equation (1)

where $R_0 (\lambda) = |r_0(\lambda)|^2$ is the intensity of the reflected light from the air-substrate system and $R(\lambda) = |r(\lambda)|^2$ is that from the air–sample–substrate system, as shown in figure 1(a).

Figure 1.

Figure 1. (a) Schematic for multi-layer interference. The refraction at the interfaces is not depicted since the interest is on near-normal incidences. The silicon layer is considered to be semi-infinite. (b) Raman spectra of 1–4 layer WS2 regions and reference spectrum from the Si/SiO2 substrate under 532 nm excitation. (c) WS2 PL spectra up to four layers. (d) PL intensity normalized by the intensity of the monolayer peak at 1.94 eV as a function of the number of layers.

Standard image High-resolution image

Equation (2)

Equation (3)

Here, $r_{ij} = \frac{\tilde{n_i}-\tilde{n_j}}{\tilde{n_i}+\tilde{n_j}}$ with n0, $\tilde{n_1}$, $\tilde{n_2}$ and $\tilde{n_3}$ are the refractive indices of air, 2DM, SiO2 and Si respectively and the phase acquired by the beam while passing a distance d through the media with refractive index n is given by $\phi_{i} = \frac{2\pi n_i d_i}{\lambda}$.

Experimentally, the value of contrast in equation (1) can be obtained from pixel values of the optical image. The digital images captured using the charge coupled device detector in the microscope classifies the entire white light into three windows, red (590–720 nm), green (520–590 nm) and blue (435–520 nm) channels using the Bayer color filter. The reflected intensity from each of these wavelength ranges gets stored as the pixel values of R, G and B channels of the digital image. The contrast can be found from equation (4):

Equation (4)

where R0 is the pixel value corresponding to the substrate and R is that of the sample in the specific channel.

3. Methods

The WS2, MoS2 and Graphene flakes were mechanically exfoliated from bulk crystals using scotch tape [1] and were transferred to a 500 µm thick Si wafer with a 285 nm thermally oxidised layer (Si/SiO2). The white light images were captured using an Olympus microscope mounted with the Raman Spectrophotometer (Horiba) used for characterization. The image processing sequence and contrast calculations were performed on these white light images. The illumination intensity was fixed during capture by keeping an approximately constant substrate pixel intensity for all the samples since contrast strongly depends on the illumination intensity, see section 9 in the supplementary information. All the images were captured using a 100X objective and range from 1–2 megabytes in memory size. The dimensions were 2048 × 1536 pixels for the uncropped image. The images under physical color filters were taken using an Olympus microscope equipped with a differential interference contrast (DIC) filter. These images were only used for confirming the wavelength range in which the samples show higher layer-dependent visibility. The Raman and PL spectra were taken using a Horiba HR Evolution LabRam spectrophotometer using a 532 nm excitation laser. A 5 s acquisition time and 0.1% (∼10 µW) power are used for capturing the PL spectra. Laser power is kept minimum to avoid heating induced variations in the PL spectra. ImageJ was used to obtain the intensity profiles. All the image processing steps were performed using Python 3. The GUI was developed using PyQt5. The full version of the python script and sample images are available in GitHub repository (details in the supplementary information file).

4. Results and discussion

4.1. Tungsten disulphide (WS2)

Figures 1(b) and (c) show the Raman and the PL spectra of 1–4 layers of mechanically exfoliated WS2 flakes. E$_{2g}^1$ peak at ∼357 cm−1 and A$_{1g}$ peak at ∼418 cm−1 are the most prominent peaks in the Raman spectra [19]. As evident from the figure 1(b), the Raman spectrometry allows the unambiguous identification of the monolayer WS2 (E$_{2g}^1$ peak intensity greater than the A$_{1g}$ peak intensity). The identification of higher layer numbers is not that straight forward as the Raman spectra using 532 nm laser excitation does not contain any feature that significantly varies with an increase in the number of WS2 layers [19]. In contrast to the Raman spectra, the PL intensity sharply decreases when the thickness of the sample increases from monolayer to a few layers. The asymmetric monolayer PL peak (∼1.94 eV) starts to decrease in intensity as the layer number increases (see figure 1(c)). Hence the PL intensity of different layers normalized by the monolayer intensity can be used to identify the number of layers of the WS2 samples, see figure 1(d) [19]. Thus, the technique for layer number determination depends on the 2D material and it is essential to identify a method and process flow which can be used universally on all 2D materials.

Optical contrast is a universal feature as it clearly varies with the thickness of the sample irrespective of the material of interest, see equation (1). But the contrast of a supported 2DM is also dependent on the wavelength, the substrate material and its thickness. The method also suffers from its sensitivity to the image capturing conditions like the illumination intensity, the light reaching the detector from other sources leading to noisy pixels, slight variations in the sensitivity of individual pixels in the detector etc. Since there are a large number of parameters, instead of looking for a universal function that connects the contrast with the layer number, we will focus on developing a workflow that can be applied to any set of customized conditions. In this study, the substrate material and thickness are fixed by the choice of SiO2/Si wafer with 285 nm oxide layer. To rule out the variations in contrast due to different illumination of different samples, a fixed pixel intensity was maintained for the substrate for all the images used for processing. Figure 2 shows images of WS2 samples with different number of layers in different wavelengths realized using a microscope equipped with a DIC filter. Here the contrast clearly varies with wavelength and the red end of the visible spectrum gives better layer-dependent contrast. This can be attributed to the low optical absorption towards the red end of the visible spectrum since the bandgap of WS2 ranges from 1.3 eV in bulk to 2.03 eV in the monolayer limit [14]. In the Fresnel's multilayer interference model, the contrast is determined by the reflected rays from various interfaces and thus the less absorbed regions of the spectra, in general, show higher contrasts. We also captured optical images of WS2 samples under white light illumination (the preferred light source in optical microscopy). Then converted these images into red, green and blue color channels to see the contrast. The optical images and intensity profile of the WS2 flakes with different layer numbers in red, green and blue color channels of the white light image, given in figure 3, also confirm that the layer-dependent contrast of WS2 thin layers is more evident in the red end of the spectrum. This establishes the equivalence between color channels in the digital image and images under different wavelengths. This has the additional advantage of mitigating the problem of taking images under different wavelength filters. Even now, the intensity profiles along the line cuts in figures 3(d)–(f) show slight variations in intensity across those regions that appear to be uniform regions for our eyes. This is detrimental to contrast calculations since these are low contrast systems and the contrast is very sensitive to the pixel values.

Figure 2.

Figure 2. WS2 image containing (a)–(c) monolayer (d)–(f) bilayer (g)–(i) 3,4 layer samples near red, green and blue wavelengths, respectively.

Standard image High-resolution image
Figure 3.

Figure 3. WS2 images under different color channels in the digital image (a) red (b) green (c) blue with the corresponding intensity profiles (d)–(f) along the line cut.

Standard image High-resolution image

In order to eliminate these variations while calculating contrast, we have implemented an image segmentation based workflow that can be used to obtain a custom fit for contrast as a function of layer number for each laboratory. The number of layers of WS2 samples on SiO2/Si substrates are identified using this workflow. Figure 4(a) shows a representative image of a region containing monolayer and bilayer flakes of WS2. The first step is to crop the relevant region of interest with some part of the substrate in the background, see figure 4(b). This helps in reducing the computation time associated with the segmentation step that follows. Cropping also helps to avoid unwanted effects due to non-uniformities in illumination intensity. Large height variations between flakes in the same image can also lead to unfocussed regions in the image if not cropped. The cropped image is converted to different color channels and the red channel with the highest layer-dependent intensity variation is selected for further processes (see figure 4(c)). A bilateral filter is applied to this image to reduce the effect of noisy pixel values that are very different from the neighbouring pixels, see figure 4(d). The bilateral filter replaces each pixel value with a weighted average of the nearby pixel values, thereby preventing the stray pixels from forming individual clusters during segmentation. Since the filter also accounts for the intensity differences between pixels, the edges are better preserved compared to the Gaussian blur, another common filter used for noise reduction. The denoised image still contains small non-uniformities in pixel values for various regions. In order to have uniform pixel values across various regions in the image, we perform image segmentation on the denoised image. Image segmentation transforms the image to several distinct regions of uniform pixel values. We choose the k-means clustering algorithm for segmentation since k-means has the simplest implementation and is faster than other clustering algorithms like the Gaussian Mixture Model and Mean-shift [32]. k-means is an unsupervised Machine Learning algorithm commonly used for solving problems in computer vision. KMeans from the Sci-kit Learn library of python is used to implement the algorithm [33]. KMeans chooses the number of clusters, k as input and initializes with k random points as centroids. Then it assigns each pixel value to its closest centroid. The new centroid is found by calculating the median of the pixel values associated with each centroid. The difference between the new and the old centroid is calculated. This procedure is repeated until the new centroid is different from the previous centroid by less than some threshold, giving rise to k clusters in the dataset and k different regions in the image.

Figure 4.

Figure 4. (a) Optical image of WS2 flakes on Si/SiO2 substrate under white light illumination. (b) The cropped image containing mono-layer, bi-layer and few-layer regions with bare substrate. (c) The red channel image of the cropped region. (d) Red channel image after applying the Bilateral filter. (e) Elbow plot—WCSS as a function of the number of clusters, showing that optimal number of clusters is five. (f) Segmented image with five clusters. (g) Intensity profile along the line cut in (f).

Standard image High-resolution image

Choosing the number of clusters is then of prime importance. The Within Cluster Sum of Squares (WCSS or Inertia) for the image is a measure of the compactness of the clusters formed [34]. The Inertia with different numbers of clusters for the image, starting from 1 to 11, is calculated. The value of k which shows a significant reduction in WCSS starting from k = 1, and after which further increase in k shows only a small reduction in Inertia, is chosen to be the optimal number of clusters (see figure 4(e)). When the number of clusters crosses 5, there is no significant decrease in WCSS indicating that 5 is the optimal number of clusters for the image. However, the choice between k = 4 and k = 5 requires attention. At first glance, the image has only four different regions in it. A careful inspection shows that there are residues in the image that forms a fifth separate region in the color space. If this region is not considered as a separate cluster, the centroids of the regions of interest will shift, leading to fictitious contrast values. Therefore, in addition to minimizing WCSS, it is also important to choose a k that reproduces all the different regions of interest in the original image. More information about this can be found in figure S1 in the supplementary information. Figure 4(f) shows the image segmented into five clusters. Figure 4(g) shows the intensity profile of the clustered image along the line cut in which the pixel values are uniform across the different regions in the image. The pixel values corresponding to various regions of the sample and the substrate are extracted from the segmented image to find the contrast of each layer using equation (4).

With the layer number from the PL measurements and contrast from the segmented image, the contrasts is plotted against layer number and the data shows a nearly linear dependence within the range of measured layer numbers, see figure 5. Hence, the data is fitted with a linear function and the fitting parameters are recorded. After obtaining the fit parameters, freshly prepared samples with unknown layer numbers are tested and the method reliably determines the layer number. The test results are given in table 1. In order to test the number of layers of new samples, the images are pre-processed and segmented. The contrast is then calculated by picking the substrate and sample pixel values. This contrast is provided as an input to the obtained linear function to determine the layer number.

Figure 5.

Figure 5. Contrast of WS2 samples as a function of layer number fitted with a linear function.

Standard image High-resolution image

Table 1. Test results on random WS2 and MoS2 samples, the observed contrast is used with the fit in figures (5 and 9) to obtain the layer number. This is compared to the layer number obtained from characterization.

SampleObserved contrastLayer number from fitLayer number from fit (rounded)Layer number from characterization
NJ_WS_05_09_a0.11870.95 1 1
NJ_WS_05_10_a0.20912.4 2 2
NJ_WS_05_10_b0.36574.93 5 >4
NJ_WS_04_02_A0.2332.79 3 3
R_MoS2_06_010.27312.37 2 2
R_MoS2_06_080.13961.07 1 1
R_MoS2_06_170.13451.03 1 1
R_MoS2_06_190.14151.09 1 1

The complete workflow to obtain a custom fit for the image–substrate pair and its execution is summarized in figure 6. The workflow contains four stages. The first stage, 'Pre-processing', includes cropping the image to the desired region, converting it to the channel with the strongest layer- dependent contrast for the material and applying the bilateral filter. The second stage, 'Image Segmentation', involves finding the required number of clusters for the pre-processed image and segmenting the image using k-means. In the third stage, the contrast is extracted for several images, which were treated through the first two stages to obtain a fit for contrast as a function of layer number for the particular 2DM–substrate pair. In the final stage, the fit obtained is used to determine the number of layers from the contrast calculated for new images after subjecting them to pre-processing and segmentation. After the successful implementation of the ML-assisted optical contrast workflow on WS2 samples, we tested its applicability on other TMDs, graphene and h-BN.

Figure 6.

Figure 6. The four stages of the image processing workflow.

Standard image High-resolution image

4.2. Molybdenum disulphide (MoS2)

Unlike WS2, the layer number of MoS2 flakes can be unambiguously confirmed using the Raman spectra, see figure 7(a). The high-frequency Raman peaks around 385 cm−1 and 405 cm−1 are labelled as E$_{2g}^1$ and A$_{1g}$, respectively. As the layer number increases the peak position of E$_{2g}^1$ mode (Pos E$_{2g}^1$) decreases and that of A$_{1g}$ (Pos A$_{1g}$) mode increases [35]; see figure 7(b). Hence, the frequency difference between E$_{2g}^1$ and A$_{1g}$ peaks can be used to determine the layer number of MoS2. In addition to this, the Raman active shear modes (S) and layer breathing mode (LB) also show layer dependence in MoS2, see figures 7(a) and (d) [36]. These modes are absent in the monolayer. The frequency of the shear mode increases as the layer number increases from bilayer to bulk. We utilised both the low and high frequency modes of MoS2 to determine the layer number accurately.

Figure 7.

Figure 7. (a) MoS2 Raman spectra (532 nm) up to four layers. (b) Peak position of the E$_{2g}^1$ and A$_{1g}$ modes as a function of the number of layers. (c) Difference between the peak positions of the E$_{2g}^1$ and A$_{1g}$ modes as a function of the number of layers. (d) Raman peak position of the low energy shear mode as a function of layer number.

Standard image High-resolution image

The optical images and intensity profile of MoS2 samples under different physical filters and samples under different color channels are shown in supplementary information figures S2 and S3. The layer dependence of intensity is more prominent towards the red end of the spectrum. This is in agreement with the calculations in previous reports claiming that MoS2 shows the highest contrast in the red wavelengths [21]. Figure 8 shows the MoS2 images taken through the same workflow that was implemented on WS2 samples. The contrast observed fitted well with a second-degree polynomial as a function of layer number, see figure 9.

Figure 8.

Figure 8. (a) MoS2 image with 2,3,4 layer regions. (b) Image cropped to the region of interest. (c) The red channel of the cropped image. (d) Red channel image after applying the Bilateral filter. (e) Elbow plot—WCSS as a function of the number of clusters. (f) Segmented image with five clusters. (g) Intensity profile along the line cut in (f).

Standard image High-resolution image
Figure 9.

Figure 9. The contrast of MoS2 samples as a function of layer number fitted with a second-degree polynomial.

Standard image High-resolution image

Randomly chosen new MoS2 samples were used to validate the obtained fits and the results are given in table 1. The layer number determined using the above fits, figures (5 and 9) is compared to the layer number determined from standard characterization methods (Raman and PL). The results clearly show that the workflow can withstand variations and can be adopted for characterizing the number of layers reliably for different two-dimensional TMDs.

Finally, we also demonstrate that the workflow can be extended to general two-dimensional materials by applying it on graphene and h-BN samples, see supplementary information figure S4, S5, S9 and S10. The graphene samples were mechanically exfoliated and upto three layers were confirmed using Raman spectra. The intensity ratio of the G and 2D peak were used for identifying the layer number of the samples [37]. The images of graphene samples were taken under different physical filters and the layer-dependent contrast is found to be a maximum in the green channel in agreement with previous reports [21]. The segmentation was performed on the green channel image and a linear fit is obtained for contrast as a function of layer number for the samples tested. Thus the workflow devised is robust and can be adopted universally to determine the layer number of 2DMs 1 .

In order to make this workflow easily accessible, a GUI is developed using PyQt5. Cropped images can be loaded to the interface and all the subsequent steps in the workflow can be accessed using various push buttons. The number of layers of WS2 is displayed as an output. Currently, the red channel is kept active and the fit parameters of contrast as a function of layer number for WS2 are used in the GUI. However, the code can be easily customized for the desired material, colour channel and corresponding fit parameters. Figure S6 shows the interface of the GUI titled 'TMD Layer Number Calculator'. More information about using the GUI can be found in the supplementary information file.

5. Conclusion

We have shown that the number of layers of WS2 flakes exfoliated on top of the SiO2/Si substrate can be determined from the optical images. Instead of training a model, we have implemented an image processing workflow that utilizes the k-means clustering algorithm to partition the images to separate homogeneous regions, which makes the contrast calculation more reliable. Our work shows that the inclusion of ML-based algorithms like k-means without a full-fledged adoption of the conventional Machine Learning workflow can improve optical layer number determination without the requirement for a large training data set. We have also demonstrated that this workflow can be utilized for obtaining correlations between optical contrast and the number of layers for MoS2, graphene and h-BN, showing its generality. A GUI program titled 'TMD Layer Number Calculator' is provided to make the workflow more accessible to the research community. We believe that this method will enable researchers identify the layer number of different two dimensional materials reliably and will also complement the existing experimental techniques. Since the method relies on optical contrast of the images, it is compatible with automatic sample detection systems and can be implemented with large scale 2DMs production systems in the future.

Data availability statement

The data that support the findings of this study are available within the article and its supplementary material, and in https://github.com/NikhilJosephJoy/Layer-number-identification-of-TMDs [38].

Acknowledgments

J B acknowledges support under DST-SERB Grant CRG/2020/000615. The authors thank the Central Instrumentation Facility (CIF) and Central Micro-Nano Fabrication Facility (CMFF), Indian Institute of Technology Palakkad for the experimental facilities.

Footnotes

  • Depending on the material and the wavelength, the contrast may have a linear or nonlinear dependence with layer number. But atleast up-to ten layers the fit parameters alone can predict the layer number accurately. For higher number of layers similar contrast may appear for different layer numbers. In this case we may have to look for confocal microscopy images or AFM or other methods for height profile along with optical contrast to accurately predict the layer number.

Please wait… references are loading.

Supplementary data (4.8 MB PDF)