Thermal expansion and compressibility of single-crystal silicon between 285 K and 320 K

The absolute length of a single-crystal silicon gauge block was measured by interferometry in the temperature range between 285 K and 320 K and at different air pressures from atmospheric conditions down to 10−5 hPa. From the obtained dataset, the coefficient of thermal expansion (CTE) was determined as well as the compressibility—or the bulk modulus—of single-crystal silicon in consideration of a systematic correction of the refractometer used. As the choice of the underlying model for the evaluation is not unambiguous, a Bayesian model averaging approach was applied to take into account possible model errors in the uncertainty evaluation. The result of the CTE is not only in agreement with the recommended reference data of CODATA, but provides a standard uncertainty of less than 1 × 10−9 K−1, which is less than half the uncertainty stated so far in the relevant temperature range.


Introduction
Justified by the need for a reference material for thermal expansion measurements of high precision (e.g. [1,2] and others), there have been a variety of thermal expansion measurements on silicon over a wide temperature range in the past. Due to its diamond-like crystallographic structure, singlecrystal silicon provides isotropy with regard to thermal expansion. Moreover, the industry-driven availability of high-purity material makes silicon an ideal candidate-in contrast to technical-application specimens as described in [3]. Also the recent revision of the 'Mise en pratique for the definition of the metre in the SI' [4] refers to the lattice spacing of silicon as a basis for secondary methods of realising the metre on the Original Content from this work may be used under the terms of the Creative Commons Attribution 4.0 licence. Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. nanometre scale. In this context, knowledge of the CTE can be used in the secondary realisation of the length as described in the updated Mise en pratique. Overview reports on thermophysical properties, e.g. [5], usually refer to [6], which presents a compilation of several measurement results from different laboratories covering the temperature range from 90 K up to 850 K and has acted as a reference source of thermal expansion data of silicon ever since its publication 1 . Some of the involved data were taken from [8][9][10], and a later extension of the temperature range up to 1300 K [11] was taken into account afterwards [12]. In the scope of a programme to establish a thermal expansion standard, part of the upper temperature range above 600 K has been reinvestigated after experimental refinements [13]. The results of thermal expansion measurements between 7 K and 293 K have been presented in [14] in comparison with data from the National Metrology Institute of Japan (NMIJ) and the Jet Propulsion Laboratory (JPL) [15], indicating a systematic deviation from the CODATA reference data. Data of the highest precision in quite a limited temperature range close to 20 • C have been reported earlier in [16], accompanied by the determination of the material's compressibility. In contrast to the dilatometric measurements from other works, the results in [14,16] are derived from absolute length measurements. The work is followed up by the present study of thermal expansion extending the temperature range to 320 K with a reduced measurement uncertainty. Moreover, the simultaneous determination of the compressibility of silicon resolves the discrepancy of the previous results stated in [16] and match with data from the literature [17] (referring to [18,19]).
The results in this work originate from measurements carried out on a single-crystal silicon sample, the dimensions of which are 197 mm × 35 mm × 9 mm (see figure 1). The original high-purity crystal ingot-called #41969/05was grown by Wacker Siltronic GmbH, Germany, using the float-zone technique in an argon atmosphere without any doping [20]. The gauge block sample was cut from the dislocation-free '0-zone' region along the ⟨100⟩ direction between the axial crystal positions at 75 cm and 103 cm (see figure 2 in [20]). The two small end faces with a 35 mm × 9 mm cross section were lapped to achieve optical quality and parallel orientation perpendicular to the long axis. The sample used in [16] was cut from the same crystal ingot out of the same '0-zone' region and axial position and has the same dimensions.
The measurement data were analysed using a novel technique recently presented in [21]. Being a quantity that is computed via a derivative, the coefficient of thermal expansion turns out to be quite sensitive to the concrete model that is chosen to fit the thermal expansion measurements. In [21], a new approach, based on Bayesian model averaging, was therefore proposed that allows various models to be treated and compared in one framework. By taking model errors into account, a more coherent inference can be conducted. Furthermore, via the computation of model probabilities, a more direct comparison can be performed.

