# Hybrid metrology investigation combining Raman and ellipsometry spectroscopy applied to in line GeSbTe crystallization measurements and deep learning approaches for accurate prediction

Jon de Vecchy, Jean-Hervé Tortai, Maxime Besacier, Delphine Le Cunff, Bernard Pelissier

► To cite this version:

Jon de Vecchy, Jean-Hervé Tortai, Maxime Besacier, Delphine Le Cunff, Bernard Pelissier. Hybrid metrology investigation combining Raman and ellipsometry spectroscopy applied to in line GeSbTe crystallization measurements and deep learning approaches for accurate prediction. Microelectronic Engineering, 2025, pp.112345. 10.1016/j.mee.2025.112345. hal-05025127

HAL Id: hal-05025127

https://cnrs.hal.science/hal-05025127v1

Submitted on 8 Apr 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

Public Domain
Hybrid metrology investigation combining Raman and
Ellipsometry spectroscopy applied to in line GeSbTe
crystallisation measurements and deep learning approaches for
accurate prediction

Jon de Vecchy<fnref n="1" /> a,*, Jean-Hervé Tortai a, Maxime Besacier a, Delphine Le Cunff b, B. Pelissier a,†

a Univ. Grenoble-Alpes, CNRS, CEA-LETI, Grenoble INP, LTM, Grenoble F-38054, France

b STMicroelectronics, Crolles Cedex 38920, France

**Corresponding authors**

E-mail addresses: jdevec@gmail.com (J. de Vecchy) and bernard.pelissier@cea.fr (B. Pelissier)

KEYWORDS

Metrology, Phase-Change Materials, Neuron Networks, Ellipsometry, Raman spectroscopy

HIGHLIGHTS

* Combination of model-less ellipsometry and Raman spectroscopy using Neural Networks

* An ellipsometry/Raman dataset is acquired in cleanroom in a range of 25-410°C

* Hybridization feasibility is made, although the lack of data hinders better results

* Adding synthetic data improves the crystallization rate prediction results

<footnote n="1">¹ Permanent address : Autorité de sûreté nucléaire et de radioprotection (ASNR), PSN-RES/SEREX/L2EC, F-13115, Saint-Paul-lez-Durance, France</footnote>
GRAPHICAL ABSTRACT

Hybrid metrology investigation combining Raman and Ellipsometry spectroscopy applied to in line GST crystallisation measurements and deep learning approaches for accurate prediction

1) Model-less ellipsometry spectra and Raman-extracted crystallinity rate hybridization using neuron networks. The experimental dataset used for the learning process has been acquired in cleanroom in the 25°C-410°C range.

<figure><img src="image_1.png" /></figure>

2) Ellipsometry spectra and crystallinity rate evolutions were independently modelled to generate synthetic datasets. Different datapoint distributions were tested: homogeneously distributed on the 25°C-410°C range (Homogeneous) or equally distributed before and after the beginning of the crystallization process (Balanced). Models were trained using these datasets.

3) The synthetic models were used to make predictions for the experimental dataset. Their prediction results in the 390-410°C range overpass the models' results trained with this dataset. The best models were trained with the synthetic balanced dataset.

<figure><img src="image_2.png" /></figure>

1. ABSTRACT

To monitor the Ge-rich GeSbTe crystallization process, ellipsometry and Raman spectroscopy were correlated by Machine Learning using a Neural Network approach. Ellipsometry was selected for being a fast, non-destructive, and in-line metrology technique and Raman spectroscopy was selected for its crystallization monitoring potential.

An experimental ellipsometry/Raman spectroscopy dataset was acquired in a 25-410°C range. Assuming that the crystallization process is germanium-driven, the crystallinity rate was extracted from fitting the Raman germanium-related modes. Neural Network hybridization was performed to predict the crystallinity rate (output) from the raw ellipsometry spectra (input).

Models trained with these experimental data show poor performance, especially in the 390-410°C crystallization range. The lack of data was identified to be the main issue. To generate data, the experimental data were independently modeled using numeric temperature laws. Ellipsometry spectra were then generated and labeled with crystallinity rates at any given temperature. The synthetic datasets were then used as a training dataset, leading to a better prediction of the crystallization value, dividing the models' mean squared error by more than ten times.

Finally, the synthetic data-trained models were compared to the experimental data-trained models on an experimental test set. Synthetic data-trained models showed better performance in the crystallization range than the experimental data-trained models (~2 times lower mean squared prediction errors). The proof of concept of this study is thus validated and could lead to promising potential results in fast crystallization rate prediction of phase change material simply using optical experimental raw measurements.
## 2. INTRODUCTION

Random Access Memories (RAM) are known to be potential successors of widely used non-volatile flash memory devices, that reach their scaling limits [1, 2, 3]. Unlike flash memories, which operate with charge storage using a floating grid, these RAMs are based on their inner materials property variations. Among RAM types, Phase Change RAM (PCRAM) embeds a material with highly contrasted electrical or optical properties between its amorphous and its crystalline state. The reversible and fast transition between these two states can then be converted into 0 or 1 logic state writing and reading [1].

GeSbTe (GST) particularly fulfill these requirements [4, 5]. GST can be declined in different stoichiometries, one of the most studied and already widely used being Ge<sub>2</sub>Sb<sub>2</sub>Te<sub>5</sub> (GST 225) [6, 7, 8, 9]. This system crystallizes in Face-Centered Cubic phase (FCC) at around 150°C and shifts from FCC to hexagonal structure at 300°C [6, 10]. However, as the operating temperature of the device may be close to the amorphous to FCC structure transition, GST 225 is not the optimal candidate. Therefore, in this study, “Ge-rich” GST was used (compared to the 225 stoichiometry). Its crystallization temperature is above 380°C [11] and makes it a more promising candidate than the GST 225 [12].

The production of high-quality GST-based non-volatile memories needs reliable metrology support to monitor GST key parameters, such as the crystallinity rate of the deposited active layer. The metrology techniques must also be based on non-destructive, fast, and widely used in-line tools, to avoid breaking the process continuity. Based on these criteria, ellipsometry is well suited and has been selected as a starting point for this study. However, this technique does not give direct access to the GST crystallinity rate. Such material structure information can be obtained Raman spectroscopy, but this technique is not widely implemented as an in-line metrology method. The aim of this study is therefore to combine ellipsometry ease of in-line use with Raman spectroscopy crystallization monitoring. To do so, raw and model-less ellipsometry spectra are correlated with crystallinity rate extracted from Raman spectra using Machine Learning Neural Network (NN). Such model would allow to perform automated quality control and predict crystallinity rate of a deposited GST layer only using ellipsometry, without needing any stacking model. This full spectral response processing approach can be considered as a proof-of-concept that can potentially be adapted to more characterization techniques.

