Nondestructive Photoelastic and Machine Learning Characterization of Surface Cracks and Prediction of Weibull Parameters for Photovoltaic Silicon Wafers

Crystalline silicon is the leading substrate used in solar cells today and accounts for over 85% of the global photovoltaic (PV) market [1,2]. The increasingly competitive nature of the industry has driven solar cell manufacturers to seek new methods of reducing the overall cost of PV cells. This is accomplished in part by reducing the thickness of silicon substrates in order to obtain more wafers per ingot [3–6]. Given the brittle nature of silicon, thinner wafers are more prone to fracture, which is a common occurrence during wafer handling steps and throughout the solar cell production process. This is particularly costly if a wafer fractures during the screen printing stages, which are later in the production process, requiring an operator to stop the production line to clean off aluminum paste and wafer shards [7]. The presence of residual stress and defects within the wafers increases the likelihood of fracture and can reduce the production yield by up to 25 percent [7]. Between 5 and 15 percent of wafers fail during production due to defects such as microcracks [2,8–10].

The wafers used in this study were cut from ingots that were formed via the Czochralski process, which is the most commonly used method for obtaining monocrystalline silicon. The Czochralski process begins by melting hyperpure silicon in a quartz crucible. A seed crystal with the desired crystal orientation is lowered into the melt, rotated, and gradually pulled away from the melt. The pull rate determines the diameter of the silicon ingot and is varied so that the ingot begins with a narrow, dislocation-free neck and increases to a diameter slightly larger than the desired wafer diameter. Since impurity atoms are more likely to remain in the melt than to become part of the ingot, the density of impurity atoms increases from the top to the bottom of the ingot as the concentration of impurity atoms increases in the melt. Once the ingot has cooled, the outer diameter is reduced to the desired wafer size and the conical top and bottom sections of the ingot are sliced off. The sides of the ingot are cut off to produce a square-like cross section that allows for a higher packing density of solar cells. The ingot is then sliced into wafers by a wire saw. Wire saw slicing is the preferred method of wafering slicing silicon ingots because it simultaneously produces several hundred wafers and the use of thin wires produces less kerf than an internal diameter saw.

Microcracks primarily arise as a result of the wire saw slicing process. The saw damage layer induced by loose abrasive wire sawing (LAWS) is localized to approximately the top and bottom 20–40 μm of each wafer [4,11,12]. The depth of the microcracks depends on the size and shape of the abrasive, SiC particles [11]. Post-processing steps such as etching are used to mitigate the effects of surface microcracks which in turn increases both the solar cell efficiency and mechanical strength of the wafer [4,13].

Barredo et al. studied the relationship between reduction in saw damage layer by etching and the mechanical strength of silicon wafers. Wafer strength was observed to increase from 60 MPa in the as-cut state to 260 MPa after 30 μm were removed from each face [4]. Möller et al. noted that typical cracks have a depth of up to 10–15 μm and observed a gradual increase in wafer strength after etching away the top 5, 10, and 20 μm of the saw damage layer [11]. Pogue estimates that the saw damage layer is 24–40 (16–25) μm thick for wafers sliced by LAWS (diamond wire saw). After etching, the wafer strength remains far below the ultimate tensile strength of silicon which is due to the continued presence of microcracks. Etching reduces the depth of surface cracks and blunts crack tips but it does not remove the entire saw damage layer [9,14,15]. Typically, the top 10–20 μm of the saw damage layer is removed during the etching process [9,10].

To address the issue of yield loss due to microcracks, early in the cell fabrication process, wafers are inspected for critical cracks or those that are likely to result in wafer failure. Wafers that are found to contain a critical flaw are recycled while other wafers are fabricated into solar cells. There are many optical and acoustic approaches to detecting critical microcracks, several of which have been outlined in detail in Refs. [6,7,13]. The challenge arises in accurately identifying critical microcracks. Fractographic analysis has shown that critical microcracks are favorably oriented semielliptic surface cracks that are approximately 10 μm in length [16]. The crack plane is considered to be favorably oriented when it is approximately orthogonal to both the externally applied load and the wafer surface.