Experimental setup
The details of the Twyman-Green-type imaging interferometer used (figure 2) have already been described comprehensively in [22]. For the length measurements, the light of three stabilised lasers (wavelengths: 532 nm, 633 nm and 780 nm) is successively fed into the interferometer beam path. The gaugeblock-shaped sample as well as the refractometer cell, which is used for the interferometric in-situ determination of the air refractive index [23][24][25], are placed on the sample platform (similar to [16]). The mirror behind the refractometer cell reflects the transmitted light towards the camera so that the cell is passed twice. The full field of view is imaged onto the camera, by which a sequence of ten interferograms is recorded, enabling the averaged five-frame phase-shifting approach [22,26] applied for the evaluation of the interference phase. In the focus of the exit collimator, an aperture stop  blocks parasitic reflections from surfaces in the beam path (e.g. chamber windows) leaving the interferometer. A wedged glass plate in the exit beam path compensates for wavelengthdependent images shifts so that the beam locations are almost identical for all three wavelengths applied in the measurements [27].
In order to achieve constant thermal conditions, the interferometer is set up in a pressure-tight/vacuum chamber enveloped by a tubing system with water from a thermostat flowing through it. To monitor the temperature of the sample, three thermocouple sensors [28] are attached to its surface at different locations. The reference junction of each sensor couple is fixed to a copper block, the absolute temperature of which is measured by a Pt-25 resistance thermometer calibrated at the fixed points of water and gallium according to the ITS-90 2 [29]. This operating principle requires the thermocouple sensors to only measure small temperature differences, i.e. a few tens of millikelvins at maximum. Six additional thermocouple sensors are positioned around the sample and the refractometer cell to keep track of the temperature homogeneity of the surrounding gas medium when the chamber is not evacuated. The pressure of the gas, i.e. air, is measured by traceably calibrated sensors of the MENSOR CPT6100 or CERAVAC CTR 91 (1 Torr) type, depending on the actual pressure ranges. The air is kept dry by means of a cooling trap consisting of a tube which is dipped into a Dewar vessel filled with liquid nitrogen. The vacuum state, i.e. less than 10 −4 hPa, is monitored by a sensor of the type IONIVAC IE 20.

Measurement strategy
The investigation of the thermal expansion and the compressibility based on absolute length measurements can be realised by a convenient temporal sequence of temperature and pressure steps. A specific order is not necessary. Therefore, the parameter space of temperature T and pressure p has been sampled adequately within the intervals of 285 K < T < 320 K and 10 −5 hPa < p < 1030 hPa, respectively. The histogram in figure 3 shows the corresponding distribution of measurements which have arisen from experimental circumstances or practical constraints. For instance, keeping the chamber at temperatures above 300 K for something more than a day leads to condensation liquid appearing on the inner side of the windows of the chamber. Therefore, the measurement point density is kept smaller in this temperature region. The data acquisition took place over a period of approximately three months. After each change of the temperature or the pressure, a waiting period of at least half a day was necessary to achieve thermal equilibrium again. Length measurements were only carried out when the temperature drift was close to zero. Near room temperature, i.e. 20 • C, also the measured temperature differences between the different sensor locations were negligibly small. Away from this reference temperature, however, spatial temperature differences of up to 10 mK may occur, which have to be taken into account in the temperature evaluation and the corresponding uncertainty analysis.