To train the model, a self-made dataset consisting of raw ellipsometry and Raman measurements is acquired in a large temperature range (i.e. at different GST crystallization rates). The first part of the results presents the GST crystallization dataset building. The second part presents NN hybridization tests based on this dataset (raw and enriched).

## 3. MATERIAL AND METHODS

Two identical TiN-capped Ge-rich GST active layers were deposited on SiO<sub>2</sub>/Si 300 mm wafers at ST Microelectronics facilities. The upper TiN layer allows the electrical contact of the memory cell and plays a role to prevent oxidation of the GST active layer. Dataset was acquired using a specific 300mm setup, developed at CNRS-LTM in the scope of the IMPACT EquipEx program. Among others, the platform embedded a Raman spectrometer and an UVISEL2 ellipsometer equipped with an in-situ heating system (see Fig. 1). Automated vacuum transfers are made through a Pfeiffer Interface, compatible with 300 mm wafers using a robot. This system prevents wafers from air exposure during the whole study (heating, ellipsometry and Raman measurements).
<figure><img src="image_3.png" /><figcaption>Fig. 1: Raman spectrometer and ellipsometer cluster used for dataset acquisition. The whole cluster is under vacuum, and an in-situ heating system is included in the ellipsometer chamber.</figcaption></figure>

33 temperature data points were selected according to internal preliminary results and are described in Fig. 2. The first wafer was used to acquire the data points from 25 to 380°C, 390°C, 397.5°C and from 400 to 410°C. The second wafer was used to acquire more data points around the transition temperature range, i.e. 385°C, 392°C, 394°C, 396°C and 398°C.

<figure><img src="image_4.png" /><figcaption>Fig. 2: Chosen temperatures (in °C) at which the data points were measured for the dataset acquisition using two different wafers.</figcaption></figure>

Each data point is acquired following this sequence: (1) Heating to chosen temperature in the ellipsometry chamber under partial pressure (1 mbar of N₂) at 20°C/min with 120 seconds of stabilization. (2) ellipsometry measurement in the 0.6 to 3 eV energy range with 0.025 eV steps. (3) Natural cooling down in the ellipsometry chamber. (4) Vacuum transfer to Raman chamber. (5) Raman spectroscopy measurement recorded from 70 to 400 cm⁻¹ with a 532 nm green laser. The specified laser power is 20 mW but was kept at 5% to prevent any induced crystallization. Acquisition time was chosen to maintain the same intensity level and two accumulations were recorded to prevent energy spikes. (6) Vacuum transfer to ellipsometry chamber to heat the wafer at the next temperature of the ramp.

The ellipsometry parameters used in this study are *Is* and *Ic*. They are linked to the amplitude ratio *ψ* and the phase difference *Δ* by Eq. 1 and Eq. 2:

$$Is = \sin(2\Psi) \sin(\Delta) \qquad \text{Eq. 1}$$

$$Ic = \sin(2\Psi) \cos(\Delta) \qquad \text{Eq. 2}$$

Raman spectra were smoothed using a Savitzky-Golay filter directly implemented in Python's Scipy package. A baseline was also used. We selected a baseline suited for Raman spectra and that can be implemented directly on Python code [13].
# 4. RESULTS AND DISCUSSION

## 4.1. Building the GST crystallization experimental dataset

Ellipsometry and Raman datasets are shown in Fig. 3 and Fig. 4 respectively. A focus on the region of interest, where the GST transition is clearly visible, is also shown for each spectrum. Although easily detectable, the precise variation of raw ellipsometry Is and Ic spectra as a function of the temperature appears to be non-linear. A NN model would then allow to correlate every Is and Ic datapoints variation with crystallinity rate evolution, which may be tedious to perform using another method.

<figure><img src="image_5.png" /></figure>

Fig. 3: Raw ellipsometry (a) Is and (c) Ic dataset, with (b,d) focus between 0.7 and 1.2 eV. Temperature datapoints highlighting the start of the GST transition are marked (390-394°C). The same colorbar is applied to all graphs.
<figure><img src="image_6.png" /><figcaption>Fig. 4: (a) Raman spectra dataset. Vibration modes are highlighted. All spectra are smoothed and standardized to GeTe mode at around 80 cm⁻¹. (b) A focus on the c-Ge mode. Temperature datapoints highlighting the start of the GST transition are marked (390-394°C). The same colorbar is applied to both graphs.</figcaption></figure>

Five modes can be identified on the Raman spectra. The first one, at around 80 cm⁻¹ is assigned to E mode of GeTe₄ tetrahedra [14, 15]. To compare the evolution of each mode with temperature, all spectra are standardized to this mode as its intensity variation with temperature is small. At early temperature stages, when GST is still amorphous, a second mode is identified at around 150 cm⁻¹, assigned to A²₁g of SbTe₃ pyramids [16]. These two modes are related to the GST structure and are not specific to Ge-rich stoichiometries. Two additional modes are identified at around 220 and 270 cm⁻¹. These modes were assigned to amorphous germanium Longitudinal Optical (LO) and Transversal Optical (TO) vibrations [17]. In Ge-rich GST, germanium in excess tends to segregate from the GST structure [11], which can then approach a more stable stoichiometry before the crystallization occurs.

A transition starts above 385°C and can be seen in both ellipsometry and Raman datasets. The good correlation between both techniques is thus promising for NN hybridization. In the Raman spectra, firstly, the SbTe₃ vibration mode gets more intense (relatively to the GeTe 80 cm⁻¹ mode). This might be surprising considering that GST crystallizes and that all related-mode intensities should increase with structural rearrangement. This can be explained by the fact that SbTe₃ mode has higher polarizability than the GeTe mode [15]. As the transition occurs, an additional mode appears at around 125 cm⁻¹ and is assigned to A₁ mode of GeTe₄·nGeₙ corner-sharing tetrahedra (n=1,2) [16, 17, 18]. This mode was observed in GST 225 and would support the assumption that the final GST stoichiometry is close to the latter when crystallization begins. Segregated germanium also evolves at this temperature: a well-defined mode appears at 300 cm⁻¹ and can be assigned to a Ge-Ge vibration mode in crystalline germanium [17]. At the same time, intensity of both amorphous germanium-related modes decreases, which shows crystallization of pure germanium. Since the germanium depletion and the GST crystallization are concomitant, we assume that the GST crystallization is germanium-driven and can be monitored by the evolution of the germanium-related Raman modes.

