**Towards the Development of Rapid and Low-Cost Pathogen Detection Systems Using Microfluidic Technology and Optical Image Processing**

#### **Abdelfateh Kerrouche 1,\*, Jordan Lithgow 1, Ilyas Muhammad <sup>1</sup> and Imed Romdhani <sup>2</sup>**


Received: 29 February 2020; Accepted: 1 April 2020; Published: 7 April 2020

**Abstract:** Waterborne pathogens affect all waters globally and proceed to be an ongoing concern. Previous methods for detection of pathogens consist of a high test time and a high sample consumption, but they are very expensive and require specialist operators. This study aims to develop a monitoring system capable of identifying waterborne pathogens with particular characteristics using a microfluidic device, optical imaging and a classification algorithm to provide low-cost and portable solutions. This paper investigates the detection of small size microbeads (1–5 μm) from a measured water sample by using a cost-effective microscopic camera and computational algorithms. Results provide areas of opportunities to decrease sample consumption, reduce testing time and minimize the use of expensive equipment.

**Keywords:** pathogen detection; microfluidics; image processing; computational algorithms

#### **1. Introduction**

#### *1.1. Water Borne Pathogen*

Due to the increasing concern for water scarcity, it is expected that by 2025, half of the world's population will endure 'water-stressed' areas [1]. This is due to climate change affecting both the developing and developed world. Already, 1.1 billion people lack access to clean water and 2.7 billion have a water shortage for at least one month of the year, exposing waterborne diseases such as cholera, typhoid fever and hepatitis A [2]. As a result, approximately 525,000 children under the age of five die annually from diarrhea diseases alone. According to the World Health Organization (WHO), lower respiratory infections and diarrheal diseases are the top two causes of deaths in low-income countries in 2016. From this, 11–20 million people get an illness due to typhoid fever, and 128,000–161,000 people die from this annually [2]. Typhoid fever is caused by the Enterobacteriaceae Salmonella Typhimurium, which typically is around 2–5 μm in size. In addition to this, the developed world suffers from pathogens such as Cryptosporidium, which can survive months within waters [3,4]. This particular pathogen causes gastrointestinal and respiratory illnesses, although the death rate is much lower than typhoid fever. Detection protocols such as the U.S. Environmental Protection Agency (EPA) method 1623.1 [5] stipulate a procedure for *Cryptosporidium* detection. This method consists of several steps involving filtration, elution, centrifugation, immuno-magnetic separation (IMS) and staining with fluorescent dyes, followed by microscopic examination for the manual identification and enumeration of oocysts.

Bridle [6] has previously established that research has been conducted to identify advances in obtaining results of waterborne pathogen movements. In compliance with WHO, waterborne pathogens have to be carefully monitored and analyzed to implement the safety of water intake. Considering this, microfluidic devices (a set of micro-channels etched into a material e.g., glass [7]) can take a very small sample and conduct a process to analyses the safety of water through an experimental monitoring method.

#### *1.2. Microfluidic Technology*

Newly introduced technology is continuously being miniaturized, hence the attractiveness of microfluidic devices. Microfluidic technology was introduced around 20 years ago and is maturing as a technology, despite the initial lack of success on the market [8,9]. The device was influenced from microelectronic chips; however, the use and analysis of fluid movement was changed from the use of pathways for electrons. Microfluidic devices accommodate more mixers and separators rather than transistors and other microelectronic hardware [10]. This increases microfluidics applications to be developed to perform in applications, such as drug testing in nanosystems, biochemistry immunoassay and gene sequencing etc., to proceed at a faster, cheaper rate, whilst decreasing the sample and reagent consumptions [11]. From these benefits, other applications use the devices to provide aid in their success, such as cosmetic formulations, biology cell culture and energy plasma confinement [12]. The science behind microfluidic technology is to manipulate and control small levels of fluid from microliters (μL) to picoliters (pL) within small micrometer channels, allowing for the precision of experiments to be maximized and a higher probability of detection to be available. For the purpose of this paper, the main focus will be to detect and monitor 'acting' pathogens via a smart designed system. In the past few years, microfluidic systems have been intensely researched to advance applications in the fields of chemistry, biology, genomics, proteomics, pharmaceuticals, biodefense [13] emerging as a powerful useful technology. Large processes that once were carried out in bulky, slow, expensive equipment are now processed on a single miniaturized chip. Microfluidic chips are purposely designed to monitor the movement of cells, pathogens or nanoparticles within the microchannels, allowing easy detection and removal methods under controlled conditions. In addition, microchannels typically have a width of 100–800 μm which allows advantages of sample consumption to be reduced, whilst also increasing the speed of results. This allows increasing efficiency and an inexpensive process [12,13].

#### *1.3. Microfluidic Flow Cytometry*

Microfluidic flow cytometers techniques are one of the most powerful approaches for a high resolution of cell analysis. This technique was developed from a traditional flow cytometer, which aids in biological applications and clinic research [14]. Flow cytometers are popular due to the ability to sort, count and detect individual cells, while also measuring scientific characteristics and the handling of cell populations [14,15]. However, flow cytometers do not effectively detect particles and cells smaller than 0.5 μm in diameter via light scattering [16]. Light scattering manipulates singles out pathogens based on size and shape. Microfluidic flow cytometry technique aims to process a sample by 3 stages: external pumping force of the sample, focusing of particles within the sample and sorting of particles by separation. This allows particles to be manipulated and analyzed faster, whilst being less expensive than a conventional flow cytometer [17,18]. The limitations of a microfluidic flow cytometer are that, due to a small flow rate of the external pump force, it will consist of a non-steady flow, which contributes to installing extra fabrication steps and limited range of flow rate, increasing the production cost [19].

In addition, optical detection within activated flow cytometry for particle detection is well advanced for microfluidic techniques. For instance, fluorescent dye is the most common optical detection technique due to its high sensitivity, which illuminates targets within samples under a specific wavelength of a light source [6]. This allows reflection of the particle to be transmitted and detected from the optical point of a microscope. The development and integration of different controlling units (e.g., pumps, valves and switches) into microfluidic devices are challenging for an automated process. This is due to the response time needed to redirect the particles into the desired output channels [18]. New advances for optical detection and the separation of cells or particles in the modern flow cytometer defeat these drawbacks and reduce the bulkiness of the system. Integrating optical fibers and waveguides to within the device around the microfluidic channels allows results to be optimized and analyzed without the use of a microscope or camera, as Lin et al. [19] proved with the detection of Phalaenopsis pathogens using optical fibres.