Failure to correctly identify wafers that contain critical microcracks results in processing flawed wafers while incorrectly classifying wafers that do not contain critical microcracks results in rejecting wafers that would survive the fabrication process. These are respectively known as type I and type II errors [17,18]. Type I errors occur when cracks are not properly resolved, resulting in wafers being processed that have higher than acceptable probabilities of fracturing. Type II errors arise when wafers are incorrectly characterized as having critical flaws, leading to unnecessary material loss. We aim to minimize type I errors by utilizing a high-resolution photoelastic imaging system to detect defects on the same size scale as critical microcracks and to minimize type II errors by using a multistep validation process which each defect must pass before it is considered to be a microcrack.

Each microcrack experiences a bulk applied stress that consists of thermal residual stresses from the ingot cooling process and those that arise from the wire sawing process. This residual stress gives rise to a locally elevated crack tip stress field. In order to characterize each microcrack, we must decouple the bulk and local stress fields observed in photoelastic measurements. Similarly, this decoupling of local and bulk residual stress fields has been accomplished for subsurface inclusions in multicrystalline silicon wafers [19].

In this work, we characterize 18 Czochralski grown, etched, silicon (100) wafers according to fracture strength under a virtual uniaxial tensile load. The virtual load refers to a uniform theoretical tensile load that is computationally superimposed with the wafer’s bulk residual stress field and acts upon the distribution of observed microcracks. We use the distribution of fracture strengths to calculate the Weibull parameters of the wafer set and then demonstrate the Weibull parameters’ dependence on external loading conditions. To do so, we use a machine-learning algorithm to locate potential surface microcracks, decouple the local and bulk stress fields near each microcrack, and characterize each microcrack by comparing the local stress field to the theoretical crack tip stress field. The distribution of microcracks in each wafer is subjected to a virtual uniaxial tensile load and each wafer is characterized according to the stress at which the first crack is expected to propagate. Since the direction along which the virtual uniaxial load is applied is variable, it can be used to compare the expected strength of a wafer when loaded along and perpendicular to the wire speed direction. The method could be used to improve upon current manufacturing practices for crystalline photovoltaic wafers.

The IR-GFP in this work uses a narrow band circularly polarized light source centered at approximately 1150 nm. When circularly polarized light experiences a linear retardation as described in Eq. (1), the light gains an ellipticity proportional to the magnitude of shear stress in the wafer. The intensity of light is monitored by a charge-coupled detector which resides behind a rotating analyzer. Thus, the image is brightest when the analyzer aligns with the major axis of the polarization ellipse and dimmest when the analyzer aligns with the minor axis of the polarization ellipse.

The degree of ellipticity and orientation of the ellipse are determined via a lock-in algorithm which monitors the intensity of each pixel at 16 distinct analyzer angles. Since the degree of ellipticity is proportional to the magnitude of the shear strain, this can be used to calculate the shear strain that lies in the x–y plane, the shear strain that lies in the $±45−deg$ plane, and the maximum shear strain. A more in-depth summary of this instrument is given in Ref. [20]. Because the measured ellipse is the result of the net retardation through the thickness of the wafer, when any variation in shear stress through the wafer is present, the variation in the z-direction must be accounted for in Eq. (1). This is the case for many microscopic flaws such as cracks or voids where the shear stress is greater near the flaw. As such, when observing microscopic flaws, we refer to the net retardation caused by x–y shear stress (δ₀) and $±45−deg$ shear stress (δ₄₅) instead of the measured shear stress.

The spatial resolution of the IR-GFP can be increased to approximately 5.5 μm with the addition of a 5× objective. This allows us to resolve bow tie features in the δ₀ and δ₄₅ images which result from locally elevated stresses. Examples are shown in Fig. 1 for a defect in a typical PV-grade crystalline silicon wafer with (001) orientation.

Rotating $σ¯¯$ by α, the angle between the crack plane and [100] direction on the wafer, aligns the crack tip stress field coordinate system with the x–y coordinate system of the GFP. We will refer to the rotated stress field as σ. Then, substituting Eqs. (7)–(9) into Eqs. (2)–(5) provides the form for the δ₀, δ₄₅, and δ_M fields shown in Fig. 2.