Based on the latter assumption, the crystallinity rate, which will be used as the output in the NN approach, is directly extracted from germanium-related modes (all located from 200 to 320 cm⁻¹). Three lorentzians were used, centred at 220 cm⁻¹, 270 cm⁻¹ and 300 cm⁻¹, to fit amorphous Ge LO mode (a-Ge₁), amorphous Ge TO mode (a-Ge₂) and crystalline Ge mode (c-Ge) respectively. Examples of fitted spectra recorded at temperature under and above the beginning of crystallization are shown in Fig. 5a. The integrated area ratio ($A_{cryst}$) between amorphous and crystalline-related modes was calculated using Eq. 3. Making the final assumption that crystallization is complete at 410°C, $A_{cryst}$ is directly converted into crystallinity rate ($\tau_{cryst}$), following Eq. 4. Calculated crystallinity rate in the transition temperature range is shown in Fig. 5b.
$$A_{cryst}(T) \approx \frac{A_{c-Ge}(T)}{A_{a-Ge1}(T) + A_{a-Ge2}(T)} \qquad \text{Eq. 3}$$

$$\tau_{cryst}(T) \approx 100 * \frac{A_{cryst}(T)}{A_{cryst}(410^{\circ}\text{C})} \qquad \text{Eq. 4}$$

<figure><img src="image_7.png" /><figcaption>Fig. 5: (a) Fitted Raman spectra taken at 390 and 392°C (i.e. before and after the beginning of crystallization). Amorphous Ge LO and TO modes are named "a-Ge₁" and "a-Ge₂" respectively. Crystalline Ge mode is named "c-Ge". (b) Crystal/amorphous integrated area ratio and direct conversion to crystallinity rate. (c) Crystallinity rate evolution as a function of the temperature.</figcaption></figure>

At this stage, a complete experimental dataset of 33 datapoints (recorded at 33 different temperatures) is built. Each datapoint contains Is & Ic raw spectra (inputs) associated with a crystallinity rate (output).

## 4.2. NN hybridization

The dataset was randomly separated into train and test splits (60% and 40% respectively). The split randomness was controlled by putting random state at 0. In a very standard way, the model is trained on the train set and tries to predict the values of the test set. The metric used is the Mean-Squared Error (MSE) as it highly penalizes large prediction errors. This choice counterbalances the fact that the biggest part of the output dataset is close to 0% (datapoints measured at T ≤ 390°C) and can lead to high R² despite poor predictions in the transition temperature range (datapoints measured at T > 390°C). A maximum 5% mean error tolerance is chosen for the crystallinity rate value, meaning that the MSE target value would not exceed 25. To really focus on the model performance inside the crystallization temperature range, another "crystallization" MSE was also calculated including only data points in the 390-410°C range. For visualization convenience, model results are shown as crystalline rates as a function of the temperature, even though the temperature is not directly considered during the learning process.
The data pre-processing method is shown in Fig. 6. *Is* and *Ic* spectra in their entirety are concatenated and are used as input variables to guess the crystallinity rate (the output). All variables were scaled (using StandardScaler). Different network architectures were tested using Keras Python package.

Model loss curves were used to check the learning process. They show the MSE of the model predictions for both train and test sets at each iteration (epoch). In that case, the MSE is calculated on scaled data, so it is lower than the MSE calculated for the final prediction after unscaling.

<figure><img src="image_8.png" /><figcaption>Fig. 6: Neural Network-based approach with (a) input data pre-processing and (b) global NN strategy.</figcaption></figure>

### 4.2.1. Results using the experimental dataset

The very first tested architecture is constructed as follows:

- The input layer (194 neurons – one for each concatenated *Is*+*Ic* energy point).

- Two hidden layers (first layer of 200 neurons and second layer of 100 neurons).

- The output layer (1 neuron – the crystallinity rate).

The output activation function was “linear”, and the hidden activation layer was “relu”. Adam optimizer was used with a learning rate fixed at 0.01. Fig. 7 shows the NN test set predictions.
<figure><img src="image_9.png" /><figcaption>Fig. 7: Neural Network results using two hidden layers of 200 and 100 neurons with preliminary data scaling (using StandardScaler). Random State was fixed at 0. a) Model loss curve on train and test set (using Mean Squared Error) and b) final test set predictions (unscaled plot for convenience). Unscaled MSE is also mentioned. The insert shows a focus on the crystallization temperature range.</figcaption></figure>

The model Loss curve decreases as a function of epochs for both train and test sets. This indicates a learning process, which shows that the correlation is possible. However, one must be careful with overfitting (because test set performances are poorer compared with train set ones) or underfitting. One can easily notice that the model predictions for datapoints in the crystallization range (390-410°C) are inaccurate and lead to a very high MSE. To check NN consistency, train/test split was made without using any random state value. The same NN architecture was kept. Results are not shown here but there is still an inconsistency, associated with MSEs in the 150-200 range at best.

Numerous architectures (two or three hidden layers) were also tested, again fixing the random state value to 0. Each architecture was tested five times in a row and the model was deleted after each run. MSE was calculated each time and a final mean MSE and the standard deviation were finally calculated at the end of the five runs. The architectures are defined as [N₁,N₂] or [N₁,N₂,N₃], where N₁, N₂ and N₃ are the number of neurons of the first, the second and the third hidden layer respectively. The input and the output layers remained the same as the very first tested architecture. Fig. 8 shows results using all experimental datapoints and considering only the 390-410°C datapoints.
<figure><img src="image_10.png" /><figcaption>Fig. 8: Model results fitting the raw experimental dataset. Different architectures are tested using two or three hidden layers including a variable number of neurons. The target MSE is mentioned in red. Light blue results are calculated for the complete dataset whereas dark blue results are calculated for 390-410°C datapoints only (i.e., inside the crystallization temperature range). In all cases, the error bar represents the standard deviation.</figcaption></figure>

None of the tested architectures can produce a model with prediction performances below the target MSE. The modeling in the crystallization temperature range is inaccurate, leading to MSE > 750 when considering only the 390-410°C datapoints. The lack of data may still explain the inconsistency. However, measuring new data points in the transition temperature range is experimentally extremely difficult to achieve because of the heating precision and the achievable temperature stability during ellipsometry measurement (there is a 2°C oscillation inside a temperature step).

Therefore, to overpass the experimental issues, data enrichment is performed by adding synthetic data to the dataset. Since the study is focused on ellipsometry and crystallinity rate correlation, it is necessary to simulate both separately.

### 4.2.2. Building the synthetic datasets

The aim of this section is to be able to independently generate the inputs (the ellipsometry spectra) and the output (the crystallinity rate). Our approach is to fit both as a function of the temperature. This will enrich the dataset with synthetic inputs/output generated at any desired temperature.

Crystallinity rate was defined using germanium-related Raman modes, as explained in Fig. 5 and Eq. 4. The temperature evolution of their integrated area was fitted using a sigmoid function defined as Eq. 5:

$$f(T) = \alpha + \frac{\gamma}{1 + e^{\delta(T - \varepsilon)}} \qquad \text{Eq. 5}$$