#### **2. Sample Preparation**

#### *2.1. Methodology*

To monitor the movement of a waterborne pathogen, a micro-pump (The Aladdin single-syringe infusion) drives a sample of 3 mL through a microfluidic device. To simulate the movement of a pathogen inside the microfluidic channel, a polymer fluorescent green microsphere beads are used with a range size of 1–5 μm.

A USB microscopic camera (The Dino-Lite AM4515T8 Edge 700x~900) was set up and positioned at a perspective angle. The camera emits an ultraviolet light source, causing the fluorescent microspheres to reflect to the camera. As shown in Figure 1a, the camera was secured in a position using a clamp to obtain a perspective view of the microfluidic device. The camera's focus point was altered to output a clear view of the microfluidic channels on the software. Five hundred micrograms of polymer fluorescent green microspheres were diluted and mixed into a 50 mL test tube containing clean fluid (Propane-2-Ol, "IPA"). This allowed the sample to be withdrawn easily with the use of a BD Plastipak 3 mL max syringe. The tubing was connected to the input port of the microfluidic device. Tubes were connected to the output ports of the microfluidic device to release the wastage processed. These green polymer beads are used in this paper because they have the same size and the same color as the waterborne pathogen: Cryptosporidium. Figure 1b shows an image of a microbead inside the microfluidic channel.

**Figure 1.** (**a**) The system setup includes a syringe pump and a USB camera with 500–900 times zoom; (**b**) image detection of the microbeads of 1–5 μm within the microfluidic channel.

The tests were performed in the lab with controlled ventilation. The temperature was set to an ambient room temperature of 22 degrees Celsius, and the humidity was measured for an indoor level of 43%. These two parameters of ambient temperature and relative humidity can affect optical image quality taken by the camera and can also affect the fluid movement inside the channel, as the microfluidic is made from a clear plastic.

#### *2.2. Velocity Tests*

Velocity tests of microbeads have been carried out in this research to study the movement of beads inside the microfluidic channel. The micro-pump injected the sample at a controlled flow rate through the microchannel of a custom designed microfluidic chip via capillary tubes. Velocity measurements will help to identify the limitation of this type of microscopic camera in terms of the minimum frames per second needed for this application with a selective flow rate. This is important, as lower frame rates can result in choppy or a broken movement.

Flow rates were tested from 5–235 μL/min (in increments of 10 μL/min) to find the most probable flow paths for a pathogen cell. Videos from the camera stored the movement of the fluorescent bead within the microchannel at a rated flow rate in increments of 30 μL/min. A highlighted bead expresses the motion of which it travels at a specific flow rate. Flow rates between 5–35 μL/min resulted in no movement of particles therefore were not included. The movement is processed from the images using Matlab, allowing analysis of where most beads flow in terms of the 0.8 mm diameter channel, as shown in Figure 2.

**Figure 2.** Detection of microsphere movement at 45 μL/min.

In order to find the velocity of the fluorescent beads at different flow rates, a distance is fixed at 250 mm for all flow rates tested using the software provided by the camera, as shown in Figure 3, The general equation for finding speed is Equation (1), where v is the velocity, d is the distance and t is the time. The time is found using the video recorded by the Dino Lite camera—this is also highlighted in Figure 3.

$$\mathbf{v} = \mathbf{d}/\mathbf{t} \tag{1}$$

**Figure 3.** Dino Lite 2.0 software with evidence of fixed distance of 250 mm and a time interval of when microsphere travels the distance.

An average time is recorded from multiple tests for microspheres travelling in 3 sections and these are inside, middle and outside, as described in Figure 4.

**Figure 4.** Drawing illustrates the three sections of the microfluidic channel.

#### **3. Results**

#### *3.1. Fluorescent Microspheres Detection*

The main objective of this study is to monitor the movement of waterborne pathogen cells using a microscopic camera and a light source integrated into a microfluidic device. The micro-pump is configured at a very low flow rate of 45 μL/min, to investigate the possibility of detection. After the operating equipment is correctly configured, 500 μg of a solid bead powder is diluted into a 50 mL test tube containing clean fluid (Propane-2-Ol, "IPA"). This master sample is mixed using ultrasonic and steel ball merging to gain a homogeneous bead distribution. This allowed samples to be withdrawn easily, by extracting 5 μL from the master sample and added it into 10 mL and 100 mL of IPA. The first images are obtained using the method in Section 2.1 to produce a high reflection rate of fluorescent microspheres. Therefore, a diluting process is performed to receive more realistic images.

Figure 5a demonstrates the mass of fluorescent green beads within the master sample. Figure 5b shows image with diluted sample, which contains a less dense concentration of fluorescent beads. Therefore, a further dilution process is performed to result the image in Figure 5c, where only a very small concentration of beads are detected, expressing a more realistic result of pathogen movement.

**Figure 5.** Diluting process images illustrated by (**a**–**c**)**.**

#### *3.2. Image Processing Using Computational Algorithm*

The microsphere beads used in this research are inert, nontoxic and do not contaminate the device, so there is no requirement to move it to a dedicated biology laboratory. However, the manufacturer (Cospheric) does not specify the size distribution of the beads, hence the number of beads per milligram can vary (calculated with the given diameter range) between 1.2 <sup>×</sup> <sup>10</sup><sup>7</sup> and 1.5 <sup>×</sup> <sup>10</sup><sup>9</sup> (using ideal spheres and the material density 1.3 g/cc). This makes a measurement by volume and weight not suitable.