The δ fields of several thousand surface cracks have been measured over a total of 18 photovoltaic silicon wafers. Each δ field is analyzed in order to characterize each crack by its stress intensity factors, the geometric factor which relates applied stress to the stress intensity factors, and the ratio of bulk principal stresses acting on each crack (η). δ^(∞) and η will be used to determine the bulk biaxial stress field acting on each crack. Finally, in accordance with the weakest link theory, the strength of each wafer is calculated based on the most detrimental surface crack.

The probability of survival of the ith volume element is represented by $ps(i)$ and the probability of survival of the system is represented by P_S. This statistical approach to predicting the probability of survival is known as the weakest link theory. Each volume element can be thought of as a link in a chain, and when the weakest link breaks, the system will fail. This approach applies well to silicon wafers since failure can be assumed to occur when the first crack begins to propagate under tension [24].

We base the survival of each volume element, an 800-by-600 μm rectangular region, on the bow tie feature with the greatest δ^(loc) value. Volume elements that contain a bow tie feature caused by a microcrack are characterized by the uniaxial tensile stress required to increase the mode I stress intensity factor of the crack to equal the fracture toughness of silicon. Thus, when the mode I stress intensity factor equals the fracture toughness of silicon, the probability of survival of the volume element is zero. Like the stress optic coefficient, the fracture toughness has also been shown to be anisotropic, varying with the orientation of the crack plane [26].

This section discusses the methods used to collect δ fields, process images, identify and analyze microcracks, and calculate wafer strengths. The PV wafers used in this study, provided by a leading solar cell manufacturer, are Czochralski grown, have been sliced via a loose abrasive wire saw and have had 13 μm removed from each surface during the etching process. Etched wafers were selected for this study because they have a spatially sparse population of microcracks and a slightly higher signal-to-noise ratio when compared to as-sawn wafers. The sparseness of microcracks in etched wafers facilitated the process of identifying crack tip stress fields since it is less likely that two different locally elevated stress fields will overlap spatially. The increased signal-to-noise ratio may be attributed to the slightly smoother surface of etched wafers which produces less scattering of the light source.

The IR-GFP used in this work is capable of recording one 640-by-480 pixel array for each δ field measurement. A full field image of a wafer has a spatial resolution of approximately 325 μm. As can be seen in Fig. 1, a resolution of at least 80 μm is required to resolve spikes in the δ field associated with the presence of a surface crack. The addition of the 5× objective provides 5.5 μm spatial resolution which permits not only the resolution of locally elevated δ values but also the fluctuations in intensity around each bow tie pattern. The increased spatial resolution comes at the cost of a smaller field of view, approximately 3.5 by 2.6 millimeters. Thus, over 3000 images are required to obtain a full δ field of a 156 by 156 mm wafer at 5 − μm resolution. The task of measuring each wafer’s δ field has been automated by using an x–y stage, retrofitted with stepper motors, to raster scan each wafer across the optic path.

There are hundreds of resolvable bow tie defect features on each etched wafer and thousands on each as-cut wafer. There are many more which can be resolved at higher magnification, but we are only concerned with cracks that may initiate wafer fracture, so it is sufficient to only observe the population of cracks exhibiting the largest δ^(loc) field.

To identify most of these bowties, a support vector classifier (SVC) is used to classify the profile of δ₀ and δ₄₅ fields, near locally elevated retardation levels, as a bow tie or non-bow tie. The profile of each δ field is measured by sweeping a 20 − μm linescan in a circle around the point of maximum retardation. The average value of each linescan is recorded and then scaled to range from 0 to 1 as shown for $δM(loc)$ in Fig. 4.

The SVC is trained by a set of δ profiles consisting of 500 bowties and 650 non-bowties that are manually identified from wafers that are not included in the 18 wafers evaluated here. Each δ profile consists of a list of 100 retardation values where each value in the list is an attribute of the δ profile. When training the SVC, each δ profile is considered to be a 100-dimensional vector and a kernel is used to map each delta profile point so that the clusters of bow tie and non-bow tie points are distinctly separate from each other and a hyperplane can be placed between them. All future δ profiles can then be classified as bow tie or non-bow tie depending on which side of the hyperplane they reside.