Where $\alpha$, $\gamma$, $\delta$ and $\varepsilon$ are constants to be optimized with the Python curve_fit method. Please note that this is a purely numerical approach. Fig. 9(a-c) shows all three fittings with their related expressions. As the c-Ge mode is the main contribution in the crystalline rate, it counterbalances the a-Ge modes poor fitting. Finally, the crystallinity rate evolution as a
function of the temperature was calculated and is shown in Fig. 9(d). Finally, a random noise is added to the final crystallinity rate model to consider numerical approximations for the machine learning process. Another random noise is added to the temperature value to consider experimental temperature measurement uncertainty.

<figure><img src="image_11.png" /><figcaption>Fig. 9: Fitting of (a) a-Ge 220 cm⁻¹ mode, (b) a-Ge 270 cm⁻¹ and (c) c-Ge mode as a function of temperature. The fit expression and associated MSE are also shown for each mode. (d) Final crystallinity rate ($\tau_{cryst}$) as a function of temperature. Its expression is also shown.</figcaption></figure>

For generating the ellipsometry spectra as a function of the temperature, one must look at the ellipsometry theory. Is and Ic spectra are generated from optical indices n and k of each layer of a chosen stack model. The optical indices are calculated from imaginary and real parts of permittivity $\epsilon$, generated using Tauc-Lorentz oscillators (with adjustable A, E, C parameters). Principle and equations are given in supplementary information (A) and are based on the Handbook of ellipsometry [19] and previous works [20, 21, 22]. GST bandgap values are around 0.45 eV and 0.20 eV for the amorphous and crystalline structure respectively. These values are slightly lower than values reported in the literature [23, 24], but as the present GST does not have the same stoichiometry, it is difficult to compare.

Once all spectra were fitted, the Tauc-Lorentz A, E and C parameters, GST bandgap, dielectric constant and thickness were extracted to model their evolution as a function of temperature. The pure numerical approach was used again, similarly to the Raman area peak modelling. A “biased” sigmoid function (Eq. 6) was used for each parameter:

$$f(T) = \alpha + \beta * T + \frac{\gamma}{1 + e^{\delta(T-\epsilon)}} \qquad \text{Eq. 6}$$

Where $\alpha$, $\beta$, $\gamma$, $\delta$ and $\epsilon$ are constants to be optimized with the curve_fit method. An example of a single Tauc-Lorentz oscillator fitting (A, E and C respectively strength, energy centre and peak broadening of the Tauc Lorentz oscillator) is shown in Fig. 10.
<figure><img src="image_12.png" /><figcaption>Fig. 10: Tauc-Lorentz oscillator parameters modelling as a function of the temperature. MSEs are also given.</figcaption></figure>

Having $A(T)$, $E(T)$, $C(T)$, $E_g(T)$ and $\epsilon_r(\infty)(T)$, synthetic $Is$ and $Ic$ spectra can be generated at any desired temperature using the method explained above. Thanks to the independent temperature evolution modelling of the inputs (the $Is$ and $Ic$ spectra) and the output (the crystallinity rate), a datapoint can be generated at any desired temperature. Two different datasets were generated and are shown in Fig. 11:

<figure><img src="image_13.png" /><figcaption>Fig. 11: Balanced (a) $Is$ and (b) $Ic$ and Homogeneous (c) $Is$ and (d) $Ic$ synthetic spectra.</figcaption></figure>

The “homogeneous” set is made of 5507 data points homogeneously distributed from 25°C to 410°C (with a 0.07°C temperature step). The “balanced” set is made of 4937 data points equally distributed before and after the GST transition (which means a 0.15°C temperature step in the 25-390°C range and a 0.008°C temperature step in the 390-410°C range). These datasets were finally duplicated with a temperature value noise. This was made for considering the experimental temperature measurement uncertainty. These synthetic datasets were then used for NN learning. For convenience, homogeneous and balanced set-trained models will be simply labeled homogeneous and balanced models respectively. Experimental set-trained models will also be labeled experimental models.
### 4.2.3. Results using the synthetic datasets

The previously described datasets were split into training and testing sets. The measured dataset was kept as a validation set and will only be used for the final model validation. Inputs and outputs were scaled using `StandardScaler`.

Examples of homogeneous and balanced test set predictions are shown in Fig. 12 and Fig. 13 respectively. An overfitting appears for both train and test sets, where model loss suddenly increases and is no longer improved with additional epochs (see Fig. 12(a)). This is corrected by adding a stopping function to the learning process: when the model loss is not improved after 10 consecutive epochs, the learning process stops, and the best neuron weights are kept.

<figure><img src="image_14.png" /><figcaption>Fig. 12: Neural Network results on "homogeneous" synthetic dataset using two hidden layers of 200 and 100 neurons with preliminary data scaling (using `StandardScaler`). Random State was fixed at 0. a) Model loss curve on train and test set (using Mean Squared Error) and b) final test set predictions (unscaled plot for convenience). Unscaled global MSE is also shown. The insert shows a focus on the crystallization temperature range.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Fig. 13: Neural Network results on "balanced" synthetic dataset using two hidden layers of 200 and 100 neurons with preliminary data scaling (using `StandardScaler`). Random State was fixed at 0. a) Model loss curve on train and test set (using Mean Squared Error) and b) final test set predictions (unscaled plot for convenience). Unscaled global MSE is also shown. The insert shows a focus on the crystallization temperature range.</figcaption></figure>

For both datasets, one can notice a significant performance improvement when compared with the experimental models: MSEs are approximatively ten times lower and even under the 25-target. The crystallization modelling must be improved though, especially around the starting and the ending of the crystallization. This observation is more striking for the homogeneous model. As for the experimental models, the contrast between good MSE and perfectible crystallization modelling is due to the excess of pre-crystallization data (at T<390°C), unbalancing the final MSE. The second MSE calculated using only the 390-410°C datapoints is thus kept for model evaluation.
The NN consistency is checked using numerous architectures for both datasets, by performing five runs per architecture, as already presented for experimental models. Fig. 14(a) and (b) show the homogeneous and balanced model results respectively (mean MSE and standard deviations). Results using all datapoints and only 390-410°C datapoints are once again discriminated.

<figure><img src="image_16.png" /><figcaption>Fig. 14: Model results fitting (a) the homogeneous synthetic dataset and (b) the balanced synthetic dataset. Different architectures are tested using two or three hidden layers including a variable number of neurons. The target MSE is mentioned in red. The error bars represent the standard deviation. For better comparison, the same y scale is kept for both graphs. Consequently, arrows indicate when a MSE is out of bond.</figcaption></figure>