Different computational methods and algorithms for biological microscopic image analysis have demonstrated a comprehensive contribution for efficient and quantitative study of microorganism cells [20,21]. Experiments have been conducted in this paper to detect and quantify fluorescent microbeads (1–5 μm) in a measured sample of a clean water. Detailed steps are outlined as: (i) taking 2 μL of sample via micropipette; (ii) placing one drop on top of a microscopic glass slide; (iii) captureing an image of the full drop under microscopic camera; (iv) capturing zoomed images of the same drop frame by frame and finally (v) analyzing captured images by using computational algorithms using Matlab. The Circular Hough Transform (CHT) algorithm [22,23] based code in Matlab, is used to find circular/spherical shapes in the captured image as one frame of a drop, as discussed above. The CHT algorithm is used to mark/circle specific sized microbeads for quantification. CHT algorithm accuracy rate is affected by different parameters/factors, such as minimum and maximum radius of the target circular shape, object polarity in the sense of dark or bright, and computation method. The results using CHT algorithm are approximately 80%, as reported in Table 1.


**Table 1.** Values of different parameters used during different iterations.

Figure 6 shows resulting images of different iterations, which refer to the number of runs of the CHT algorithm applied in each image capture. First, the image is converted to grayscale, as color images are not giving a good segmentation. Then, minimum and maximum radius are defined as 1 and 600 pixels respectively. Finally, the object polarity is set to dark.

**Figure 6.** Original image and resulting images of different iterations, as shown in Table 1.

#### *3.3. Microsphere Detection with Change of Flow Rate*

Table 2 illustrates a sample of data collected in terms of number of beads over a fixed period of time (30 s). This time is fixed by Matlab to count the number of the total beads passed within the microfluidic channel, as shown in Figure 3. The data shown in Table 2 demonstrates an overview of the highest probable chances microspheres are going to travel inside the channel at specific flow rates. The evidence of the movement detection allows an analysis of where most beads flow in terms of the 0.8 mm diameter channel displayed in Table 2. The microfluidic channel diameter is then divided into two sections on Table 2 (outside to middle, middle to inside) to demonstrate where the bead is more likely to travel.


**Table 2.** Number of beads counted after 30 s. Terms of specific flow rates against specific areas of microfluidic channels at the same time for each test.

As it could be seen in Table 2, between the outside of the meander 0.8 mm to the middle sector at 0.4 mm, there is a higher probability that the bead will travel between the flow rates of 45–105 μL/min.

Figure 7 shows results from the data collected over three tests for each flow rate to determine the recovery rate, which is the number of beads detected in the region in question over the first 100 beads detected within the channel. As demonstrated in Figure 4, 69% of beads can be recovered from the sample between the outside of the meander 0.8 mm to the middle sector at 0.4 mm. On the other side, from the middle of the channel at 0.4 mm to the inside 0.0 mm the probability that the bead will flow between these points is at the flow rate between 135–225 μL/min. This provides a recovery rate of 88%, conveying that the increase of flow rate increases the consistency of where the microsphere will travel. However, this will not allow an easy detection, as the bead's velocity will increase too, as explained in Section 3.4.

**Figure 7.** Recovery rates within the microfluidic channel.

From these results, fabrication of a microchannel can be inserted for the pathogen to be removed from the sample. Being able to do this will allow microfluidic technology to essentially replace most filtration methods used to remove waterborne pathogen such as EPA *Cryptosporidium* detection method [5].

#### *3.4. Velocity Measurements*

As this method detects the movement of a fluorescent bead, it allowed further testing to find a single flow rate, where a pathogen cell will consistently flow on the same path. This provides an opportunity to locate the pathogen within the sample, which helps with placing the camera at the right position for detection. The microfluidic channel diameter is divided into three sections to demonstrate where the bead is more likely to travel, as shown in Figure 4.

Figure 8 illustrates the average velocity of the microsphere travelling at the representing flow rate within the microchannel. Results demonstrate that the gradual increase of flow rate up to 145 μL/min slightly increases the velocity of the microsphere. However, after this flow rate, the velocity of the beads increased significantly in all three parts of the microfluidic channel. This is due to the extremely high velocity of the microsphere that is unable to be detected by the camera—and only some 'slower' microspheres could be detected as altering the average velocity.

**Figure 8.** Gradual increase of velocity for the microsphere travelling within the microfluidic channel with respect of flow rate increase.

#### **4. Discussion and Conclusions**

This research aims to design a smart monitoring system to detect a waterborne pathogen cell using microfluidic technology and a cost-effective microscopic camera. The sample was driven through the designed microfluidic channel at a selected flow rate on a micro-pump. The pathogen movement was tracked by a high-quality camera at a perspective angle. With repeated experiments and alterations of the input flow rate, the output is hypothesized to find a consistent 'track line' that the pathogen travels within the channel, so a separation channel can be fabricated to separate the pathogen from the sample.

The prominent challenges for the pathogen detection methods in practice are: (i) physical characteristics of different groups of pathogens, (ii) large volume of water containing low concentration of target pathogens and (iii) sample preparation protocols. The key advantages for the results delivered by this new method allows sample consumption to be decreased, the reduction of testing time and reducing the use of expensive equipment, providing an attractive cost. This will benefit existing methods such as Flow Cytometry and EPA Cryptosporidium detection approach by replacing procedures with a less complicated process that does not need specialized persons to operate. This method provides new possibilities for further analysis and research to improve results on finding precise intersection points for high recovery rates of waterborne pathogens with different sizes such as Giardia or *Escherichia coli* (*E. coli*). This provides vital intersection points for a microchannel to be fabricated to remove the waterborne pathogen.

A key objective for this paper is to investigate the design of an intersection point on the microchannel to separate the waterborne pathogen. To achieve this, the micro pump flow rates were changed slightly to find a consistent path of which the microsphere will travel. Figure 3 displays the movement of microspheres around a meander. This gives a general observation that the slower the flow rate is (45 μL/min), the microsphere movement is towards the outside of the microchannel. Comparing this to 225 μL/min flow rate, it portrays that the movement of the microsphere is on the inside of the microchannel. From the videos taken by the camera, Table 2 gathered data from the movement of the microsphere. This shows that, at slower flow rates of 45–105 μL/min, the microsphere has a higher probability of 69% to be recovered at 0.4–0.8 mm of the microchannel. Whereas at faster flow rates of 135– 225 μL/min, the recovery rate has a probability of 88% at 0.0–0.4 mm of the microchannel. However, as stated in Section 3.4, due to the fast velocities of microspheres at these flow rates, the detection became more difficult. This could slightly alter the probability rate of recovery, due to the possibility of some microspheres not been detected.

