Analysis of photoelastic properties of monocrystalline silicon

Photoelasticity is considered a useful measurement tool for the non-destructive and contactless determination of mechanical stresses or strains in the production of silicon wafers. It describes a change in the indices of refraction of a material when the material is mechanically loaded. As silicon has a diamond lattice structure, the stress-dependent change in the refractive indices varies with the loading direction. In this work, an anisotropic stress-optic law is derived, and the corresponding stress-optical parameters are measured using a Brazilian disc test. The parameters were determined to be (π11−π12)= 14.4 · 10−7 MPa−1 and π44 = 9.4 · 10−7 MPa−1. The results of this work are compared to previous works found in the literature, and the deviations are discussed.


Introduction
The photoelastic measurement of mechanical stresses is based on the birefringence caused by mechanical stresses (or strains). Birefringence describes the ability of materials to split an incident electromagnetic wave into two refracted waves instead of one (Zinth and Zinth, 2009). These two refracted waves of light show different coefficients of refraction and different states of polarization depending on mechanical stresses. This change can be measured by a polariscope, and mechanical stresses can be deduced from the measured change (Wolf, 1961;Ramesh, 2000).
Photoelasticity was first described by David Brewster in the early 19th century (Brewster, 1815(Brewster, , 1816, and the first analytical models to describe it were published some years later (Neumann, 1841;Pockels, 1889Pockels, , 1906. In the early 20th century, the first industry applications for experimental stress measurements by means of photoelasticity can be found (Wolf, 1961). Photoelasticity was a key method to determine stresses by building models out of photoelastically active materials. This was extensively used until numerical methods, e.g. the finite element method, became more powerful and convenient for the evaluation of stresses (Ramesh, 2000). However, photoelasticity is still used as an in-line measurement method for glasses to directly measure me-chanical stresses (Vivek and Ramesh, 2015). In the same manner, it could be applied to production processes for silicon wafers, as silicon is transparent with respect to infrared light. Some research, e.g. Lederhandler (1959), Brito et al. (2005), Ganapati et al. (2010), Jagailloux et al. (2016) and Herms et al. (2019), has been carried out in this field in recent years.
A mechanical stress applied to a material susceptible to stress-induced birefringence results in a change in the difference of the two polarization states and, hence, a change in the indices of refraction n 1 and n 2 : This change can be measured using a polariscope by measuring the phase difference δ between the two refracted light waves that increases with material thickness t: For mechanically and photoelastically isotropic materials, the phase difference δ is proportional to the difference in the first and second principal stresses (σ I , σ II ) for plane stress conditions: Figure 1. Comparison of different reported models for photoelasticity in silicon from the literature expressed as the difference n between the coefficients of refraction (the orientations of differently defined coordinate systems have been aligned).
The stress-optical coefficient C generally depends on the material and wavelength. It relates the difference between the principal stresses σ I and σ II to the phase difference δ.
Silicon has an inherent mechanical and photoelastic anisotropy because of its lattice structure. Therefore, a corresponding model is required to determine the photoelastic properties of silicon. For mechanically and photoelastically anisotropic materials, Eq. (3) does not hold true. Silicon's lattice structure leads to a direction-dependent behaviour. In the literature, there are several different models that aim to describe the photoelastic properties of silicon. Some studies explicitly state the model used (Liang et al., 1992;He et al., 2004;Zheng and Danyluk, 2002), whereas others only include a brief description of their approach (Giardini, 1958;Ajmera et al., 1988;Krüger et al., 2016). However, those that explicitly state their derived model lead to different analytical expressions, even when based on the same approach, as can be seen in Fig. 1.
Therefore, in this work, photoelasticity in silicon is revisited and a new model based on the same phenomenological approach for birefringence and photoelasticity as that in Neumann (1841) and Pockels (1889Pockels ( , 1906 is derived to address the disagreement among the different models in the literature. Stress-optical parameters are determined by measurement using a polariscope and a (100)-silicon wafer in a Brazilian disc test. Results from this work are then compared to works from the literature, and the deviations are discussed.

Theoretical analysis: photoelasticity in {100} silicon
Indices of refraction due to the birefringence of an unstressed material can be described using Maxwell's equations (Hecht, 2018;Zinth and Zinth, 2009). However, this approach leads to impractically long terms. In order to simplify this, an analogy is used that consists of an ellipsoid which expresses the material properties. In the literature, this ellipsoid is commonly called an index ellipsoid or an indicatrix. It is based on an imaginary light ray falling into the centre of the indicatrix. Perpendicular to this ray, a plane is constructed. The intersection of the plane and the indicatrix results in an ellipse. The lengths of the two half-axes of the ellipse correspond to the two indices of refraction. This is shown in Fig. 2. The indicatrix is a quadratic surface that generally has six independent parameters: in which B ij is a 3 × 3 symmetric second-order tensor that is material-dependent. It is called the impermeability tensor. Vector x i represents a Cartesian coordinate system which is defined in this work so that x 1 → x, x 2 → y and x 3 → z. Using the Einstein summation convention, Eq. (4) is expanded as follows: B 11 x 2 +B 22 y 2 +B 33 z 2 +2B 12 xy+2B 13 xz+2B 23 yz = 1. (5) By rotating the indicatrix to align with the (x, y, z)coordinate system used, it can be expressed using just three independent values: Here, B I , B II and B III are principal values of B ij . By definition, the indicatrix can also be expressed as the reciprocal of the squared indices of refraction (n I , n II , n III ): or, in a more general form, as The phenomenological approach by Neumann (1841) and Pockels (1889Pockels ( , 1906 to describe photoelasticity in stressed materials links the change in impermeability B ij to the mechanical stress tensor σ ij and the mechanical strain tensor ε ij , respectively, using a fourth-order tensor. In terms of the analogy used, this means that mechanical stresses or strains deform the indicatrix by changing the impermeability matrix B ij by B ij : Here, π ij kl and p ij kl are the stress-optical and strainoptical tensors, respectively. Both can be expressed by one another considering Hooke's law for linear elasticity (Narasimhamurty, 1981). Therefore, in this work, only the stress-optical relationship is considered. The change in impermeability is expressed as the difference between impermeability matrices B ij and B o ij in a stressed and unstressed state. Using Eq. (8), they can be expressed by the refraction indices n ij and n o ij : By converting them to a common denominator, this can be rewritten as follows: The change in impermeability is assumed to be small compared with the unstressed impermeability, meaning that n ij ≈ n o ij . Thus, two simplifications can be made: With the above-mentioned simplifications, Eq. (11) can be expressed as follows: This links a change in the indices of refraction to the change in impermeability B ij and, therefore, to mechanical stresses using Eq. (9). In the following, the approximation sign is omitted, although it is still an approximation that is only valid for small changes in impermeability. Rearranging Eq. (13) yields the following: For an indicatrix whose half-axes are aligned to the coordinate system, this leads to a set of three equations: With respect to Eq. (9), the change in the indices of refraction can be calculated for stress state σ ij using stressoptical tensor π ij kl . Because the impermeability and the stress tensor are both symmetric second-order tensors, the stress-optical tensor π ij kl also has to show certain symmetries (Narasimhamurty, 1981). This allows it to be written as a 6 × 6 matrix in Voigt notation. In this form, indices ij and kl of the tensor π ij kl are reduced to 11 → 1, 22 → 2, 33 → 3, 23 → 4, 13 → 5 and 12 → 6. In the following, Voigt notation will be indicated by (V) . Due to the diamond structure of monocrystalline silicon, there are only three independent parameters of the stressoptical tensor (Narasimhamurty, 1981): π 11 π 12 π 12 0 0 0 π 12 π 11 π 12 0 0 0 π 12 π 12 π 11 0 0 0 Attention has to be paid to accounting for the appropriate coefficients while converting from tensor notation to Voigt notation (Narasimhamurty, 1981). Accordingly, the stressoptical parameters π (V) ij in Voigt notation are related to the stress-optical coefficients π ij kl in tensor notation by the following equation: In order to account for different crystalline orientations of silicon, a rotation matrix R ij is introduced. The rotated stressoptical tensor π ij kl accounts for different orientations of the silicon lattice structure by applying a rotation matrix R ij to it as follows: For simplicity, in the following only mechanical stresses inplane with the (100) plane of silicon are discussed. The [100] direction is further assumed to be parallel to the z axis of the coordinate system. A situation in which these simplifications arise is represented by a (100)-silicon wafer, as shown in Fig. 3. In this case, rotation matrix R ij (φ) describes a clockwise rotation around the z axis by an angle φ, which is the angle between the [010] direction and the x axis: As the (100)-silicon wafer is considered sufficiently thin, only plane stresses are evaluated: Inserting rotation matrix R ij (φ) into Eq. (18) and applying both to Eq. (9) yields the change in the impermeability tensor B ij for a plane stress state: If the incident light ray falling onto the indicatrix is parallel to the z axis of the chosen coordinate system, Eq. (15c) can be neglected. Therefore, the change in the indices of refraction can be expressed by subtracting Eq. (15b) from Eq. (15a). Further, the natural birefringence is comparably small against stress-induced birefringence, meaning that n o ij = n o (Krüger et al., 2016), yielding the following: To obtain the principal values B I and B II , a simple eigenvalue analysis on the impermeability tensor of Eq. (21) can be performed: in which a = 2 (π 11 + π 12 ) (σ 11 + σ 22 ) ,  or, in more detail, n = (n o ) 3 4 2 (π 11 − π 12 ) 2 + π 2 44 4σ 2 12 + (σ 11 − σ 22 ) 2 + 2 (π 11 − π 12 ) 2 − π 2 44 −4σ 2 12 + (σ 11 − σ 22 ) 2 cos 4φ From this equation it can be seen that there are only two material-dependent parameters, as the difference between (π 11 − π 12 ) cannot be separated.

Measurement of stress-optical coefficients
In order to measure the coefficients of π ij kl , a Brazilian disc test was used (as depicted in Fig. 5). The stress-induced retardation was measured using a grey field polariscope, as shown in Fig. 6. In parallel, mechanical stresses σ 11 , σ 22 and σ 12 were determined using the finite element method with Ansys Mechanical APDL. This allows for a full mapping of mechanical stresses onto the measured retardation values with consideration of the anisotropic mechanical response of a (100)-silicon wafer. Retardation values and stresses were correlated using Eq. (26) to determine the parameters (π 11 − π 12 ) and π 44 . The grey field polariscope uses circularly polarized light in the near-infrared spectral range. With this device, retardation values up to a quarter of the wavelength utilized can be measured without manual determination of fringe values. The measurement principle is described in detail elsewhere (Horn et al., 2005). Linked to the polariscope is a computer that calculates the retardation and the orientation of the optical axes from measured light intensities and saves both in ASCII format.
The Brazilian disc test consists of two diametrically placed clamping jaws that load the wafer on its thin sides. Force on the jaws is manually applied by a screw and is measured by a load cell. To mitigate high contact stresses, several layers of adhesive tape were placed in the contact areas of the jaws. This leads to a reduction in the applied force due to the relaxation of the tape material, which was taken into account by waiting for about 5 min before each measurement. For all measurements, the difference between start and end force was smaller than 0.1 N for an applied force of 10 N. This load is sufficient for a distinct measurable photoelastic signal without buckling the wafer. Buckling was observed to start at approx. 20 N.
The wafer was manually inserted between the clamping jaws. In total, 28 measurements were carried out for angles from −60 to 60 • with respect to the [011]-wafer direction. Angles were varied in increments of approx. 2.5 • . The orientation of the wafer was measured graphically from video data from the polariscope as the angle between the clamping jaw edge and the wafer flat.
For the finite element simulation of the Brazilian disc test, 19 200 elements with a quadratic displacement function were used in a linear elastic simulation. For each measurement, the mechanical stresses were simulated considering the specific applied load and the wafer orientation. Figure 7 shows the measured retardation values and the simulated stresses as an example of loading along the [011] direction of the wafer. From these, mechanical stresses and retardation measurements were taken to determine the stress-optical coefficients using a least-squares fit algorithm for Eq. (27). The stress-optical coefficients were determined to be (π 11 − π 12 ) = 14.4 · 10 −7 MPa −1 and π 44 = 9.4 · 10 −7 MPa −1 . The fitted model of Eq. (27) shows a coefficient of determination (R 2 ) of 0.81 for the angle-dependent stress-optical coefficients, and it is deemed to be a good fit. The levels of confidence were estimated by 10 repeated measurements for wafer orientations of 0 and 45 • . Confidence intervals for the measurement of wafer angles were assumed to be ±1 • at a 95 % confidence level. To estimate the confidence level of the fitted parameters, a Monte Carlo simulation for the non-linear fit of Eq. (27) was carried out. In    27) is shown for a normalized wafer thickness t and a normalized stress difference (σ I − σ II ) with measured stress-optical coefficients. Also, the angle-dependent isotropic coefficients C(φ) based on Eq. (3) for each measurement are plotted. Confidence intervals for the determined wafer orientations are smaller then the plotted dot in Fig. 8 and are therefore omitted.

Discussion
In comparison to models given in the literature, the stressoptical law derived in this work shows a close resemblance in shape to the work of Liang et al. (1992) and He et al. (2004). This can be seen in Fig. 9 for a normalized stress difference and a normalized wafer thickness. The work by Zheng and Danyluk (2002) shows a higher degree of anisotropy then other models. Quantitatively, this work shows a stronger stress-optical effect then studies in the literature (Liang et al., 1992;Zheng and Danyluk, 2002;He et al., 2004). This is reflected by the determined stress-optical coefficients which are higher then previously reported values, as shown in Table 1. Table 1. Stress-optical coefficients for a (100) silicon from experiments in this study and from the literature.
The sign of the stress-optical parameters depends on the definition of the sign of tensile and compressive stresses. Tensile stresses are considered positive stresses in this work. However, due to the measurement principal, the sign of the measured retardation value can not be determined. Therefore, the stress-optical parameters are assumed to be positive. The signs of the stress-optical parameters from the literature in Table 1 are as reported in the studies listed.
Although the model has a similar shape to those in the studies by He et al. (2004) and Liang et al. (1992), the mathematical expressions differ between the models. This is despite the fact that all models stem from the same general approach for photoelasticity given by Neumann (1841) and Pockels (1889Pockels ( , 1906. Only for the isotropic case with parameters (π 44 = π 11 − π 12 ) are all models except for the model from Zheng and Danyluk (2002) identical. For this case, the models show the expected direction-independent stress-optical constant, as is shown in Fig. 10. The model from Zheng and Danyluk (2002) only becomes isotropic for a parameter combination of π 44 = 2(π 11 − π 12 ).
For more detailed comparison, the rotated stress-optical tensors π ij kl and π (V) ij derived in this work and in He et al. (2004), respectively, are identical, considering that the definition of the angle φ in He et al. (2004) is the angle between the Figure 10. Comparison of the model derived in this study to models from the literature for the isotropic case of π 44 = (π 11 − π 12 ).
[010] direction of the crystal lattice and the principal axis of the stress tensor. For Liang et al. (1992) and Zheng and Danyluk (2002), this correspondence could not be established. Additionally, Liang et al. (1992) derived the stress-optical law assuming that the direction of incident light is lying in plane to the plane stress. In the work of He et al. (2004), a simplification is made by stating that the principal direction of the stress tensor σ ij coincides with the principal direction of the impermeability tensor B ij . This allows for an easier theoretical determination of B I and B II . In practice, this means that at the principal direction of the impermeability tensor must first be measured (e.g. by measuring the isoclinic angle), and the stress-optical tensor π ij kl must then be rotated accordingly. However, this assumption leads to a deviation of up to 4 % (8 % using the parameters determined by He et al., 2004) by ignoring off-diagonal terms of B ij . This deviation does not occur upon determination of B I and B II in an eigenvalue analysis in this work. The model from He et al. (2004) for n can be reduced to an expression where only the difference in principal stresses σ I − σ II occurs. The model derived in this work cannot be reduced to that. For this model, n is always dependent on the stresses σ 11 , σ 22 and σ 12 . Figure 7 shows a measured retardation for a loading direction of φ = 45 • . Measurements with slight asymmetric re-tardation maps were also found. This can be caused by misalignment of the clamping jaws, i.e. the two opposite forces compressing the silicon wafer do not lay on the same axis. This can lead to a deviation in the stress field. However, the comparison of the different experimental set-up results showed no strong influence on the experimental results here. In general, numerical analyses indicated that even a small misalignment influences the values of π 44 more strongly than (π 11 − π 12 ).

Summary and conclusion
Using the approach from Neumann (1841) and Pockels (1889Pockels ( , 1906 for describing stress-induced birefringence, a model was derived for the stress-optical effect in silicon. This model contains two independent stress-optical parameters (π 11 −π 12 ) and π 44 . The stress-optical parameters were determined to be (π 11 −π 12 ) = 14.4·10 −7 MPa −1 and π 44 = 9.4 · 10 −7 MPa −1 . using a Brazilian disc test and a finite element simulation. In total, 46 measurements were used to determine the parameters by loading a (100)-silicon wafer under different angles.
In appearance, the model derived shows a similar shape to those from Liang et al. (1992) and He et al. (2004). However, they differ with respect to their mathematical expression. The derivation of the stress-optical law from He et al. (2004) could be reproduced in part, while a simplification made by He et al. (2004) that introduces a certain error was not necessary in the model in this work. It determines the stress-optical law for a (100)-silicon wafer and, hence, allows for a more precise characterization of the mechanical stresses.
The determined parameter values are fairly large compared with those from the literature. However, considering the different sources of error in the experiment and simulation, the values here tend to increase rather than decreasing toward a better quantitative match with literature parameters. To get a closer match, either the measured retardation values need to be higher or the simulated stresses need to be lower. Unfortunately, neither the experimental or numerical steps indicate an approach toward this behaviour. Hence, this mismatch and the application of the model derived are the focus of further research.
Code and data availability. The underlying measurement data and code used in this work can be requested from the authors if required.