For both datasets, considering all datapoints, on-target MSEs are reached multiple times using different neuron network structures, which is very promising, although it is difficult to conclude about finding an optimal neural network structure. Homogeneous models have overall better results. This trend is totally inverted when only considering the crystallization temperature range (390-410°C datapoints). In that case, balanced models have significant lower and less dispersed MSEs. As the same y axis was kept in Fig. 14(a) and (b) for comparing both datasets, these observations are graphically translated by the larger error bars and some out-of-window MSEs on the homogeneous models' side. Training with a balanced dataset is then more effective to model the crystallization process.

Finally, homogeneous and balanced NN models, statistically representative of each architecture results, were stored. These stored models were used for final validation predicting the raw experimental dataset.

### 4.2.4. Final model validation

Although building models from fully synthetic data seem promising, one must validate these synthetic models using the raw dataset and compare them with experimental models trained with these raw data. In this purpose, the same raw dataset train/test split was kept, and all synthetic models were used on the test set. They were compared to the best experimental models trained on this test set. All experimental models, homogeneous models and balanced models' results are shown on Fig. 15 a), b) and c) respectively. The same y axis was kept for easier graphical comparison. The overall synthetic models' performances are better than the experimental ones trained on this set, thus validating the promising improvements highlighted in the previous part. Once again, when comparing the synthetic models, the balanced models reach lower and less dispersed MSEs, as for the synthetic data case.
However, the overall results of the synthetic models on this experimental set are significantly less efficient than on the synthetic set, as any of the stored model can reach the target MSE. There is still a gap between the synthetic and the experimental data and modelling must then be improved, especially the numeric laws in temperature.

<figure><img src="image_17.png" /><figcaption>Fig. 15: Stored (a) experimental, (b) homogeneous, and (c) balanced models' predictions on the experimental test set. On each graph, results are shown considering all the test set and considering only datapoints in the crystallization temperature range (390-410°C). Models are classified by their neuron architecture. The target MSE is mentioned in red. For better comparison, the same y scale is kept for all graphs.</figcaption></figure>

For a better crystallization modelling visualization, models' predictions are also shown point by point on Fig. 16 a), b), and c) for the experimental models, the homogeneous models, and the balanced models respectively. Aggregated results are shown in Fig. 17.

<figure><img src="image_18.png" /><figcaption>Fig. 16: Experimental test set predictions using stored a) experimental, b) homogeneous, and c) balanced models. The solid lines indicate the mean models' performance (i.e. gathering all the neuron architectures results for each type of model). The coloured ranges are bound by the maximum and the minimum crystallinity rate prediction considering all the models' predictions for each datapoint. The inserts show a focus on the crystallization temperature range.</figcaption></figure>
<figure><img src="image_19.png" /><figcaption>Fig. 17: Test set predictions using all stored models on the a) whole experimental test set and b) focused on the crystallization temperature range. The solid lines indicate the mean models' performance (i.e. gathering all the different neuron architectures results for each type of model). The coloured range is bound by the maximum and the minimum crystallinity rate prediction considering all the same type models' predictions for each datapoint.</figcaption></figure>

<figure><img src="image_20.png" /></figure>

It appears a modelling contrast between pre-crystallization datapoints (25-390°C) and in-crystallization range datapoints (390-410°C). For the pre-crystallization datapoints, the synthetic models are less efficient than the experimental ones. This trend reverses for the in-crystallization range datapoints. In the latter case, balanced models are the most efficient models. This is statistically verified, as shown in Table 1, where all models' results are gathered into single MSE/StD global and crystallization-focused values, for each model type.

<table><caption>Table 1: Experimental test set predictions' results for each model type, gathered into single MSE/StD global and crystallization-focused values.</caption><thead><tr><th></th><th>Experimental</th><th>Homogeneous</th><th>Balanced</th></tr></thead><tbody><tr><td>Global MSE ± StD</td><td>226,7 ± 54,9</td><td>139,8 ± 45,6</td><td>95,7 ± 20,1</td></tr><tr><td>“Crystallization” MSE ± StD</td><td>618,5 ± 162,2</td><td>361,6 ± 129,5</td><td>175,2 ± 31,6</td></tr></tbody></table>

The overall modelling of the crystallization temperature range is still significantly dispersed and inaccurate, especially for 405-410°C datapoints, thus showing the need of more experimental datapoints in this range, although being technically challenging. This may improve the synthetic dataset building and allow to train even more efficient synthetic models and consistently reach the target MSE.

## 5. SUMMARY AND CONCLUSION

To make a proof of concept of a GeSbTe (GST) crystallization hybrid metrology using Raman spectroscopy and ellipsometry, a complete experimental dataset was acquired. 33 data points were quasi in-situ measured in a temperature range of 25-410°C.

The amorphous GST is a mixture of a pure GST phase and an amorphous germanium phase formed by germanium segregation. Crystallization process is dual: segregated germanium forms a pure crystalline germanium phase and Ge-depleted GST crystallizes aside. Assuming that the crystallization process is germanium driven, the crystallinity rate was extracted from Raman germanium-related modes using Lorentzian fitting and amorphous/crystalline ratio calculated from the integrated area.

Correlation between raw ellipsometry spectra and Raman-extracted crystallinity rates was achieved using Neural Networks. Train/test split randomness was controlled, and different hidden layer architectures were tested. The very first results show that a learning process is involved, highlighted by a decrease of both train set and test set Mean Squared Error (MSE) as a function of the epoch, proving the correlation feasibility. However, all tested
architectures suffer from inconsistency and poorly model the crystallization range (reaching
MSE > 250 in this range).

To check if the lack of data is the cause of the previous issues, the dataset was enriched with synthetic data. To do so, ellipsometry spectra were modeled using Tauc-Lorentz (TL) oscillators. A purely numeric law in temperature was used to generate TL oscillators and then ellipsometry spectra at any desired temperature. In parallel, the same numeric approach was used for Raman Lorentzian fitting parameters, also allowing to generate a crystallinity rate at any temperature. Two synthetic datasets were generated: a homogeneous dataset – made of data points equally distributed in the whole 25-410°C range – and a balanced dataset – made with the same amount of data points before (25-390°C range) and during the crystallization (390-410°C range).

Neural Network models trained with these synthetic datasets then led to promising results
with ~10 times lower MSEs than the raw data-trained models. Statistically relevant synthetic
models were stored for every tested neuron architecture and were used to predict crystallinity
rates of the experimental test set. Synthetic models were even more accurate than the
experimental models in the crystallization range. At the same time, they poorly predict pre-
crystallization points. The gap between raw and synthetic data, due to inaccurate laws in
temperature, may explain this behavior.

Despite this modelling issue, the synthetic model performance overpassing the experimental
one in the crystallization range is very promising. The hybrid metrology proof of concept can
be considered a success and could be applied to experimental cases with more subtle
spectral variations. It also validates the approach of modelling experimental datasets and
enriching them with synthetic data that are experimentally out of reach.