The detection accuracy of the beads with the CHT algorithm used in this paper can be improved using different platforms such as ImageJ or Python. The processed image of a full drop can be performed using sequential steps to study different parameters such as target object size, the distance between two objects and threshold.

Microfluidic chips can be manufactured using different types of materials such as PDMS (PolyDiMethylSiloxane), Silicon, Glass or PMMA (PolyMethylMethAcrylate). These materials can be investigated to improve the images taken by the camera in order to improve recovery rates. Alongside this, the system described in this paper can be reduced even more into a single automated unit, to be able to optimize results that requires no manual intervention.

From all factors discussed in this paper, the idea of detecting movement of a waterborne pathogen using a cost-effective camera is very promising. This benefits engineering problems by minimizing the need for expensive laboratories and equipment to commence procedures on tackling waterborne pathogens.

**Author Contributions:** Conceptualization, A.K. and J.L.; methodology, A.K.; software, A.K. and I.M.; validation, A.K. and I.R.; formal analysis, A.K.; investigation, A.K; resources, A.K.; data curation, A.K.; writing—original draft preparation, A.K.; writing—review and editing, A.K.; visualization, A.K.; supervision, A.K.; project administration, A.K.; funding acquisition, A.K. and J.L. Please turn to the CRediT taxonomy for the term explanation. Authorship must be limited to those who have contributed substantially to the work reported. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors would like to acknowledge the support from Helen Bridle, Heriot Watt University.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

### *Article* **Dependence-Analysis-Based Data-Refinement in Optical Scatterometry for Fast Nanostructure Reconstruction**

**Zhengqiong Dong 1, Xiuguo Chen 2,\*, Xuanze Wang 1, Yating Shi 2,\*, Hao Jiang <sup>2</sup> and Shiyuan Liu <sup>2</sup>**


Received: 16 August 2019; Accepted: 27 September 2019; Published: 30 September 2019

**Abstract:** Optical scatterometry is known as a powerful tool for nanostructure reconstruction due to its advantages of being non-contact, non-destructive, low cost, and easy to integrate. As a typical model-based method, it usually makes use of abundant measured data for structural profile reconstruction, on the other hand, too much redundant information significantly degrades the efficiency in profile reconstruction. We propose a method based on dependence analysis to identify and then eliminate the measurement configurations with redundant information. Our experiments demonstrated the capability of the proposed method in an optimized selection of a subset of measurement wavelengths that contained sufficient information for profile reconstruction and strikingly improved the profile reconstruction efficiency without sacrificing accuracy, compared with the primitive approach, by making use of the whole spectrum.

**Keywords:** optical scatterometry; inverse problem; profile reconstruction; dependence analysis; data refinement

#### **1. Introduction**

Nano-metrology is the only effective method that ensures the reliability and consistency of nano-manufacturing [1–3]. Compared with other techniques such as scanning electron microscopy (SEM), atomic force microscopy (AFM) [4], and near-field scanning optical microscope (NSOM) [5], optical scatterometry [6,7], also known as optical critical dimension metrology or optical critical dimension (OCD) metrology, is more suitable for monitoring, assessing, and optimizing the nano-manufacturing processes due to its advantages of being non-contact, non-destructive, low in cost, and easy to integrate, etc. Recently, optical scatterometry has been applied in many fields with great success, such as the process control for back-end-of-the-line (BEOL) [8], the in-chip critical dimension (CD), overlay metrology [9], and in-situ measurement of pattern reflow in nanoimprinting [10].

In general, optical scatterometry involves two procedures: the forward optical modeling of sub-wavelength structures and the reconstruction of structural profiles from the measured signatures [11]. Here, the general term signatures means the scattered light information from the diffractive grating structure, which can be in the form of reflectance, ellipsometric angles, Stokes vector elements, or Mueller matrix elements. The forward optical model describes the light-nanostructure interaction by solving the complex Maxwell's equations. There are many reliable forward-modeling techniques such as the rigorous coupled-wave analysis (RCWA), the finite element method (FEM), the boundary element method (BEM), or the finite-difference time-domain (FDTD) method. The profile

reconstruction process conducted under a fixed measurement configuration is an inverse problem with the objective of optimizing a set of floating profile parameters (e.g., CD, sidewall angle, and height) whose theoretical signatures can best match the measured ones through regression analysis or library search [12–14]. The measurement configuration is defined as a combination of specially selected wavelengths and incident and azimuthal angles.

As a typical model-based method, abundant measured data is usually collected under the measurement configurations with fixed incidence and azimuthal angles, but a wide waveband for spectroscopic scatterometry, or with fixed measurement wavelength and azimuthal angle, but incidence angular range for angle-resolved scatterometry. Taking the inverse structural profile reconstruction in spectroscopic scatterometry as an example, the measurement sensitivities of unknown profile parameters are affected by the fixed incident and azimuthal angles, which can be improved by traversing their possible values; the size of fitting data usually directly depends on the number of selected wavelength points within the waveband. A large sized measured data set should be used to guarantee the fidelity of the fitting results, but on the other hand, too much redundant information will significantly degrade the efficiency during the inverse problem solving process due to the frequent iteration of the forward model. Moreover, the employment of the redundant information, which is insensitive to the measurands in solving inverse problems, leads to the loss in precision as well. Thus, an appropriate measurement configuration with high sensitivity and a small amount of fitting data will benefit the measurement with both precision and speed.

Many instrument optimization and spectrum denoising methods have been proposed with the objective of improving the quality of measured data [15–17]. As the measurement precision mainly depends on the sensitivity of model parameters, several approaches based on local sensitivity analysis (LSA) also have been developed to determine optimal incident and azimuthal angles for spectroscopic scatterometry [18–21]. We have previously provided an optimization method based on global sensitivity analysis (GSA) to determine an optimal combination of the fixed incident and azimuthal angles corresponding to the best measurement precision [22]. However, the present measurement-configuration optimization methods only focus on the improvement of measurement precision, the speed of the reconstruction process has not been considered. Usually, a vector consisting of several data points at multiple measurement wavelengths or incidence angles is calculated by a forward model to match the measured spectrum. The forward model's computing time linearly depends on the number of selected data points and it will be called frequently, therefore, it is highly desirable to find an appropriate strategy that can select the measurement configurations containing sufficient information for structural parameter reconstruction to enhance the measurement speed without sacrificing the accuracies, especially for the nanostructures whose forward model is very complex and time-consuming.