Data evaluation
The principle of the interferogram analysis, i.e. the evaluation chain from the camera images to the resulting length values, has been described in detail in [22]. A typical interferogram is shown in figure 4. The interference phase is calculated from a sequence of phase-shifted interferograms so that the phase difference ϕ between regions of interest on the sample's front surface in relation to the attached reference plate can be determined. For each wavelength of the light used, the length of the sample is then calculated by in terms of the measured integral and fractional orders of interference N and ϕ/2 π, respectively, the vacuum wavelength λ and the refractive index n of the surrounding (gas) medium. An estimate for the value of N is obtained from a separate mechanical measurement with an uncertainty in the order of 1 µm in combination with an extrapolation to the actual temperature and pressure conditions during the measurement. The method of exact fractions is then applied to the measurement results from the different laser wavelengths to find the correct order of interference N [23]. If the pressure is below 1 hPa, the refractive index n is evaluated by the modified Edlén formula [30] involving the measured gas parameters, which yields sufficient accuracy in this pressure range [31]. At higher pressures, the phase difference between regions of interest corresponding to the beam path through the evacuated refractometer cell and the path along the cell is determined analogously to the described length evaluation. This is done in order to measure the actual refractive index in situ. In this case, the result from Edlén's formula is only considered as the required prior knowledge. As the windows of the refractometer cell implicate non-uniform geometric conditions depending on the location in the field of view, a correction of the accompanying optical path length difference has to be applied. This is realised by determining the optical path length difference with the evacuated refractometer cell being in the evacuated vacuum chamber. The resulting value is considered as a corrective offset for the measurements in air. As these vacuum measurements were carried out repeatedly during the entire measurement period, the correction is monitored continuously and no systematic dependence on the actual temperature is found. Moreover, measurements at pressures above 50 hPa require an additional correction due to stress-induced optical path length changes in the windows of the refractometer cell [24,25,32]. This is because the optical paths through and along the glass tube are affected differently. Since-unlike the pressure-induced stress-temperature-dependent changes of the refractive index of the cell windows are sufficiently homogeneous across the windows, there is no need for a temperature-dependent correction. At this point, the subsequent processing of the collected length data at different temperatures and pressures combined with a coherent evaluation of the corresponding measurement uncertainty is in the foreground. The approach used in the following sections is founded on the absolute length values and their individual uncertainty budgets, the latter of which are important due to significant pairwise correlations. Two example uncertainty budgets representing different experimental conditions are listed in tables 1 and 2. These examples represent the ideal case of 20 • C and vacuum, which provide very stable experimental conditions, and a worse case at elevated temperature and atmospheric pressure causing larger temperature gradients. Some of the involved uncertainty contributions are negligibly small (e.g. the CO 2 content of the air) and, therefore, omitted for clarity, here. Other quantities contribute significant correlations to the evaluation of the CTE and the compressibility. These correlations of different length measurements are, for instance, due to the common calibration of the temperature sensors or commonly applied corrections. It is to be noted that the uncertainty contributions of the corrections considering the influence of surface roughness, optical phase shift and wringing contact are intentionally set to zero. Even though, in fact, they contribute a dominating portion (i.e. several nanometres) to the total uncertainty of each individual absolute length value, the uncertainty variation within the parameter space of temperature and pressure become clearer without them. Zeroing these contributions is justified, because the uncertainty of the results of the thermal expansion and the compressibility are not affected. The reason is, that both quantities are differentially deduced from the length data, which are correlated with regard to the mentioned corrections.

Model equation for L(T, p)
In order to consider the two-dimensional parameter space of the length data L(T, p), an appropriate model function has to be applied in the evaluation via a least squares regression. Therefore, we start with the total differential of L(T, p) which is given by Considering a homogeneous material like pure single-crystal silicon, the coefficient of linear thermal expansion at constant pressure p is and the isothermal compressibility κ as the reciprocal of the isothermal bulk modulus K is defined as With these relations, equation (2) can be rearranged as While α l is assumed to change little with pressure 3 , the temperature influence on κ is expected to be negligible in the covered temperature interval. Hence, a linear dependence of α l on the pressure is considered by where α(T) is the coefficient of linear thermal expansion at p = 0 andα is a constant parameter. Actually, considering the temperature dependence of the compressibility [34], the same approach could also be applied to κ, but the principle of which can be approximated by a Taylor expansion of the exponential function: Finally, there remains the choice of a suitable expression for α(T).