6. ACKNOWLEDGEMENTS

This work has been supported by the European ECSEL project "Madeln4".
It has also been partially supported by the French Government program "Investissements
d'Avenir" managed by the National Research Agency (ANR) under the contract number
ANR-10-EQPX-33 "IMPACT" Project, and by the French RENATECH network as well.
7. REFERENCES

[1] C. Zambelli, G. Navarro, V. Sousa, I. Prejbeanu and L. Perniola, "Phase Change and Magnetic Memories for Solid-State Drive Applications," Proceedings of the IEEE, vol. 105, no. 9, pp. 1790-1811, 2017.

[2] K. Takeuchi, "Scaling challenges of NAND flash memory and hybrid memory system with storage class memory & NAND flash memory," in Proceedings of the IEEE 2013 Custom Integrated Circuits Conference, San Jose, CA, USA, 2013.

[3] R. Annunziata, P. Zuliani, M. Borghi, G. De Sandre, L. Scotti, C. Prelini, M. Tosi, I. Tortorelli and F. Pellizzer, "Phase Change Memory technology for embedded non volatile memory applications for 90nm and beyond," in 2009 IEEE International Electron Devices Meeting (IEDM), Baltimore, MD, USA, 2009.

[4] E. Morales-Sánchez, E. Prokhorov, J. González-Hernández and A. Mendoza-Galván, "Structural, electric and kinetic parameters of ternary alloys of GeSbTe," Thin Solid Films, vol. 471, no. 1-2, pp. 243-247, 2005.

[5] M. Wuttig and N. Yamada, "Phase-change materials for rewriteable data storage," Nature Materials, vol. 6, pp. 824-832, 2007.

[6] C. Liu, Q. Tang, Y. Zheng, B. Zhang, J. Zhao, W. Song and Y. Cheng, "The origin of hexagonal phase and its evolution process in Ge2Sb2Te5 alloy," APL Materials, vol. 10, no. 2, p. 021102, 2022.

[7] Y. Zhang, Q. Duan, X. Yan, Q. Zhang and Y. Lu, "Nonvolatile multilevel adjustable optical switch based on plasmonic slot waveguide and GST segmented structure," Optics Express, vol. 32, no. 9, pp. 16548-16562, 2024.

[8] R. Duan, J. Sun, Y. Zhang, H. Li, Y. Li and Z. Liu, "A Non-volatile Quasi-Continuous All-Optical Fiber Programmable Platform Based on GST-Coated Microspheres," ACS Photonics, vol. 9, no. 4, pp. 1180-1187, 2022.

[9] P. Bartlett, S. Benjamin, C. (Kees) de Groot, A. Hector, R. Huang, A. Jolleys, G. Kissling, W. Levason, S. Pearce, G. Reid and Y. Wang, "Non-aqueous electrodeposition of functional semiconducting metal chalcogenides: Ge2Sb2Te5 phase change memory," Materials Horizons, vol. 2, pp. 420-426, 2015.

[10] P. Guo, A. Sarangan and I. Agha, "A Review of Germanium-Antimony-Telluride Phase Change Materials for Non-Volatile Memories and Optical Modulators," Applied Sciences, vol. 9, no. 3, p. 530, 2019.

[11] M. Agati, M. Vallet, S. Joulié, D. Benoit and A. Claverie, "Chemical phase segregation during the crystallization of Ge-rich GeSbTe alloys," Journal of Materials Chemistry C, vol. 7, no. 28, pp. 8720-8729, 2019.

[12] M. Luong, M. Agati, N. Ratel Ramond, J. Grisolia, Y. Le Friec, D. Benoit and A. Claverie, "On Some Unique Specificities of Ge-Rich GeSbTe Phase-Change Material Alloys for Nonvolatile Embedded-Memory Applications," Physica Status Solidi RRL, vol. 15, p. 2000471, 2021.

[13] V. Korepanov, "Asymmetric least-squares baseline algorithm with peak screening for automatic processing of the Raman spectra," Journal of Raman Spectroscopy, vol. 51, no. 10, pp. 2061-2065, 2020.

[14] R. De Bastiani, E. Carria, S. Gibilisco, A. Mio, C. Bongiorno, F. Piccinelli, M. Bettinelli, A. Pennisi, M. Grimaldi and E. Rimini, "Crystallization of ion amorphized Ge2Sb2Te5 thin films in presence of cubic or hexagonal phase," Journal of Applied Physics, vol. 107, no. 11, p. 113521, 2010.

[15] V. Bragaglia, K. Holldack, J. Boschker, F. Arciprete, E. Zallo, T. Flissikowski and R. Calarco, "Far-Infrared and Raman Spectroscopy Investigation of Phonon Modes in Amorphous and Crystalline Epitaxial GeTe-Sb2Te3 Alloys," Scientific Reports, vol. 6,
p. 28560, 2016.

[16] S. Kozyukhin, M. Veres, H. Nguyen and A. Ingram, "Structural Changes in Doped Ge2Sb2Te5 Thin Films Studied by Raman Spectroscopy," *Physics Procedia*, vol. 44, no. 2, pp. 82-90, 2013.

[17] L. Prazakova, E. Nolot, E. Martinez, F. Fillot, D. Rouchon, N. Rochat, M. Bernard, C. Sabbione, D. Morel, N. Bernier, A. Grenier, A.-M. Papon, M.-C. Cyrille and G. Navarro, "Temperature Driven Structural Evolution of Ge-rich GeSbTe Alloys and Role of N-doping," in *E2021 - European Phase-Change and Ovonic Symposium*, Oxford, United Kingdom, 2021.

[18] R. Golovchak, Y. Choi, S. Kozyukhin, Y. Chigirinsky, A. Kovalskiy, P. Xiong-Skiba, J. Trimble, R. Pafchek and H. Jain, "Oxygen incorporation into GST phase-change memory matrix," *Applied Surface Science*, vol. 332, pp. 533-541, 2015.

[19] H. Tompkins and E. Irene, Handbook of Ellipsometry, William Andrew Publishing, 2005.

[20] G. Jellison and F. Modine, "Parameterization of the optical functions of amorphous materials in the interband region," *Applied Physics Letters*, vol. 69, pp. 371-373, 1996.

[21] T. Levert, A. Zakhtser, J. Duval, C. Raguenez, S. Verdier, D. Le Cunff, J.-H. Tortai and B. Pelissier, "Hybridization of ellipsometry and XPS energy loss: Robust band gap and broadband optical constants determination of SiGe, HfON and MoOx thin films," *Microelectronic Engineering*, vol. 283, p. 112117, 2024.

[22] J. Resende, D. Fuard, D. Le Cunff, J.-H. Tortai and B. Pelissier, "Hybridization of Ellipsometry and Energy Loss Spectra from XPS for bandgap and optical constants determination in SiON thin films," *Materials Chemistry and Physics*, vol. 259, no. 1, 2021.