In this paper, we propose to refine the measured signatures by identifying and removing the measurement configurations with redundant information based on dependence analysis before the reconstruction process. As for the mathematical inverse problem, Twomey's pioneer work on atmospheric measurements has demonstrated that in a chosen set of data the most linearly independent measurements contain all the information needed for an inverse problem solution [23–25]. Then, Assaad estimated additional reflectance values using the several acquired independent measurements to obtain a larger set of measurement data in optical scatterometry, which is desired to reduce the uncertainty range of the parameter estimates [26]. Inspired by these analyses, we conducted information content analysis of the measured signatures based on dependence-analysis theory for optical scatterometry. The analysis was conducted by the eigen-decomposition of the Jacobian matrix **JJ**<sup>T</sup> in the linear model to allocate the most independent measurements in the acquired data set. The few remaining independent measurements were then adopted to reduce the reconstruction time, which is crucial for real-time and in-process applications.

It should be noted that this paper does not intend to discuss the existence and uniqueness of the inverse problem in optical scatterometry with respect to measurement conditions, which is a tough and still open issue in mathematics [27–29]. The main purpose of this paper is to provide an effective approach to choose a subset of measurement configurations from the whole optical wavelength spectrum or angular range that contains sufficient information for efficient structural profile reconstruction without sacrificing accuracy. The selected subset of measurement configurations may not be unique; its size may not be the smallest, but it should provide a higher reconstruction accuracy than its randomly selected counterparts and a higher reconstruction efficiency than the primitive approach that analyzes the whole optical wavelength spectrum or angular range.

#### **2. Method**

The inverse problem in optical scatterometry is described as an objective to minimize a least-square function, which can be generally expressed as:

$$\chi^2 = \sum\_{i=1}^{N} \left\{ \left[ \mathbf{y}\_i - \mathbf{f}(\mathbf{x}\_i \ \mathbf{a}\_i) \right]^T \mathbf{w}\_i \left[ \mathbf{y}\_i - \mathbf{f}(\mathbf{x}\_i \ \mathbf{a}\_i) \right] \right\}, \tag{1}$$

where, **x** = [*x*1, *x*2, ... , *xm*] <sup>T</sup> represents a column vector of *m* profile parameters (e.g., CD, sidewall angle, and height) and **f**(**x**, **a***i*) is a column vector containing the calculated signature according to the forward optical model under the *i*th measurement configuration **a***i*. The measurement configuration **a***<sup>i</sup>* is a combination of fixed azimuthal angle ϕ, fixed incident angle θ, and the *i*th measurement wavelength λ*<sup>i</sup>* for spectroscopic scatterometry, namely **a***<sup>i</sup>* = [ϕ, θ, λ*i*] T, and is a combination of fixed azimuthal angle ϕ, fixed measurement wavelength λ, and the *i*th incident angle θ*<sup>i</sup>* for angle-resolved scatterometry, namely, **a***<sup>i</sup>* = [ϕ, θ*i*, λ] T. Correspondingly, the *N* in Equation (1) denotes the number of measurement configurations, which can be the number of measurement wavelengths for spectroscopic scatterometry or the number of incidence angles for angle-resolved scatterometry. The column vector **y***<sup>i</sup>* = [*yi*1, *yi*2, ... , *yil*] <sup>T</sup> consists of the corresponding measured signature under the same measurement configuration **a***i*. Here, the value of *l* depends on the specific types of the measured signature. If the measured signature is reflectance, *l* = 1, while if the measured signature is a Mueller matrix, *l* = 16 (or *l* = 15 for the normalized Mueller matrix). The weighting matrix **w***<sup>i</sup>* is a *l* × *l* positive definite matrix, which is usually chosen to be the inverse of the covariance matrix of the measured signature under the *i*th measurement configuration **a***i*. In this case, Equation (1) relates to the commonly used multivariate chi-square statistics χ2. Writing the right side of Equation (1) in a matrix expression, the solution **x**ˆ of the inverse problem can be obtained by:

$$\hat{\mathbf{x}} = \underset{\mathbf{x} \in \Omega}{\text{arg min}} \Big[ \left[ \mathbf{y} - \mathbf{f}(\mathbf{x}, \mathbf{a}) \right]^T \mathbf{w} \big[ \mathbf{y} - \mathbf{f}(\mathbf{x}, \mathbf{a}) \big] \Big] \tag{2}$$

here, Ω is the associated profile parameter domain.

The use of χ<sup>2</sup> is built on top of the belief that the measurement errors are normally distributed with zero mean, namely, the measured signature **y** can be expressed as:

$$\mathbf{y} = \mathbf{f}(\mathbf{\hat{x}}, \mathbf{a}) + \mathbf{e}\_{\prime} \tag{3}$$

where ε is a vector of multiple random and independent variables with each element subjected to a normal distribution with zero mean. Assuming a profile parameter vector **x** is close enough to the true or nominal parameter vector **x**\* under a measurement configuration **a***i*, and the function **f**(**x**, **a***i*) is sufficiently smooth, then the function value **f**(**x**, **a***i*) can be expanded in the vicinity of **x**\* using the first-order Taylor expansion formulation [30]:

$$\mathbf{f}(\mathbf{x}\_{\prime}|\mathbf{a}\_{i}) \approx \mathbf{f}(\mathbf{x}\*\_{\prime}\mathbf{a}\_{i}) + \mathbf{J}(\mathbf{x}\*\_{\prime}\mathbf{a}\_{i}) \cdot \Delta \mathbf{x},\tag{4}$$

where **J**(**x**\*, **a***i*) is the Jacobian matrix with respect to **x** at **x** = **x\*** and Δ**x** represents the error in **x** and is given by Δ**x** = **x** − **x**\*. The differences between the measured signatures and the expected values can be expressed as:

$$
\Delta \mathbf{y}\_i \approx \mathbf{J}(\mathbf{x} \star \mathbf{a}\_i) \cdot \Delta \mathbf{x} + \varepsilon\_i. \tag{5}
$$