Model equation for α(T)
As for the selection of an appropriate model for the CTE α(T), there are two approaches which are considered in the following. In [14], the evaluation of the CTE of single-crystal silicon covering a temperature range from 7 K to 293 K has been reported. Because of the nonlinear behaviour of the CTE of single-crystal silicon in the respective large temperature interval (particularly at cryogenic temperatures), a simple polynomial fitting was not adequate. Alternatively, a piecewise evaluation of smaller temperature intervals by polynomials as used in [22] was at the very least not so pleasant due to the incoherent treatment of the data. Moreover, the required polynomial degrees may vary and an extrapolation beyond the limits of the temperature range covered by the actual measurements is only reasonable based on a physically motivated model. Therefore, the appropriate approach from [35]already introduced in [36]-was applied, in which the CTE is modelled by an empirical nonlinear function based on socalled Einstein terms. However, in the current case of a limited temperature range between 285 K and 320 K, both model approaches appear to be applicable, which raises the question of how the uncertainty is affected by the decision for or against one particular model. Therefore, the method of Bayesian model averaging has been adapted and applied to the evaluation of the CTE and its uncertainty based on measured absolute length data [21]. This approach (in contrast to [14] or [37]) does not require a commitment to a specific model but allows several competing models to be incorporated and compared in one framework. Starting from a priori probabilities for each considered model, which can for instance all be chosen equally, these probabilities are updated in the light of the data using Bayes' theorem. As an outcome, one obtains a posteriori probabilities for each model, which allow a direct comparison. Moreover, deduced estimates and uncertainties in Bayesian model averaging avoid inconsistencies between models and take model errors into account. This is of particular importance in the context of CTE determination as the quantity of interest is often obtained via the derivation of fitted curves-a badly behaved operation that might lead to distinct discrepancies between two models yielding practically indistinguishable fits on the measurement data. (10) with the parameters ζ 1 , ζ 2 and ζ 3 as well as by the nonlinear Einstein-term expression

Accordingly, α(T) is represented by the polynomial expression
with the parameters ξ and θ in the evaluation of the results in the following sections. Here, the explicit proportionality of α and κ (cf. [36]) is omitted for simplicity and implicitly considered by the parameter ξ. The respective number of polynomial degrees and Einstein terms is chosen so that the model functions can sufficiently represent the experimental data without over-or underfitting them.

Results and uncertainty
With the two selected variants of α(T), the model prototype L(T, p) from equation (9) finally becomes and respectively. The reference temperature is set to T ref = 293.15 K and the quantities L 20 ≡ L(T ref , 0 hPa) and and ζ ′ 3 = L 20 · ζ3 3 are treated as fit parameters. The parameters of both models, (12) and (13), are estimated via a weighted least squares regression of the measured data to derive the thermal expansion and the compressibility in a second step. The resulting fit parameters are with the corresponding covariance matrix  Please note that the physical units in the entries of the covariance matrices are omitted for simplicity only. Moreover, the order of the entries is in accordance with the order of the listed results of the fit parameters. It must be emphasised that weighting of the optimisation is done based on the full covariance matrix of the entire dataset-in contrast to simply using the reciprocal squared standard uncertainties as can be found in most applications. Doing so coherently takes into account the significant mutual correlation between the data points and yields the correct output uncertainties of the fit parameters of both models.
As an example, the result of the Einstein-term model is shown in figure 5 for the two cases of vacuum and atmospheric pressure. Since the pressure dependence of the length is significantly smaller compared to the effect of temperature changes on the chosen scale, the two plotted lines are hard to distinguish. The corresponding standard uncertainty is shown in figure 6 together with the related residuals. The present results have been evaluated based on the measured temperatures according to the ITS-90. A correction with regard to thermodynamic temperatures [38] has not been applied.

