Development and comparison of the techniques for solving the inverse problem in photoacoustic characterization of semiconductors

In this work, theoretically/mathematically simulated (TMS) model is presented for the photoacoustic (PA) frequency response of a semiconductor in a minimum volume PA cell. By analyzing of the TMS model, the influences of thermal diffusivity and linear coefficient of thermal expansion on silicon sample PA frequency response were investigated and two methods were developed for their estimation. The first one is a self consistent inverse procedure for solving the exponential problems of mathematical physics, based on regression. The second one, a well trained three-layer perceptron with back propagation, based upon theory of artificial neural networks, is developed and presented. These two inverse problem solving concepts are applied to thermo-elastic characterization of silicon, compared and discussed in the domain of semiconductor characterization.


Introduction
Photoacoustics is a non-destructive experimental technique, being extensively developed over the last forty years and finding its application in the characterization and imaging of various materials and devices (Bialkowski 1996;Tam 1986;Vargas and Miranda 1988) The oldest and the most widespread PA technique is frequency gas-microphone method (Galovic et al. 2014;Rosencwaig and Gerscho 1976;Rousset et al. 1983), which has proven to be very efficient in determination of optical, thermal and elastic properties of semiconductors, semiconductor structures, metals, composite materials, ceramics, polymers, etc. (Balderas-Lopez and Mandelis 2001;Mansanares et al. 1991;Markushev et al. 2018;Nesic et al. 2018;Todorović and Nikolic 2000).
In frequency photoacoustics, the sample is illuminated by a monochromatic EM beam of harmonically modulated intensity. Part of the excitation energy is absorbed in the sample, and part of this absorbed energy is converted into heat through various de-excitation and/or relaxation processes, thus generating a harmonic heat source with frequency equal to the modulation frequency of the excitation beam. This heat source introduces a periodical perturbation in thermodynamic state of the sample, causing periodical temperature variations at sample surfaces, periodical surface displacements and sample bending. Due to adiabatic state change, thermally and mechanically induced periodical pressure build-up is formed in the surrounding gas, enclosed in the PA cell, which is recorded by a microphone (Rosencwaig and Gerscho 1976;Rousset et al. 1983) The detected signal can be associated with optical, thermal and elastic properties of the sample by means of establishing a TMS model which accounts for all described physical processes, with sample properties being represented as model parameters. It is these parameters which are to be obtained from the recorded signal by means of inverse solving, which classifies photoacoustics into modeldependent experimental techniques. In order to obtain satisfactory level of accuracy in the determination of physical properties, in model-dependent techniques, it is of utter importance to develop an optimal TMS model. Herein, the theoretical model of frequency PA response in a minimum-volume PA cell relates amplitude and phase of the variation of the recorded acoustic signal to two optical properties of the sample (reflection and absorption coefficient) and to four thermoelastic properties (coefficients of thermal diffusion and thermal conduction, linear coefficient of thermal expansion and heat velocity), while in semiconductors there are five additional properties to follow (energy gap, ambipolar diffusion coefficient, lifetime of photo-generated carriers and velocities of surface generation and recombination). Basing the research on PA frequency measurements, up to eleven physical properties of the material can be assessed! However, the listed properties affect PA response in very complex and non-linear manner, and, besides that, each influence is dependent on modulation frequency. In addition, small variations in model parameter values result in large variations of the model solution, meaning that inverse PA problem is non-linear, ill-conditioned problem of mathematical physics (Power 2002). With this in mind, this task needs to be carefully assessed from a multitude of aspects.
The analysis of the model for various materials-metals, semiconductors and polymers-has indicated different physical processes to be dominant in the experimentally obtainable frequency range. However, in this case, the influence of thermal diffusivity and linear coefficient of thermal expansion are the ones which affect the PA response the most.
In this paper, a SCIP is developed for thermoelastic characterization of semiconductors, based on previous experience in polymer characterization (Nesic et al. 2019(Nesic et al. , 2018, with overall good results in a limited measurement frequency range, i.e. in the range where the influence of the measurement chain can be neglected, where only one PA response component is dominant. Besides that, an ANN based procedure is developed, tailored to a limited measurement frequency range, based on previous experience in semiconductor characterization performed in the whole frequency range and on measurements previously preprocessed in order to be deprived of measurement chain influence (Djordjevic et al. 2020b(Djordjevic et al. , 2019. Finally, the results of these two approaches are compared. The aim of this research is to discover the methodology inverse PA problem solving with maximum precision, reliability and applicability on various samples and in various transmission PA experiments, with different microphones and other measurement chain elements.
In Sect. 2, the model of semiconductor frequency PA response is presented and its sensitivity to the variations of thermal diffusivity and linear expansion coefficient is discussed. In Sect. 3, both techniques of inverse problem solving are described: the self-consistent procedure and the neural network. In Sect. 4, the results of both techniques are presented and discussed. Finally, in Sect. 5, the most important conclusions are drawn.

TMS modeling-direct problem solving
The core of the experimental set-up consists of the light source, a minimum volume PA cell, and the microphone (Rabasovic et al. 2009). The light source is either a LED or a laser diode, its power ranging from 0.8 to 8 mW and its wavelength being 660 nm, intensity-modulated by a driver. This way produced electromagnetic radiation is projected on a sample, which is mounted on a minimum volume PA cell in the manner of a simple supported plate. The cell consists of a microphone and a photodiode, both installed in an Aluminum shield-quite inexpensive concept promising good signal-to-noise ratio (Balderas-Lopez and Mandelis 2001;Perondi and Miranda 1987). The electret microphone diameter (ECM60) is 9.8 mm, its sensitivity being 2.5 mVPa −1 ; the light source is intensity modulated and frequency-dependent measurements of PA response are performed using a lockin amplifier. A schematic representation of such a transmission configuration is given in Fig. 1: The measured open-cell PA response S exp (f ) can be represented in the following form: where A exp (f ) and exp (f ) are its respective amplitude and phase, i is the complex unit, f is the modulation frequency, G (f ) is the instrumental complex transfer function and p total (f ) is the pressure variations signal. To find the sample properties, one has to fit S exp (f ) to the p total (f ) function, determined by the well-known composite piston theoretical model, which, for semiconductors as the n-type Si considered in this work, establishes that Todorović et al. 2013Todorović et al. , 2014: where and are thermal-diffusion p TD (f ) , thermoelastic p TE (f ) and plasmaelastic p PE (f ) components of the PA signal. Here, γ is the adiabatic ratio, T 0 and p 0 are the ambient temperature and pressure, respectively, V 0 and l c are the PA cell volume and its length, l is the sample thickness, R s is the radius of the sample D g is thermal diffusivity of the gas in the PA cell, d n is the coefficient of electronic deformation, and ñ p and ̃s are complex representatives of the density of photogenerated excess charge carriers and of the temperature variations, respectively.
The complex representative of the dynamic component of the photo-generated minority carrier density across the sample ñ p (x) is given by: (1 + i ) , and where L p = √ D p τ is the excess carrier (holes) diffusion length, characterized by their diffusion coefficient D p and lifetime τ , is optical absorption coefficient and the integration constants Ã ± are defined by: According to Eqs. (1)-(11), the adjustment of S exp (f ) to p total (f ) can lead to the determination of the sample properties, such as thermal diffusivity D T , the coefficient of thermal expansion a T , etc.
In Fig. 2, amplitude and phase data of PA response of Silicon are presented (experiment and theory): Comparing the prediction of the abovementioned model to the experimental measurements on 150 µm thick Silicon sample, it is shown that TMS model matches the experiment only in the (50 Hz-2 kHz) frequency range. Notable difference at low and at high frequencies is attributed to the influence of the noise and to the influence of the measurement chain, so the utilization of the whole frequency measurement range demands detailed analysis and the elimination of these influences (Djordjevic et al. 2020a;Jordovic-Pavlovic et al. 2020;Markushev et al. 2020;Pavlović et al. 2019).
In this paper, inverse solving procedures are developed on a limited frequency range, aiming at the determination of silicon properties with dominant influence in the mentioned range, namely the coefficient of thermal diffusion and linear thermal expansion coefficient. Literature values of used parameters are, among others, given in Table 1.
As can be seen in Fig. 2a, the amplitude of PA response in the range (50-2000 Hz) changes several orders of magnitude of its value. In order to eliminate insufficiently known optical excitation and amplification influences, signal normalization at a thinner sample of the same kind, previously well characterized, is introduced (Balderas-Lopez and Mandelis 2003).
In Fig. 3, normalized amplitude and phase data is presented for the D T ± 0.1D T of thermal diffusivity change, where D T = 9 × 10 -5 m 2 s −1 .
In Fig. 4, normalized amplitude and phase data is presented for the a T ± 0.1a T of linear thermal expansion coefficient change, where a T = 2.6 × 10 −6 K −1 .
As can be seen from Fig. 3, thermal diffusivity affects both amplitude and phase characteristics of PA response in the whole frequency range. The influence of the coefficient of thermal expansion is more noticeable in amplitude characteristics and at high frequencies, as can be observed in Fig. 4. Contrary to the case with polymers, where the dominance of the PA component (Eq. 4) in the observed frequency range results in linear coefficient of thermal expansion not affecting the normalized PA response, in the case of semiconductors, in this frequency range, this influence is visible, especially in amplitude characteristics.

Self-consistent inverse procedure-SCIP
Analytical considerations of the model revealed that the phase of the PA response exhibits stronger dependence upon the parameter of thermal diffusivity, D T , and thus it can be estimated only from this part of the PA response with satisfactory accuracy. After the estimation, D T is put into the expression for the amplitude response thus improving the reliability of fitting the coefficient of thermal expansion.
In brief, the SCIP consists of two iterative single-parameter fitting processes: the first one aimed at thermal diffusivity determination and the second one, concerning the estimation of linear coefficient of thermal expansion (see Fig. 5).
In each iterative process, the parameter value is taken from its literature range and it is varied by the addition (or the subtraction) of 0.1% of its initial value, calculating the corresponding PA response in each iteration and comparing it to the corresponding PA response of the numerical experiment. The numerical experiment is obtained by the calculation of amplitude and phase characteristics in the limited frequency range (ROI) from 50 Hz to 2 kHz (Fig. 2) with preset values of D T and a T randomly taken from the corresponding literature range, given in Table 1. The calculated values are normalized to amplitude and phase characteristics of a well characterized 150 µm thick semiconductor sample ( D T = 9.0 × 10 −5 m 2 s −1 and a T = 2.6 × 10 −6 K −1 ). The comparison method is based on simple MSE computation scaled on the ROI. These iterations are repeated until one of the conditions is met: either a satisfactory match is achieved between the mentioned preset PA response and the calculated one (i.e. minimum value of the SME on the observed ROI), or the literature range boundary value is reached (security check). Then, the parameter value from the ongoing iteration is recorded as the closest one to the value corresponding to the preset PA response and it is considered the SCIP result. At the same time, the relative error for each parameter is calculated (in percentage), i.e. as the difference of value obtained by the SCIP iterations and the correct (i.e. preset) value, divided by the preset value. These error values represent the characteristic of the SCIP.
The self-consistency of the inverse procedure is reflected in the fact that, firstly, D T is extracted from the phase data of the PA response, and, secondly, this result is utilized in the calculation process of the determination of a T from the amplitude data of the PA response.

SCIP results
In Fig. 6, the results of the SCIP after thirty numerical experiments are presented. The results are sorted in ascending order of preset D T (Fig. 6, up) and preset a T (Fig. 6, down). The preset values are presented in red. Alongside, corresponding estimates are presented in green (400 µm thick silicon PA response normalized on the PA response at 150 µm) and in blue (900 µm thick silicon PA response normalized on the PA response at 150 µm). As can be noted in Fig. 6a, the estimated values of D T in the case of 400/150 normalization (green dots) are closer to the expected values of the numerical experiment (red dots, sorted in ascending order) in comparison to the ones based on 900/150 normalization (blue dots). With regard to theoretical considerations, presented in Fig. 3, this result seem to be unexpected, since normalized phase difference is lower in thinner samples (Fig. 4b), and the accuracy of D T results is expected to be lower in thinner samples.
In Fig. 6b the accuracy of a T estimates is also better in thinner samples. This result is also unexpected, based on theoretical considerations presented in Fig. 5, since the difference in normalized amplitudes is larger in thick samples (Fig. 4b).
In Fig. 7 relative estimation errors gained in numerical experiments are presented as the function of preset D T (up) and a T (down) values, sorted in ascending order, where green lines are indicative of 400/150 normalization, while blue lines are for the 900/150 one. It can be noted that relative errors of the procedure do not exert any strong dependency on the preset D T and a T values, although slight increase is observed in a T relative error for high values of preset a T , (blue line in Fig. 7b).
Maximum and mean values of relative error, regarded as the characteristics of the SCIP and gained after the performance of 30 numerical experiments, are presented in Table 2.
The analysis of Table 2 shows that maximum as well as mean values of relative error are lower for both D T and a T presets in thinner samples (below 1%), while the increase of sample thickness brings these values to around 2.5%. Also, the corresponding error values are lower in the case of D T in comparison to the case of a T estimation, i.e. the developed procedure estimates D T with greater precision then a T .

Inverse solution of photoacoustic problem by neural networks
Solving the inverse problem of photoacoustics in frequency domain using some normalization procedure still attracts considerable interest. Especially for the silicon samples of n-type which have the widest industrial application. That is why the necessity to find the procedure of quick and precise experimental determination of silicon individual properties applying some normalization procedure is an ultimate goal which can be achieved, besides the other, by using the artificial neural networks (ANN).
Using regression as a supervised learning technique with simple structured ANN one can solve the inverse problem estimating typical semiconductor parameters from the normalized photoacoustic signals as predicted output values of the network (for example, thermal diffusivity, thermal expansion, thickness, etc.).
Here, neural networks of the simple structure were formed with one input, one hidden and one output layer, represented in Fig. 8. The data used in the input layer are the amplitude-phase characteristics of the normalized photoacoustic signal magnitude defined with 40 modulation frequency values each, so 80 neurons establish the input layer. The output layer consists of two neurons, corresponding to the sample thermal parameters each, thermal diffusivity D T and thermal expansion a T . The hidden layer is formed by 50 neurons.  Neural networks presented here are trained in the specific modulation frequency range from 50 Hz to 2 kHz, using a large base of 284 theoretical photoacoustic signals obtained with Eqs. (2-5) for sample thicknesses of 400 µm and 900 µm, changing the thermal diffusivity and thermal expansion values within the range of D T = (8.1 − 9.9) × 10 −5 m 2 s −1 and a T = (2.34 − 2.86) × 10 −6 K −1 , respectively.
The amplitudes of the input signals are normalized with the function 20 log A i A i (150 μm) , i = 1,2,3..,40. In such a way their values are reduced to the values that correspond to the changes of the phase, taking the amplitude of 150 µm thick sample as a referent one with constant D T = 9.0 × 10 −5 m 2 s −1 and a T = 2.6 × 10 −6 K −1 . Phases are normalized with the function i − i (150 μm) , i = 1,2,3..,40, taking the same referent sample thickness.
The results of the network training using back propagation method for both levels of sample thickness of 400 µm and 900 µm normalized with 150 µm sample thickness are shown in Fig. 9. The first neural network ANN(400/150) was trained in 297 epochs activating the termination criteria with the performance of 5.7857 × 10 −8 (Fig. 9a). The second neural network ANN(900/150) was trained in 334 epochs activating the termination criteria with the performance of 2.9984 × 10 −8 (Fig. 9b).

Results of the neural network predictions
The results of the independent test with 5 signals from the large base, which are not used in network training, were presented for both networks in the Tables 3 and 4, with % relative prediction error (maximum and mean) obtained as a difference between PA signal parameters D T an a T and their ANN prediction.  The results presented in Tables 3 and 4 indicate high network accuracy in the predictions of the sample parameters. Table 5 shows that the errors are much lower than 0,1%, but even that is enough to confirm the theoretical model of photoacoustics used for their training. Namely, the errors for D T are lower than the errors for a T in both cases. Also, the errors for D T decrease while the errors for a T increase with decreasing thickness. This is in accordance with the theoretical model of photoacoustics (Eqs. 2-5) which predicts the dominance of thermal diffusion component of the signal (used to determine D T ) for thinner samples in the entire  Hz frequency range (other signal components are negligible).

Discussion and comparison of SCIP and ANN results
The characteristics of the two presented procedures aimed at inverse solving of PA problem and the estimation of two unknown semiconductor properties, in the frequency range limited to the part where measurement chain influence can be neglected, are described by maximum and mean relative error values, given in Tables 2 and 5. From these tables, it is clearly visible that ANN based procedure is characterized by a hundred (100) times lower level of error. This result can be explained by the fact that certain approximations, introduced in the SCIP in order to lower its processing time and based on previous theoretical considerations, are rough and seem to introduce significant level of error. The most illustrative one is the assumption that the influence of a T on phase characteristics can be neglected, since it is barely noticeable and present at high frequencies only. This approximation increases the error in a T estimation by considering that D T is estimated accurately enough in the first iterative process of the SCIP.
Another difference in methodology between these procedures is the fact that SCIP processes amplitude and phase data separately, while ANN-based procedure processes them simultaneously, effectively discovering (D T , a T ) pairs which result in the most similar PA response amplitude and phase at the same time.

Conclusions
The determination of two or more material properties from a single PA response measurement represents an important problem, the solving of which should result in a broader application of photoacoustics in various areas of science and technology. In this paper, two procedures aimed at solving this kind of problem are presented. Characterization parameters, namely maximum and mean determination relative error, are also proposed. SCIP is characterized by error values lower than 3%, making this approach a satisfactory choice in the case of rough estimation of D T and a T measured properties. ANN approach is characterized by error values under 0, 1%. The result this good establishes it as the method of choice in very accurate D T and a T determination from PA measurements, especially in the case where the properties obtained this way are used in further sample characterization, allowing accompanying procedures to bring in their determination errors and still keeping the obtained result within 1% uncertainty boundary.
The comparison of the presented procedures in terms of maximum relative error criterion shows that ANN based methodology is significantly better that the SCIP. The differences in methodologies are discussed and possible causes of significant errors in the SCIP are reported.
In future explorations, the developed and presented ANN-based procedure can be used as reference method in inverse solving of the PA problem while the SCIP can and should be modified in order to significantly lower the level of relative error.