To determine the degree of independence among *N* measurements, Equation (5) is multiplied by an arbitrary factor η*<sup>i</sup>* and is summed over all *i*, to give:

$$\sum\_{i=1}^{N} \eta\_i \Delta \mathbf{y}\_i = \sum\_{i=1}^{N} \eta\_i \mathbf{J}(\mathbf{x} \star\_i \mathbf{a}\_i) \cdot \Delta \mathbf{x} + \sum\_{i=1}^{N} \eta\_i \varepsilon\_{i\prime} \tag{6}$$

where, *N i*=1 η2 *<sup>i</sup>* = 1. Then, the difference Δ**y***<sup>j</sup>* for a certain measurement configuration *<sup>j</sup>* can be expressed as [23]:

$$\Delta \mathbf{y}\_j = -\eta\_j^{-1} \sum\_{i \neq j} \eta\_i \Delta \mathbf{y}\_i + \eta\_j^{-1} \left\{ \sum\_{i=1}^N \eta\_i \mathbf{J}(\mathbf{x} \star\_i \mathbf{a}\_i) \cdot \Delta \mathbf{x} + \sum\_{i=1}^N \eta\_i \varepsilon\_i \right\}. \tag{7}$$

The first term of Equation (7) on the right side is completely dependent on the other differences and represents the predictable part of Δ**y***j*; the bracketed expression on the right side is unpredictable and independent of the other measurements. When the first term in the bracketed expression does not exceed the second random measurement error part, Equation (7) provides a way for predicting Δ**y***<sup>j</sup>* from other measurements, in other words, the Δ**y***<sup>j</sup>* is redundant information and can be obtained by others; this is based on:

$$\left| \sum\_{i=1}^{N} \eta\_i \mathbf{J}(\mathbf{x} \star\_i \mathbf{a}\_i) \cdot \Delta \mathbf{x} \right| \le \left| \sum\_{i=1}^{N} \eta\_i \varepsilon\_i \right|. \tag{8}$$

It follows that the smaller the left term of Equation (8), the better the prediction. The question of mutual independence can thereby be examined by considering the minimal value of <sup>η</sup><sup>T</sup> ·**J**(**x**∗, **<sup>a</sup>**) 2 , which can be achieved by a suitable selection of η*i*. Here, η is a *N* dimensional vector consisting of elements η*<sup>i</sup>* (*i* = 1, 2, *..., N*).

According to the Schwartz inequality, the minimum value of <sup>η</sup><sup>T</sup> ·**J**(**x**∗, **<sup>a</sup>**) <sup>2</sup> is equal to the smallest <sup>ζ</sup>min of the eigenvalues of the matrix **<sup>J</sup>**(**x**∗, **<sup>a</sup>**)·**J**T(**x**∗, **<sup>a</sup>**) when the <sup>η</sup>*<sup>i</sup>* values are equal to the *<sup>N</sup>* components of the corresponding eigen-vector. If the condition by the Equation (8) holds, the difference Δ**y***<sup>i</sup>* that corresponds to the smallest η*<sup>i</sup>* of the eigen-vector for ζmin is the most linearly dependent. Consequently, the value of Δ**y***<sup>i</sup>* can be calculated by the first term on the right side of Equation (7) from the other measurement differences with no loss of information. Therefore, we can identify the redundant measurement configuration by:

$$\mathbf{a}\_{i}^{\text{redundant}} = \operatorname\*{argmin}\_{\mathbf{a} \in \mathbf{P}} [\eta(\mathbf{a}\_{i})],\tag{9}$$

here, P is the associated configuration domain.

The optimization procedure can be summarized as follows:

Firstly, give a slight deviation from the nominal dimensions of the investigated sample by calculating the signature difference Δ**y***<sup>i</sup>* = [Δ*yi*1, Δ*yi*2, ... , Δ*yil*] <sup>T</sup> (*i* = 1, 2, ... , *N*) under *N* measurement configurations according to Equation (5);

Secondly, for each element Δ*yij* (*j* = 1, 2, ... , *l*) in Δ**y***i*, carry out the eigen-analysis according to Equations (6)–(9) to identify and eliminate the measurement configurations with redundant information one-by-one; since at least *m* measured data points are required to determine the *m* unknown parameters for any mathematical inverse problem in theory [21–24], the achieved set *Sj* will contain *m* non-redundant measurement configurations;

Thirdly, considering that different *Sj* may contain identical measurement configurations, the union of *Sj* for all *j* = 1, 2, ... , *l* is the final set of non-redundant measurement configurations, namely,

$$\mathcal{S}\_{opt} = \bigcup\_{j=1}^{l} \mathcal{S}\_{j}. \tag{10}$$

After that, the inverse problem described in Equation (2) should be solved under the above achieved measurement configurations, namely, **a** ∈ *Sopt*.

It is noted that the above refinement of measurement configurations involves the calculation of the Jacobian matrix **J**(**x**\*, **a***i*), which is typically only valid at the vicinity of **x**. To ensure that the refinement results of measurement configurations are valid in the changes of profile parameters, it is necessary to repeat the above refinement process of measurement configurations at *k* different values of **x**. In this case, the final set of non-redundant measurement configurations will be the union of the achieved set at each value of **x**. Here, the value of *k* depends on the types of measured signatures as well as the complexity of the nanostructure under measurement. As a rule of thumb, *k* should take a relatively large value for a complicated nanostructure with a simple type of measured signature, such as reflectance.

#### **3. Experiments**

#### *3.1. Experimental Setup*