Thermal expansion
The expression for the thermal expansion can be derived following equation (3) from both the above-mentioned fits. As outlined in section 4.2, the method of Bayesian model averaging, which has been described in [21], is applied in this evaluation to take into account competing model variants. Therefore, the two selected models are averaged considering the corresponding a posteriori model probabilities: where L P and L E are evaluated with the parameters specified in the previous section. For the computation of posterior probabilities, we followed the same strategy as in [21]. Using data from [16] as prior knowledge for scaling, we found posterior model probabilities P E = 65.9 % and P P = 34.1 % for  (13) at 0 Pa and 1000 hPa, respectively. The blue dots represent the residuals of each measured value considered in the weighted least squares regression. Please note that, although plotted versus temperature only, the data points represent the entire pressure range from vacuum to atmospheric pressure. the data presented in this work. These are close to the values found in [21]. The resulting course of the thermal expansion depending on the temperature based on the averaged model is shown in figure 7 and the difference of the models α P (T, p) and α E (T, p) is negligibly small compared to the corresponding uncertainties. These standard uncertainties of α P (T, p), α E (T, p) andα(T, p) are plotted in figure 8. While in the middle section of the temperature interval, the uncertainties of α P (T, p) and α E (T, p) are almost equal, they differ clearly at the interval edges. This is taken into account by the Bayesian model averaging in the resulting uncertainty ofα(T). For convenience, some values are listed in table 3.
In [6,7], the compilation of several sources has been prepared which-to our knowledge-presents the lowest uncertainties (in our relevant temperature range mainly based on the data from [9]). Taking this as a reference for the results from the current work enables a comparison which is plotted in figure 9. All data agree well within their corresponding uncertainties over the whole temperature interval, while the current data provide a standard uncertainty, which is smaller by a factor of at least about five compared to the reference data.  It must be noted that the current result is related to the current temperature scale, ITS-90, while the reference data refer to the former IPTS-68 4 . A correction would imply a shift by less than 10 mK without a significant effect on the evaluated CTE in the given temperature interval.
The possible pressure dependence of the thermal expansionα is taken into consideration in the model.  . The light grey area shows the standard uncertainty of the currently most precisely known reference data of the coefficient of thermal expansion of silicon in the given temperature range according to [6]. The dark grey area represents the difference of the results from the present work according to equation (14) with regard to these reference data including the corresponding standard uncertainty. this work, this example does not reveal a significant pressure dependence with regard to the relevant measurement uncertainty. As a consequence, in the application of interferometry on gauge blocks, it is also permissible to apply the coefficient of thermal expansion determined in vacuum for measurements in atmospheric conditions without any further correction.

Compressibility
In [16], the compressibility κ was derived from measured length data at different pressures between vacuum and atmospheric pressure while keeping the temperature very close to 20 • C. In the present case, the length measurements are distributed widely in the parameter space of temperature and pressure, and the model equations (12) or (13) take both parameters inherently into account. Proceeding analogously to section 5.1 one gets the model averageκ of the compressibility of singlecrystal silicon. The numerical result at 20 • C and atmospheric pressure isκ = 1.018(5) × 10 −11 Pa −1 (15) or, consequently, noted as the bulk modulus Obviously this value of κ markedly differs from the previous result published in [16] and, moreover, is in close agreement with earlier low-uncertainty results from the literature quoted there [18,19]. The suspected explanation that material differences cause the deviation may not be the real reason. Instead, from today's point of view, there might be another main cause.
In the evaluation in [16], the required correction of stressinduced optical path length changes in the refractometer windows [24,25] was unknown and could therefore not be taken into account. This interpretation is supported by the fact that the result of this work would also be close to the old (wrong) value, if the correction of the window stress were omitted. A temperature-dependent change of the bulk modulus (cf. [34]) could not be seen from the present results (figure 10) with sufficient significance.

Summary
The absolute length of a silicon gauge block was investigated in a series of measurements at different temperatures and air pressures. From the collected dataset, the coefficient of thermal expansion α as well as the compressibility κ of silicon were derived. The evaluation is founded on a Bayesian model averaging approach to take different models and their possible errors additionally into account. While the uncertainties of the results reported previously [14] were estimated rather conservatively, the careful consideration of correlations lead to a reduction of the uncertainty presented in this work. For the CTE, standard uncertainties of less than 1 × 10 −9 K −1 were achieved. These uncertainties are up to one order of magnitude smaller compared to the results reported previously [14]. In view of the data from [14] and this work, an update of the reference data recommended so far should be taken into consideration. The result of the compressibility is in agreement with data obtained by different measurement techniques described in the literature [17][18][19]. Moreover, it enables the resolution of the discrepancy in values reported in [16]. Ultimately, considering the achieved uncertainties, no significant dependence of the CTE on the actual pressure can be determined in the covered parameter space of temperature and pressure. As a consequence, it is permissible to apply the CTE value measured in vacuum for thermal expansion corrections of length measurements carried out in atmospheric pressure.