A set of 224 manually identified bow tie and non-bow tie features, separate from the training set, are used to test the accuracy of the SVC and determine the optimal combination of kernel method, regularization parameter, and gamma parameter. The regularization parameter dictates how much error is allowed when placing the hyperplane between the two classes of data points. A large regularization parameter allows for very little error and focuses on classifying all of the data points correctly whereas a small regularization parameter allows for a larger margin of error and produces a less complex hyperplane that may misclassify a larger portion of the training data points. The gamma parameter decides which data points are considered when forming the hyperplane. A large gamma value only takes into consideration points nearest to the plane while a small gamma value considers points both close to and far from the plane.

To optimize the SVC, a grid search was conducted using a variety of kernel methods, regularization parameters, and gamma parameters. The optimal combination achieved 93% accuracy, used a radial basis function kernel, a gamma value of 0.0001, and a regularization parameter of 100, where larger regularization parameter values yielded equally accurate results.

Using the SVC to classify a single bow tie is relatively computationally efficient, but classifying every pixel in each image would be prohibitively costly. To reduce computation time, instead of checking every pixel in each image, we divide each δ field image into 16 sub-regions and only the location of maximum retardation in each sub-region is checked for a bow tie. For every bow tie that is identified by the SVC, we record the bowties’ $δM(loc)$ field where the δ_M field is calculated by inserting the measured δ₀ and δ₄₅ fields into Eqs. (2)–(5). We also record the bulk first principal stress direction $(β(∞)(1))$ and the magnitude of $δM(∞)$⁠. $β(∞)(1),δM(∞)$⁠, and the profile of $δM(loc)$ are used to characterize each microcrack by its crack plane orientation, geometric factor, and the bulk biaxial stress acting on the crack as described in the following sections.

Next, we characterize the profile of $δM(loc)$ which will be compared to the measured profile of $δM(loc)$⁠. An array of $η=σ(∞)(1)/σ(∞)(2)$ and α input values is used to check which combination produces the best match between the theoretical and measured $δM(loc)$ profile. The resolution of the input array used when calculating the theoretical $δM(loc)$ profile is 2 deg with respect to α and 0.02 with respect to η.

The local retardation field is derived for a surface crack as a function of α, $β(∞)(1)$⁠, $β(loc)(1)$⁠, wafer thickness t, angle between the location on the crack front and wafer surface ϕ, in-plane radius from crack tip ρ, and the ratio of bulk residual principal stresses η.

Values for α and η when inserted into Eq. (26) provide profiles of $δM−Norm(loc)$ for a variety of crack orientations and bulk residual principal stress ratios. Each profile is calculated for θ ranging from 0 to 2π and then scaled to range from 0 to 1. As shown in Fig. 4, computed profiles are compared to measured profiles for individual bow tie patterns.

Each microcrack is characterized by the combination of α and η that corresponds to the normalized theoretical $δM(loc)$ profile that best matches the normalized measured $δM(loc)$ profile. The best matching profile is determined by the method of least squares. Only bowties that fit the profile of a crack tip retardation field whose plane is normal to the wafer surface are accepted as crack tip stress fields. This helps to avoid misdiagnosing other defects that produce locally elevated stress fields as microcracks. Furthermore, the population of bow tie profiles accepted as microcracks by this method will consist of cracks whose plane is approximately orthogonal to the wafer surface.

Once α and η have been determined by the aforementioned process, the ratio of stress intensity factors (K_R) can be calculated from Eq. (24). Then, Eq. (19), when set equal to the measured $δM(loc)$ at a known location (ρ, θ), can be used to solve for the mode II stress intensity factor. Finally, utilizing Eq. (19), the mode I stress intensity factor is obtained.