We took spectroscopic scatterometry as an example to demonstrate the capability of the proposed method in identifying and eliminating the measurement configurations with redundant information. The experiments were carried out by a commercial Mueller matrix ellipsometer (ME-L ellipsometer, Wuhan Eoptics Technology Co., China). As schematically shown in Figure 1, the system layout of the dual rotating-compensator ellipsometer in order of light propagation is PCr1(ω1)SCr2(ω2)A, where P and A stand for the fixed polarizer and analyzer, respectively, Cr1 and Cr2 refer to the 1st and 2nd frequency-coupled rotating compensators, respectively, and S stands for the sample [31,32]. The 1st and 2nd compensators rotate synchronously at ω<sup>1</sup> = 5ω and ω<sup>2</sup> = 3ω, respectively, where ω is the fundamental mechanical frequency. With the light source used in this ellipsometer, the spectral range is from 200 to 1000 nm, covering the spectral range of 300–900 nm used in this work. The beam diameter can be changed from the nominal values of ~3 mm to a value less than 200 μm with the focusing lens. The two arms of the ellipsometer and the sample stage can be rotated to change the incidence and azimuthal angles in experiments. Except for the reflectance and ellipsometric angles of the sample under measurement, the 16 Mueller matrix elements also can be obtained with the dual-rotating compensator setting. An in-house developed MATLAB® (version R2017a, The MathWorks, Inc., Natick, MA, USA) program for analyzing the measured signatures ran on a workstation equipped with double 2.0 GHz Intel Xeon CPUs. The forward optical model in this program was developed based on rigorous coupled-wave analysis (RCWA) [33–35], and the inverse optical scattering problem was solved through the commonly used Levenberg-Marquardt (LM) algorithm [36].

**Figure 1.** Principle and instrument of the dual rotating-compensator Mueller matrix ellipsometer.

#### *3.2. Experimental Results on 2D Grating*

The first investigated sample was a 2D Si grating whose SEM cross-section image is shown in Figure 2. The optical properties of Si were taken from Reference [37]. As depicted in Figure 2, the cross section of the Si grating could be characterized by a symmetrical trapezoidal model with top width (*W*), grating height (*H*), sidewall angle (*SWA*), and period (*P*), whose dimensions obtained from Figure 2 were *W* = 350 nm, *H* = 472 nm, *SWA* = 87.63◦, and *P* = 800 nm.

**Figure 2.** Scanning electron microscopy (SEM) cross-section image of the investigated 2D Si grating sample.

In this case, the profile parameters under study were *W*, *H*, and *SWA*, while the grating period *P* was fixed at its nominal dimension, namely, *m* = 3. The measured signature contained 15 Mueller matrix elements normalized to the first element, namely, *l* = 15. As mentioned above, the measurement configuration **a** = [ϕ, θ, λ] was defined as a combination of fixed azimuthal angle ϕ, fixed incident angle θ, but a wide waveband λ for spectroscopic scatterometry. As a simple example to demonstrate our data refinement method, here, the incident and azimuthal angles were fixed at θ = 65◦ and ϕ = 30◦, the wavelength covered a fixed range of 300~900 nm with a step of 5 nm. We assumed that the dimensions of the three profile parameters under investigation had a deviation of about 1% from their nominal values. The differences of the 15 Mueller matrix elements between the actual and nominal profiles were calculated by our forward RCWA model under each wavelength point according to Equation (5). After considering all the 15 Mueller matrix elements and the wavelengths used, the enanalysis and eigen-analysis procedures described above were performed on the differences. According to Equations (6)–(9), the wavelength points with redundant information were identified and eliminated from the set of measured signatures in a repetitive manner. Since the investigated

sample shown in Figure 2 is a simple structure, and the model output Mueller matrix had 15 elements, we achieved the 3 × 15 optimized data sets by traversing the three profile parameters and the 15 Mueller matrix elements. Through merging the same wavelengths into one for the 3 × 15 optimized data sets, the remaining unions *S*opt, containing 19 points, were obtained and shown by the black circles in Figure 3, which are supposed to contain enough information needed for the profile reconstruction of the 2D Si grating.

**Figure 3.** Fitting results of the calculated and the ellipsometer-measured Mueller matrix elements with the incidence and azimuthal angles fixed at θ = 65◦ and ϕ = 30◦, respectively.

To examine the validity of the remaining information, regression analysis was conducted to reconstruct the three profile parameters from the full spectral range of 300–900 nm, the optimal spectrum containing 19 points, and the spectrum of 19 random points. The non-linear LM algorithm was applied to fit the measured 15 Mueller matrix elements with the modeled ones, which converged very quickly to a minimum when a suitable initial condition was chosen. During the fitting procedure, the root-mean-squared error (RMSE) was used to quantify the goodness of fit, which summed over all the measurement wavelengths the differences between the measured data and model-generated data. The lower the RMSE value, the better the fit or agreement between the measured and model-generated data. Figure 3 presents the fitting result of the model-calculated best-fit Mueller matrix elements (black circles) and the measured spectrum including (red lines) 19 optimal points.

Table 1 shows the extracted profile parameter values for the fits and the SEM-measured values. The uncertainties associated with the SEM-measured values were estimated by manually measuring the SEM micrographs. The uncertainties appended to the scatterometry-determined values were estimated from 30 repeated measurements. All uncertainties correspond to a 95% confidence level. As shown in Table 1, the parameters *W*, *H*, and *SWA* of the 2D grating, extracted by fitting the optimized spectrum, approximately agree with the ones obtained by using the full spectrum, while the former method is over three times faster than the later one; the time spent between the optimal and random measuring modes is slightly different, and the reason may be that the former iterations are a little

smaller than the latter during the profile reconstruction procedure as that latter has more sufficient information. Figure 3 shows that the model-calculated 15 Mueller matrix elements match well with the measured spectrum. The experiments performed on the 2D Si grating provide preliminary evidence that the advantages of the proposed method are in speeding up the reconstruction procedure and effectively improving the precision of the extracted parameters.

**Table 1.** Comparison of fitting results of the 2D grating extracted from different spectrum and SEM measurements.


#### *3.3. Experimental Results on 3D Grating*