[23] Němec, P., J. Přikryl, V. Nazabal and M. Frumar, "Optical characteristics of pulsed laser deposited Ge-Sb-Te thin films studied by spectroscopic ellipsometry," *Journal of Applied Physics*, vol. 109, p. 073520, 2011.

[24] J. Orava, T. Wágner, J. Šik, J. Přikryl, M. Frumar and L. Beneš, "Optical properties and phase change transition in Ge2Sb2Te5 flash evaporated thin films studied by temperature dependent spectroscopic ellipsometry," *Journal of Applied Physics*, vol. 104, p. 043523, 2008.
A. SUPPLEMENTARY INFORMATION: ELLIPSOMETRY SPECTRA FITTING
METHOD

The ellipsometry is a characterization technique that gives access to thickness or optical
properties of a material/layer (optical indices, bandgap). The desired properties are iteratively
fitted to make a stack model match with the measurements. As explained in the experimental
section, the amplitude ratio ψ and the phase difference Δ are measured and used to define ls
and lc, presently redefined as Eq. (A. 1) and Eq. (A. 2) respectively, for reading convenience:

$$
I_s = \sin(2\Psi) \sin(\Delta) \qquad \text{Eq. (A. 1)}
$$

$$
I_c = \sin(2\Psi) \cos(\Delta) \qquad \text{Eq. (A. 2)}
$$

A stack model is built to generate simulated *ls* and *lc* and to compare them to the measured *ls* and *lc*. In the ellipsometry theory, *ψ* and *Δ* are linked to the optical indices *n* and *k*. If the measured sample is a stack of materials, every layer of this stack must be modelled. All equations used for building *ψ* and *Δ* from optical indices *n* and *k* are extracted from the Handbook of ellipsometry [19]. The polarization state can be divided into *s* and *p* components (perpendicular and parallel to the plane of incidence respectively). In the general case, the cross polarization is considered, leading to four *pp*, *ss*, *sp* and *ps* components. In our case, isotropic layers are considered and the *sp* and *ps* components are equal to 0. The complex reflectance ratio is then defined as Eq. (A. 3):

$$
\rho = r_{pp}/r_{ss} = \tan(\Psi)\exp(i\Delta) \qquad \text{Eq. (A. 3)}
$$

The complex Fresnel reflection coefficients $r_{pp}$ and $r_{ss}$ are defined as Eq. (A. 4) and Eq. (A. 5)
respectively:

$$
r_{pp} = M_{21,pp} / M_{11,pp} \qquad \text{Eq. (A. 4)}
$$

$$
r_{ss} = M_{21,ss} / M_{11,ss} \qquad \text{Eq. (A. 5)}
$$