Measurements are in stress· thickness and can be converted to stress by dividing by the thickness over which the retardation accumulates. The horizontal line features are a result of the raster scanning method used to collect 5× images. Elevated or low values along the immediate edge of the wafer are a result of the light source reflecting off of the internal edge of the wafer at or near Brewster’s angle, thus altering the polarization. While the shape of the measured difference in principal stress field is consistent from wafer to wafer, it is not consistent with the theoretical thermal residual stress field that might be expected to arise in a wafer cut from a 200-mm diameter Czochralski grown ingot [29,30]. It is possible that the asymmetry observed in the difference in principal stress fields is a result of uneven cooling while the silicon ingot solidified resulting in a thermal gradient across the ingot. Notice also that the measured magnitude of difference in principle stresses varies by only about 20% from minimum to maximum across each wafer indicating that the background stress is actually relatively uniform. We speculate that the magnitude of the difference in the principal stress field is potentially the result of two factors. First, while the illumination source is a narrow-band light source centered at 1150 nm, it is not monochromatic and thus cannot be perfectly circularly polarized. Any deviation from circular polarization will appear as retardation in the measured data. Prior to imaging, a stress-free electronic grade wafer is measured to check how circularly polarized the light source is. A stack of waveplates is then adjusted to form a quarter waveplate to ensure the light source is as close to circularly polarized as possible. The second factor is the presence of the saw damage layers. The theoretical stress fields inferred in Ref. [29] only account for thermal residual stresses and dislocation stress fields. It is possible that the elevated difference in principal stresses observed in Fig. 5 is the result of a combination of thermal residual stresses and damage induced by the wire sawing process. Popovich et al. [15] used micro-Raman spectroscopy to measure residual stress in the saw damage layer of multicrystalline silicon wafers and under the assumption of a biaxial stress state found the stresses to be approximately 500 MPa.

Having determined the mode I stress intensity factor of each crack and the normal component of the bulk residual stress acting on each crack, we make the assumption that the relationship between K_I and $σ⊥(∞)$ remains linear until the crack geometry changes upon crack propagation. We refer to this relationship as the crack geometric factor, χ.

The form of Eq. (34) was selected to satisfy the fracture toughness values presented by Li et al. and to mimic the anisotropic behavior observed in the (100) silicon stress optic coefficient and elastic modulus [26]. Using Eq. (33), $σC(appl)$ is calculated for all of the analyzed microcracks and the lowest value is used to characterize the strength of each wafer. Finally, according to the methods described in Sec. 1.3, the Weibull parameters for the batch of wafers are calculated from the set of wafer strengths.

Figure 6 shows the distribution of wafer strengths, as calculated by the methods described earlier, on a Weibull plot. The Weibull parameters shown in Fig. 6 are dependent on ψ, the direction of the applied uniaxial tensile load.

Figure 7 shows how the characteristic strength and Weibull shape parameter of the wafer batch vary with the direction of applied uniaxial tension for three assumed saw damage layer thicknesses.

The characteristic strength of the wafer batch calculated here is greater when the wafer is loaded perpendicular to the wire speed direction than when loaded parallel to the wire speed direction. The direction of wire and ingot motion are, respectively, along $ψ=90deg$ and $ψ=0deg$⁠. As shown in Table 1, this variation in wafer fracture strengths has been observed in both slurry and diamond wire sawn multicrystalline silicon wafers [32].

When comparing the characteristic strength determined by the virtual uniaxial tensile load to the characteristic strength determined by four-line bend tests one must consider the difference in stressed area and types of defects present. Four-line bend tests only apply tension to the lower half of the wafer in the region between the inner loading pins, whereas the virtual uniaxial tensile load applies tension to the entire wafer. Since the probability of encountering a critical microcrack increases with the size of the stressed region, the characteristic strength will be lower for wafer groups fractured by uniaxial tensile stress than for wafer groups fractured by four-line bend tests. However, the characteristic strength of wafers tested by four-line bend tests is also influenced by edge defects [33]. Polarization of the light source induced when reflecting off the internal edge of the wafer at or near Brewster’s angle prevents the method presented here from accounting for microcracks near the wafer edge. The omission of edge defects from the analysis may result in an overestimation of the characteristic strength.