The second sample was a photoresistive grating consisting of a cylinder array deposited on silicon with a bottom anti-reflective coating (BARC), which is a 3D grating whose SEM cross-section image is shown in Figure 4. The 3D grating of the cylinder array was chosen for this study due to the high computation effort of its forward optical model. As depicted in Figure 4, the cross section of the 3D grating is characterized by a model with cylinder diameter (*D*) and height (*H*1), BARC height (*H*2), and periods of arrays in two directions (*P*<sup>x</sup> and *P*y). The dimensions obtained from Figure 4 are *D* = 226.7 nm, *H*<sup>1</sup> = 355.1 nm, *H*<sup>2</sup> = 104.3 nm, and *P*<sup>x</sup> = *P*<sup>y</sup> = 835 nm. In our experiments, optical properties of the photoresist were modelled by a two-term Forouhi–Bloomer model [38], whose eight undetermined parameters were determined by measuring a photoresist film deposited on the silicon substrate using the above ellipsometer and were *A*<sup>1</sup> = 0.006029, *A*<sup>2</sup> = 0.020598, *B*<sup>1</sup> = 14.195303 eV, *<sup>B</sup>*<sup>2</sup> <sup>=</sup> 14.196431 eV, *<sup>C</sup>*<sup>1</sup> <sup>=</sup> 50.523937 eV2, *<sup>C</sup>*<sup>2</sup> = = 50.537878 eV2, *<sup>n</sup>* (∞) <sup>=</sup> 1.436087, and *Eg* <sup>=</sup> 4.774050 eV. Similarly, optical properties of the BARC were modelled using the Tauc–Lorentz model [39], whose five undetermined model parameters were obtained as ε<sup>1</sup> (∞) = 1.42680, *Eg* = 3.45971 eV, *A* = 21.14955, *E*<sup>0</sup> = 9.94921, *C* = 0.98767. For more details about physical meanings of the above parameters, one can consult References [38,39], which are omitted here for brevity.

**Figure 4.** SEM cross-section image of the investigated 3D grating sample.

Here, the investigated profile parameters of the 3D grating included the *D*, the *H*1, and the *H*2, namely, *m* = 3, while the two grating periods were fixed at their nominal dimensions. In contrast to the experiment settings above, the incidence and azimuthal angle were fixed at θ = 65◦ and ϕ = 0◦, and the terms of *N*, *C*, and *S* parameters versus wavelength were used as measured signatures, which

were derived from ellipsometric Psi (Ψ) and Delta (Δ) and taken as *N* = cos (2Ψ), *C* = sin (2Ψ) cos (Δ), and *S* = sin (2Ψ) sin (Δ). Considering that the investigated 3D sample deviated largely from its nominal values, its forward model was highly complex and nonlinear, and the model output was *N*, *C*, and *S*, namely, *l* = 3, we assumed that the dimensions of the three profile parameters changed *k* = 6 times to ensure that the optimization procedure was stable, and the deviations were about ±1%, ±5%, or ±10% from their nominal dimensions, respectively. According to our proposed method, 3 × 3 × 6 optimized data sets remained, then, the union *S*opt of the remaining data sets, consisting of the 36 wavelength points shown in Figure 5, were achieved for the inverse problem solution.

**Figure 5.** Fitting results of the calculated and the ellipsometer-measured *N*, *C*, and *S* parameters at the incidence angle θ = 65◦ and azimuthal angle ϕ = 0◦.

In the same way, the regression analysis was performed to extract the three unknown profile parameters, and the non-linear LM algorithm was also applied to fit the measured signatures with the modeled ones. The extracted profile parameter values for the fits and the SEM-measured values are shown in Table 2. The uncertainties appended to the SEM-measured and scatterometry-measured values in Table 2 were estimated in a similar manner to those in Table 1. Figure 5 presents the fitting result of the model-calculated best-fit *N*, *C*, and *S* parameters and the measured spectrum (black dots) containing whole range (blue lines) and only 36 optimal points (red circles), respectively.

**Table 2.** Comparison of fitting results of the 3D photoresist grating extracted from different spectrum and SEM measurements.


From Table 2, the three average values of the 3D grating dimensions extracted using 36 optimal points are very close to the values obtained using the full range spectrum, but the extracted parameter *D* has about a 25 nm deviation from the SEM-measured value, which is probably caused by the obvious non-uniformity of the grating surface depicted in Figure 4, and the SEM-measured values include estimation errors in manually measuring the SEM micrographs. Table 2 also shows that the proposed method results in the lowest value of RMSE, the highest measurement precision, and the least time cost. It is observed from Figure 5 that the best-fit *N*, *C*, and *S* parameters using the optimized spectrum are in better agreement with the measured signatures compared to those of the full range. Consequently, we may conclude that the proposed dependence-analysis-based data-refinement method can be applied to determine several optimal measurement points without any loss in the accuracies for profile reconstruction in optical scatterometry.

#### **4. Conclusions**

In summary, we proposed a method to identify and eliminate the measurement configurations with redundant information based on dependence analysis in optical scatterometry. By assuming that the dimensions of profile parameters under investigation have some deviation from their nominal values, the differences between the actual and nominal dimensions were calculated by the forward optical model under each measurement point. A formulation was derived to identify the measurement configurations with redundant information through performance of a dependence analysis followed by an eigen-analysis. By eliminating redundant information from the measured data in a repetitive manner, a few optimal points remained and were used for the reconstruction process.

Experiments performed on a 2D Si grating and a 3D photoresist grating revealed that the reconstructed grating profiles from the optimally selected subset of measurement wavelengths according to the proposed method have a higher accuracy than those from the randomly selected counterparts with a higher efficiency than the analysis making use of the whole spectrum. This suggests that the proposed dependence-analysis-based data-refinement method can be a powerful tool to enhance the reconstruction speed of nanostructure metrology using scatterometry without sacrificing the accuracies, especially for the nanostructures whose forward model is very complex and time-consuming.

**Author Contributions:** Conceptualization, X.C. and Z.D.; methodology and formal analysis, Z.D. and Y.S.; writing—original draft preparation, Z.D.; writing—review and editing, X.C., Z.D., X.W., Y.S., H.J., and S.L.; supervision, Z.D., Y.S., and X.C.; project administration, Y.S. and X.C.; funding acquisition, Z.D., Y.S., X.C., and S.L.

**Funding:** This research was funded by the National Natural Science Foundation of China (Grant Nos. 51775217, 51727809, and 51525502), the Natural Science Foundation of Hubei Province of China (Grant Nos. 2018CFB290 and 2018CFB559), the China Postdoctoral Science Foundation (Grant Nos. 2016M602269 and 2019M652633), the National key research and development program of China (Grant No. 2017YFF0204705), and the National Science and Technology Major Project of China (Grant No. 2017ZX02101006-004).

**Acknowledgments:** The authors would like to thank Zhimou Xu from the School of Optical and Electronic Information of Huazhong University of Science and Technology, and Shanghai Micro Electronics Equipment Co., Ltd. (Shanghai, China) for preparing the samples.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


© 2019 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).

*Article*