$M_{11,pp}$, $M_{21,pp}$, $M_{11,ss}$ and $M_{21,ss}$ are part of multiplication matrices $M_{pp}$ and $M_{ss}$ which are characteristic of the layer stack. Let us consider a n-layered stack over a substrate, each layer having its own complex optical indices $\tilde{n} = n-ik$, its own thickness $d$ and being irradiated by the light (of wavelength $\lambda$) at an incident angle $\phi$ (given by the Snell's law). The multiplication matrices $M_{pp}$ and $M_{ss}$ are defined as Eq. (A. 6) and Eq. (A. 7) respectively:

$$
M_{pp} = \chi_{0,pp} \left( \prod_{j=1}^{n} P_{j,pp} \right) \chi_{sub,pp} \qquad \text{Eq. (A. 6)}
$$

$$
M_{ss} = \chi_{0,ss} \left( \prod_{j=1}^{n} P_{j,ss} \right) \chi_{sub,ss} \qquad \text{Eq. (A. 7)}
$$

Where $\chi_0$ and $\chi_{sub}$ are matrices related to the top layer and to the substrate respectively - Eq. (A. 8), Eq. (A. 9), Eq. (A. 10) and Eq. (A. 11):

$$
\chi_{0,pp} = \frac{1}{2} \begin{pmatrix} 1 & \cos(\phi)/\tilde{n}_0 \\ -1 & \cos(\phi)/\tilde{n}_0 \end{pmatrix} \qquad \text{Eq. (A. 8)}
$$

$$
\chi_{0,ss} = \frac{1}{2} \begin{pmatrix} 1 & 1/\tilde{n}_0 \cos(\phi) \\ 1 & -1/\tilde{n}_0 \cos(\phi) \end{pmatrix} \qquad \text{Eq. (A. 9)}
$$

$$
\chi_{sub,pp} = \begin{pmatrix} \cos(\phi_{sub})/\tilde{n}_{sub} & 0 \\ 1 & 0 \end{pmatrix} \qquad \text{Eq. (A. 10)}
$$
$$ \chi_{sub,ss} = \begin{pmatrix} 1/\tilde{n}_{sub} \cos(\phi_{sub}) & 0 \\ 1 & 0 \end{pmatrix} \qquad \text{Eq. (A. 11)} $$

And where $P_j$ are the transfer matrices related to the jth layer - Eq. (A. 12) and Eq. (A. 13):

$$ P_{j,pp} = \begin{pmatrix} \cos(b_j) & -i \cos(\phi_j) \sin(b_j)/\tilde{n}_j \\ i \tilde{n}_j \sin(b_j)/\cos(\phi_j) & \cos(b_j) \end{pmatrix} \qquad \text{Eq. (A. 12)} $$

$$ P_{j,ss} = \begin{pmatrix} \cos(b_j) & i \sin(b_j)/\cos(\phi_j) \tilde{n}_j \\ i \tilde{n}_j \cos(\phi_j) \sin(b_j) & \cos(b_j) \end{pmatrix} \qquad \text{Eq. (A. 13)} $$

With $b_j \ll 1$ and defined as Eq. (A. 14):

$$ b_j = 2\pi d_j \tilde{n}_j \cos(\phi_j)/\lambda \qquad \text{Eq. (A. 14)} $$

The optical indices were fitted using a Tauc-Lorentz (TL) model [20], a physically relevant dispersion law that respects the Kramers-Kronig relations. The method was already detailed elsewhere [21, 22]. The TL models the complex dielectric function of a material, defined as Eq. (A. 15):

$$ \tilde{\epsilon} = \epsilon_r + i\epsilon_i \qquad \text{Eq. (A. 15)} $$

The optical indices $n$ and $k$ are then given by the Eq. (A. 16) and Eq. (A. 17) respectively:

$$ n^2 = \frac{1}{2} \left( \epsilon_r + \sqrt{\epsilon_r^2 + \epsilon_i^2} \right) \qquad \text{Eq. (A. 16)} $$

$$ k^2 = \frac{1}{2} \left( -\epsilon_r + \sqrt{\epsilon_r^2 + \epsilon_i^2} \right) \qquad \text{Eq. (A. 17)} $$

In this model, the complex part of the dielectric function is the product of the Tauc density of state and the Lorentz oscillator part - Eq. (A. 18), Eq. (A. 19) and Eq. (A. 20):

$$ \tilde{\epsilon}_{TL} = \epsilon_{r,TL} + i(\epsilon_{i,T} \times \epsilon_{i,L}) \qquad \text{Eq. (A. 18)} $$

$$ \text{with} \quad \begin{cases} \epsilon_{i,T}(E > E_g) = A_T (E - E_g)^2 / E^2 \\ \epsilon_{i,T}(E \le E_g) = 0 \end{cases} \qquad \text{Eq. (A. 19)} $$

$$ \text{and } \epsilon_{i,L}(E) = \frac{A_L \times E_0 \times C \times E}{(E^2 - E_0^2)^2 + C^2 \times E^2} \qquad \text{Eq. (A. 20)} $$

Where $E_g$ is the bandgap of the material, $E$ the photon energy, $A_T$ the Tauc coefficient, $A_L$ the peak strength parameter, $C$ the broadening parameter, and $E_0$ the central energy peak. For a N-oscillator model, the imaginary part of the dielectric function is then defined as Eq. (A. 21):

$$ \left\{ \begin{array}{l} \epsilon_i(E > E_g) = \sum_{i=1}^{n} \frac{1}{E} \times \frac{A_{T,i} \times A_{L,i} \times E_i \times C_i \times (E - E_g)^2}{(E^2 - E_i^2)^2 + C_i^2 \times E^2} \qquad \text{Eq. (A. 21)} \\ \epsilon_i(E \le E_g) = 0 \end{array} \right. $$

The real part of the dielectric function is obtained from the Kramers-Kronig relation, which links the real and the imaginary part - Eq. (A. 22):
$$ \epsilon_r(E) = \epsilon_r(\infty) + \frac{2}{\pi} \times P \times \int_{E_g}^{\infty} \frac{\xi \times \epsilon_i(\xi)}{\xi^2 - E^2} d\xi \qquad \text{Eq. (A. 22)} $$

With P the Cauchy principal value containing the residues of the integral and $\epsilon_r(\infty)$ the real permittivity at low frequency. This integral was derived by Jellison and Modine [20] and led to the following analytical expression in the case of N oscillators - Eq. (A. 23), Eq. (A. 24), Eq. (A. 25), Eq. (A. 26), Eq. (A. 27) and Eq. (A. 28):

$$ \epsilon_r(E) = \epsilon_r(\infty) + \sum_{i=1}^{n} \left\{ \frac{A_i \times C_i \times a_{ln,i}}{2 \times \pi \times \xi_i^4 \times \alpha_i \times E_i} \ln \left( \frac{E_i^2 + E_g^2 + \alpha_i E_g}{E_i^2 + E_g^2 - \alpha_i E_g} \right) \right. \\ \left. - \frac{A_i}{\pi} \times \frac{a_{atan,i}}{\xi_i^4 \times E_i} \left[ \pi - \arctan \left( \frac{\alpha_i + 2E_g}{C_i} \right) + \arctan \left( \frac{\alpha_i - 2E_g}{C_i} \right) \right] \right. \\ \left. + \frac{4A_i \times E_i \times E_g \times (E^2 - \gamma_i^2)}{\pi \times \xi_i^4 \times \alpha_i} \left[ \arctan \left( \frac{\alpha_i + 2E_g}{C_i} \right) + \arctan \left( \frac{\alpha_i - 2E_g}{C_i} \right) \right] \right. \\ \left. - \frac{A_i \times E_i \times C_i \times (E^2 + E_g^2)}{\pi \times \xi_i^4 \times E} \ln \left( \frac{|E - E_g|}{E + E_g} \right) \right. \\ \left. + \frac{2A_i \times E_i \times C_i \times E_g}{\pi \times \xi_i^4} \ln \left[ \frac{|E - E_g|(E + E_g)}{\sqrt{(E_i^2 - E_g^2)^2 + E_g^2 C_i^2}} \right] \right\} \qquad \text{Eq. (A. 23)} $$

Where:

$$ a_{ln,i} = (E_g^2 - E_i^2) \times E^2 + E_g^2 \times C_i^2 - E_i^2 \times (E_i^2 + 3E_g^2) \qquad \text{Eq. (A. 24)} $$

$$ a_{atan,i} = (E^2 - E_i^2) \times (E_i^2 + E_g^2) + E_g^2 \times C_i^2 \qquad \text{Eq. (A. 25)} $$

$$ \alpha_i = \sqrt{4E_i^2 - C_i^2} \qquad \text{Eq. (A. 26)} $$

$$ \gamma_i = \sqrt{E_i^2 - C_i^2 / 2} \qquad \text{Eq. (A. 27)} $$

$$ \xi_i^4 = (E^2 - \gamma_i^2)^2 + \alpha_i^2 \times C_i^2 / 4 \qquad \text{Eq. (A. 28)} $$

Depending on the number N of oscillators, the TL model uses 3N+2 fitted parameters, which are the material bandgap $E_g$, its real permittivity $\epsilon_r(\infty)$ and 3 parameters per oscillator ($A_i$, $E_i$ and $C_i$). These parameters then allow to generate the material permittivity and the optical indices. Adding the layer thickness to the fitted parameters, these indices are used to generate $I_s$ and $I_c$ through the stack model light propagation calculations. A comparison between measured $I_s$, $I_c$ and generated $I_s$, $I_c$ allows an iterative model adjustment.

In our case, optical indices, and thicknesses of the Si substrate, the SiO$_2$ layer, and the TiN top layer are known on the whole energy range and were considered as fixed parameters. The temperature evolution of the optical properties of Si was considered using different optical indices (at room temperature, 50°C, 100°C, 150°C, 200°C, 250°C, 300°C, 350°C, and 400°C – taken from the RefractiveIndex.info database). The GST thickness is also known but was allowed to be adjusted. Only the optical indices and the thickness of the GST layer are therefore fit in our model using TL oscillators.

All datapoints ellipsometry spectra were fitted using DeltaPsi software with an appropriate stack model. Three Tauc-Lorentz oscillators were used. Examples of $I_s$ and $I_c$ models, compared with measured data, are shown at 25°C (before transition) and at 410°C (after
transition) in Fig. A. 1. For this study, the fitting was deemed satisfactory, and no other
models were tested.

<figure><img src="image_21.png" /><figcaption>Fig. A. 1: Ellipsometry Is and Ic spectra measured (in blue) and modelled (in red) at (a)(b) 25°C and (c)(d) 410°C.</figcaption></figure>