The characteristic strength of wafers loaded perpendicular to the wire speed direction has been observed to be 88 percent greater for diamond sawn wafers and 11 percent greater for slurry sawn wafers than when loaded parallel to the wire speed direction [32]. The set of wafers studied here has a characteristic strength which is 30 percent greater when loaded perpendicular to the wire speed direction. In the same study [32], the authors also observed increased Weibull shape modulus when loading the samples parallel to the saw damage, suggesting less scatter in the set of wafer strengths. This is also observed in our calculated shape parameter values shown in Fig. 7. When loaded perpendicular to the wire speed direction, the shape modulus is 1.6 and when loaded parallel to the wire speed direction the shape modulus is 1.9. However, varying results have been observed with respect to the characteristic strength of wafers loaded parallel and perpendicular to the wire speed direction. Bidiville et al. found diamond wire sawn wafers to have approximately double the characteristic strength when loaded perpendicular to the wire speed direction compared to along the wire speed direction [34], while Bohne et al. observed slurry sawn monocrystalline silicon wafers to have approximately the same characteristic fracture strength when loaded parallel and perpendicular to the wire speed direction [35].

The variance in characteristic strength and Weibull shape parameter appears to originate from the distribution of observed crack orientations and preferential first bulk residual principal stress direction shown in Fig. 8.

All of the microcracks represented in Fig. 8(a) are believed to be approximately orthogonal to the wafer surface. This is because Eqs. (18)–(27) for theoretical $δM(loc)$ assume the crack is orthogonal to the surface. Thus, when the theoretical and measured $δM(loc)$ profiles are compared to determine α and η, cracks that are orthogonal to the surface are the best match to the theoretical profile. Furthermore, using the detection methods outlined in Sec. 2.2.1, cracks whose planes are oblique with respect to the wafer surface are less likely to be detected when identifying bowties since they result in lower in-plane shear stress. At first glance, the clear bimodal distribution suggests that cracks oriented orthogonal to the surface typically lie along directions 35 deg and $135deg$ from the [100] orientation. However, the crack detection methods used here influence the observed bimodal distribution of crack plane orientations. This is because in each observed region, we analyze the crack exhibiting the greatest local retardation, which, as seen in Eq. (18), is proportional to the magnitude of the stress intensity factors of the crack. The stress intensity factors depend on the crack geometry and the orientation of the crack with respect to the principal bulk residual stresses acting on the crack. Cracks oriented along the first or second bulk principal stress direction will have a mode II stress intensity factor that is equal to zero. Since the residual stresses in the saw damage layer are compressive in nature [12], the mode II stress intensity factor has a strong effect on the magnitude of local retardation. As shown in Fig. 8(b), the most common first bulk principal stress direction is along the wire speed direction. Thus, with this method, the likelihood of characterizing a crack whose plane is oriented skew to the wire speed direction is greater than that of a crack whose plane is parallel or perpendicular to the wire speed direction. This suggests that this method is more sensitive to detecting cracks that are oriented skew to the direction of wire motion, not that cracks of these orientations are more prevalent. The characteristics of saw damage-induced cracks depend on several variables such as the shape and wear of the abrasive grit [36,37].

It should be noted that Weibull parameters depend on the loading conditions used to fracture each wafer. The greater the area loaded in tension, the greater the probability of encountering a critical flaw; four-line bending tests, for example, yield lower characteristic strengths than ring-on-ring tests, though this is also related to the introduction of additional failure modes along the edge or corners of the wafers [33].

In this work, the uniaxial tensile load used to characterize each wafer’s strength was in part selected because of its uniformity across the entire wafer and because we do not know which surface each crack resides on. Under a uniaxial tensile load, a crack on the top or bottom surface of the wafer will experience the same applied stress. Furthermore, unlike four-line bending where only a fraction of one surface is in tension and the other in compression, uniaxial tension stresses the entire top and bottom surface, thus increasing the number of active surface cracks that may initiate wafer failure.

Assuming a uniform distribution of flaws across the wafer and using the distribution of geometric factors and assuming that cracks are randomly oriented, it can be shown that the strength of the wafer group via four-line bend tests could be up to 80% greater than the calculated characteristic strength under uniaxial tension. However, the distribution of microcracks is typically higher near the edges than near the centers of the wafers [4].

In this work, etched Czochralski grown silicon PV wafers are inspected with a high-resolution IR-GFP in order to observe microscopic surface crack stress fields. A machine-learning algorithm was used to identify surface crack stress fields which were then compared to the ideal crack tip stress field in order to characterize each crack according to the crack plane orientation and ratio of bulk principal stresses acting on the crack. The direction of first principal stress, ratio of bulk principal stresses, and difference in bulk principal stresses were used to determine the biaxial stress state acting on each crack. The stress intensity factors for each flaw were calculated according to the bulk stress acting on the crack, the crack orientation, and the measured local retardation field. A geometric factor was derived relating the bulk normal stress acting on each crack to the crack’s mode I stress intensity factor. The distribution of crack geometric factors and crack orientations were then used to calculate the strength of each wafer under a virtual uniaxial tensile load. Finally, the distribution of wafer strengths was used to obtain the Weibull parameters for the group of wafers.

This approach brings together two new methods to allow us to nondestructively characterize individual silicon wafers according to strength, evaluate the wafer group’s Weibull parameters, and solve for the complete residual stress field in the vicinity of microcracks. The first method involves the identification and classification of microcracks in photoelastic (PE) materials via machine learning, photoelasticity, and the theoretical crack tip stress field. While this method is not limited to silicon wafers and may also find usefulness in applications for glass characterization, the method is limited by the spatial resolution of the polariscope and by the assumption of the location of each microcrack. In this work, a 5× objective was added to the optical path to obtain sufficient resolution to study the profile of the retardation fields around crack tips. This comes at the cost of a reduced field of view and thus increased wafer imaging times. As for the crack location assumption, it is supported by a significant decrease in bow-tie defect density observed between as-sawn and etched wafers. However, this assumption may not apply to silicon wafers that have been processed differently since the distribution of microcracks will vary depending on how the silicon wafer was sawn [36,37].

The second new method presented here is the use of a virtual uniaxial tensile load to calculate the strength of each wafer and utilize the spread of strengths to determine the wafer group’s Weibull parameters. This method can also be expanded on by replacing the applied uniaxial tensile stress with a known stress distribution, such as that of a Bernoulli gripper [38]. This will allow for the calculation of a different set of Weibull parameters for each specific loading condition a wafer will experience during the fabrication process. The use of a virtual uniaxial tensile load is beneficial because it engages all of the surface flaws on both sides of the wafer equally, thus providing a better estimate of the strength of each wafer than would be obtained from a bending load that only engages a fraction of the flaws on one surface. This also circumvents the obstacle of not knowing which surface a flaw resides on.

Other limitations to the methods presented include the inability to characterize cracks that reside within several millimeters of the edge of the wafer and the fact that microcracks exhibiting low locally elevated retardation fields may be overlooked. The former is due to additional polarization caused by reflections off of the internal edge of the wafer. The latter is a result of only characterizing bowties with the greatest local retardation; thus, cracks with lower stress intensity factors and cracks residing in compressive regions whose crack plane is oriented perpendicular to a compressive stress may be overlooked.

The wafer strengths calculated using this approach are in relative agreement with published wafer strengths for etched monocrystalline wafers after accounting for the increased loaded area [11,14,32,33,39]. The primary source of error in the characteristic strength calculations is in the assumed saw damage layer thickness; this effect can be observed in Fig. 7. As such, this method is useful for ordering wafers from the same batch according to their expected strength, but one must consider the effect of saw damage thickness on calculated strength when comparing wafers that have been sawn by different processes. In future applications, more definitive measurements of the depth of each saw damage layer and how the depth varies across each wafer will increase the precision of characteristic strength calculations. One potential method of isolating the thermal residual stresses from damage-related residual stresses is to pair photoelastic methods with micro-Raman spectroscopy methods [40]. Furthermore, the Weibull shape parameters are lower than those obtained experimentally but this may also be a result of the ratio of surface flaws that are engaged in uniaxial tension versus four-line bending.

The IR-GFP used in this study has been provided by StressPhotonics Inc.

• National Science Foundation (NSF) under (Grant No. CMMI13-00466; Funder ID: 10.13039/501100008982).

There are no conflicts of interest.

The data sets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.