Determination of thermal diffusivity of thermoelectric materials using a micro four-point probe method

To develop materials with higher thermoelectric efficiency, a comprehensive characterization of material properties on the relevant spatial scales is pertinent. In this study, we develop a new method based on a micro four-point probe (M4PP) to determine the thermal diffusivity of a bulk material using the phase delay of the second harmonic voltage. The method is demonstrated on two relevant thermoelectric materials i.e. skutterudite and bismuth telluride. Individual measurements on skutterudite and bismuth telluride are characterized by a precision better than 9% and 12%, with mean thermal diffusivities of 1.87 ± 0.17 mm2/s and 1.14 ± 0.13 mm2/s for 100 measurements, respectively. These values are found to be in good agreement with their independent estimates. M4PP electrical characterization of bismuth telluride shows signs of percolation network behavior, thus we believe that the new M4PP methodology applies even to materials whose electrical resistivity is nonuniform and erratic.


Introduction
Thermoelectric materials can convert waste heat into electricity, even when the heat sources are not easily accessible by other technologies (e.g. low-temperature heat sources), and the conversion proceeds without harmful emissions. Higher efficiency thermoelectric materials may contribute significantly towards a more energy-sustainable future [1]. However, current state-of-the-art materials are limited in terms of efficiency, hampering the spread of this technology. The efficiency of thermoelectric materials is defined by the figure of merit, ZT = σα s 2 T/κ, where σ is the electrical conductivity, α s is the Seebeck coefficient, κ is the thermal conductivity, and T is the absolute temperature. Hence, even though the characterization of all these properties is a laborious task, its importance for the development of new materials cannot be underestimated. The thermal conductivity often appears particularly challenging to determine, since it typically involves subjecting a sample to a thermal gradient and measuring the heat flux, which introduces significant measurement uncertainties [2]. One way to avoid such uncertainties is by measuring thermal diffusivity D = κ/ρc P instead of measuring the thermal conductivity κ, and obtaining the latter through additional characterization of density ρ and specific heat c P [3].
Commonly, thermal diffusivity is measured with the hot disk [4] or the laser flash technique [3], but these methods are limited to particular geometries and are practically inapplicable to small sample dimensions. In these cases, the two commonly used techniques are time-domain or frequency-domain thermoreflectance [5] and scanning thermal microscopy [6]. However, these techniques heavily depend on thermal models used in data analysis [7]. In addition, thermoreflectance techniques require sample preparation, which usually involves the deposition of an additional metallic layer. Another technique recently presented for surface mapping of thermal conductivity on materials is photo-thermal imaging [8]. This infrared photo-activated thermography technique has a resolution down to ~10 μm, and can measure thermal conductivity in the range 0.1-100 W/mK [8].
An emerging alternative for thermal characterization on the microscale is the micro four-point probe (M4PP), which consists of four microscopic pins (two for applying and extracting current, and the other two for recording the voltage difference). Given the small dimensions of the probes (in the μm range), M4PP is often used for spatial characterization of material properties, and results are used for process control and optimization. The mainstream uses of M4PP include the measurement of sheet resistance [9,10], electron mobility [11,12], carrier density [11,12], and magnetoresistance [13,14]. Recently, however, novel M4PP techniques have been successfully developed and applied to address the thermal properties of materials, including the temperature coefficient of resistance (TCR) in thin film metals [15], and the ratio of the Seebeck coefficient to the thermal conductivity α s /κ in bulk semiconductors [16]. In these recent developments, M4PP measurements were done at a single low frequency, and data analysis was restricted to only the amplitude of selected higher harmonics, namely the second harmonic to constrain α s /κ [16], and the third harmonic to extract TCR [15].
In this study, we show that the variation of the four-point voltage with frequency provides an even richer information about the material properties of interest. In particular, we utilize the phase delay of the second harmonic voltage to estimate the thermal diffusivity. The technique is demonstrated on two different bulk thermoelectric materials, namely (i) a skutterudite showing a predominantly 3D electrical transport, and (ii) a commercial bismuth telluride that electrically behaves like a percolation network. The presented measurement technique falls inside the realm of 2ω methods, however heat is generated at the electrode contacts within the sample (an internal heater), and thus the technique does not require the deposition of additional materials for external heat generation.

Theoretical model
During an M4PP measurement, an AC current I=I 0 sin(2πft) at a frequency f is passed between two electrodes through a sample (see Fig. 1a) [16]. The Peltier and Joule effects produce a temperature difference between the voltage pins ( Fig. 1b), which in turn gives rise to a voltage difference via the Seebeck effect. The direct contribution of the Peltier effect to the voltage reading is recorded at the same frequency as the applied signal since this effect is linearly proportional to the current. However, the direct contribution of Joule heating to the measured voltage is recorded as a DC term and at the second harmonic frequency, 2f, since the relation between current and heat generated is quadratic. Interestingly, the second harmonic voltage is not in phase with the power introduced in the sample (see Fig. 1c); this phase shift is related to the thermal diffusivity of the material.
The second harmonic voltage difference uv pq ΔV 2ω between a generic pin p and a pin q due to a current injected in pin u and extracted at pin v can be described as [17], where l k ΔT 2ω is the second harmonic temperature rise at the contact of pin k∈{p, q} due to the heating at pin l∈{u, v}, and α e is the Seebeck coefficient of the electrodes. The temperature l k ΔT 2ω can be calculated as (derived in Appendix I), where r k,l is the distance between pins k and l, t is time, I rms is the root mean squared current, R c,l is the electrical contact resistance of pin l, ω is the angular frequency (ω = 2πf), i is the imaginary unit, and D is the thermal diffusivity. The contact resistances R c,l were determined by regressing the (two-point) load resistance of different configurations as explained in Appendix II. Eq. (2) is only valid when the power is deposited near the contacts, and the heat can be approximated as flowing in the radial direction. For the heat to flow in the radial direction, the distance between electrodes must be much larger than the probe contact radii.
Finally, once the second harmonic voltage as a function of time is obtained ( uv pq ΔV 2ω ), a numerical lock-in amplifier is used to extract the phase shift, φ 2ω . The phase shift does not depend on the parameters α s , α e , I rms , and κ, so it can be used in a fitting procedure to determine D when the geometrical parameters of the probe are known.

Experimental setup
To extract the second harmonic phase at different frequencies, M4PP measurements were performed using a modified CAPRES microRSP®-A300 tool equipped with a digital lock-in amplifier module [18]. All measurements were performed with an equidistant four-point probe with a 30 μm electrode pitch. When using probes with smaller probe pitch, higher frequencies must be measured to ensure a significant phase shift. At each measuring point (engage), the electrode configuration C' shown in Fig. 1a and b (and its spatial reciprocal, configuration C, with voltage and current pins swapped) were measured with I rms = 2 mA at frequencies 3.01, 6.03, 12.06, 24.11, 48.22, 96.44, 192.88, 385.76, and 771.52 Hz. 1 Two samples, namely a commercial n-type bismuth telluride from Interm, and an n-type skutterudite (CoSb 2⋅75 Sn 0⋅05 Te 0.20 ) were subjected to an M4PP line scan with 100 points using a step size of 10 μm. The thermal diffusivity of the bulk bismuth telluride sample was 1.147 mm 2 / s, measured by laser flash. The skutterudite sample was previously measured by the same technique and showed a thermal diffusivity of 2.1 mm 2 /s [19,20]. Both samples were polished to remove surface roughness before performing the measurements.

Results
To determine the thermal diffusivity D, the pin separations (r k,l ), must be known. These can be determined experimentally with high accuracy for uniform materials [14] but are challenging for non-uniform materials. For electrode pitch much larger than the contact position uncertainty (r k,l >=30 μm, and error<500 nm), it is adequate to use the nominal electrode positions, and thus we neglect position errors.
At each probe engagement, C and C ′ configurations were used to extract D. Fig. 2aand b show the sine of the experimental second harmonic voltage phase shift (symbols) and best fit (lines) for the configurations C and C', respectively. The insets in both figures show a schematic view of the configuration used. Fittings of the sine of the second harmonic phase were performed using MATLAB. We obtain a single value of thermal diffusivity D for each probe engaged with D as the only fitting parameter. Fig. 2c shows the best-fit values of D and their uncertainties, for the 100 engages (color dots), alongside their independent reference values (color lines). The obtained mean values of D for all engages were 1.87 ± 0.17 mm 2 /s and 1.14 ± 0.13 mm 2 /s for skutterudite and bismuth telluride, respectively. Both samples showed weak third harmonic signals due to their low TCR.
To investigate the current flow, four-point resistance was measured in 6 electrode configurations, defined as R i for the three main configurations i ∈ {A, B, C} and R i ′ for their reciprocal configurations i ′ ∈ {A',B', As demonstrated in Ref. [21] it is convenient to define flow in certain paths and this ratio may significantly deviate from 1.5. Fig. 3aand b show histograms of the resistance ratio (obtained in the same engages as results reported in Fig. 2c) along with the expected probability distribution resulting from measurements on a uniform sample with normally distributed electrode position errors.

Discussion
The mean value of D for the bismuth telluride (1.14 mm 2 /s) was less than 1% lower than the laser flash measurement, but the skutterudite mean D (1.87 mm 2 /s) was 11% lower than expected (2.1 mm 2 /s). A possible reason for this underestimation is the difference in volume measured with laser flash and M4PP, respectively. In laser flash, the characteristic length scale is the thickness of the sample (in this case 2 mm), while M4PP measures on a scale of the electrode pitch (in this case 30 μm). Thus, for the macroscopic and microscopic techniques to agree, the samples must be uniform. From the M4PP measurements, this does not appear to be the case. Since the samples are non-uniform, it is difficult to evaluate the measurement precision. However, we estimate a worst-case measurement precision for the skutterudite and the bismuth telluride to be 9% and 12%, respectively. This was calculated as the standard deviation of all the D values for each material, and it includes  The contact resistances (R c,l ) are the sources of the second harmonic power deposited on the sample. They were accurately calculated as explained in Appendix II. One could expect a small overestimation or underestimation in all pins due to the assumption of a perfect contact between the nickel of the probes and the nickel of the calibration sample. However, it does not influence the fitted D since only a change in amplitude (not phase) of the measured second harmonic voltage is produced. Fig. 4 shows that an overestimation or underestimation of the contact resistances in all pins does not affect the fitted D (in red), and even if only one pin had a considerable deviation, limited deviation on D is produced (in black).
Both samples are polycrystalline and one might expect therefore that the current flow differs from an ideal radial flow. The RÃ/ RB ratios of skutterudite in Fig. 3a are as expected around 1.5, while those of bismuth telluride in Fig. 3b show values ranging more randomly from 1 to 2, indicating that the electrical current transport inside the bismuth telluride sample is confined and not radial. Also, it is not onedimensional as expected for a percolation network with a single preferred current path, in which case the ratio would be unity [21]. The current confinement may increase the challenge of measuring the electrical resistivity since consecutive measurements can provide resistance values more than one order of magnitude apart (and hence, different values of the ratio RÃ/RB). However, even though the electrical current transport is non-uniform in the case of the bismuth telluride sample, the extracted thermal diffusivity seems to be largely unaffected as compared to the skutterudite sample, cf. in Fig. 2c. This observation is likely to be the consequences of a difference in electron and heat transport. In general, electrical conductivity in solids can change by many orders of magnitude, whereas thermal conductivity only changes a few orders of magnitude [22,23]. Thus, electron transport is more likely to be confined, than phonon heat transport.

Conclusions
The origin of the second harmonic voltage response during a micro four-point probe measurement of a bulk thermoelectric material was explained, and a time-domain frequency-dependent expression of the voltage signal was developed. With the help of a numerical lock-in amplifier, the phase shift of the second harmonic voltage was extracted, and it was demonstrated to depend on the thermal diffusivity. Two thermoelectric samples were measured 100 times in a 1 mm straight line, a skutterudite and a bismuth telluride, which showed a worst-case precision of 9% and 12%, and mean thermal diffusivities of 1.87 ± 0.17 mm 2 /s and 1.14 ± 0.13 mm 2 /s, respectively. Interestingly, the thermal diffusivity obtained for the skutterudite was 11% lower than the previous macroscopic characterization by laser flash, but the bismuth telluride sample was in good agreement (less than 1% lower than the laser flash measurement). Finally, using measurements at different configurations and Monte Carlo simulations, it was justified that the bismuth telluride sample behaved electrically like a percolation network, which implies that even samples that possess an electrical resistivity difficult to be measured by the M4PP method can be thermally characterized. This study adds up to a few recent studies showing M4PP as a potential alternative to the commonly used thermoreflectance, and scanning thermal microscopy for thermal characterization at the microscopic scale. It also represents a step forward for the complete characterization of thermoelectric materials (determination of all the properties that define the figure of merit).

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Data availability
Data will be made available on request.
This appendix presents the derivations to obtain the solution of eq. (2), which is the temperature distribution in a semi-infinite material due to periodic heating on a hemispherical indentation of radius b (see Fig. 5). The periodic heating power area density is assumed to be evenly distributed over the surface of the hemispherical indentation.

∂T ∂t
where T represents temperature and t time. This expression in spherical coordinates and assuming radial symmetry reduces to, To solve the problem of Fig. S1, the following harmonic ansatz is used, and, hence, the temperature profile is separated into a spatial-dependent part U(r) and a time-dependent part oscillating at a single angular frequency ω q . Inserting (5) into eq. (4) yields to, This corresponds to the Fourier transform of the differential equation. The solution for this second-order ordinary differential equation is, where the two constants c 1 and c 2 are determined from the boundary conditions. Requiring that the temperature must remain finite at infinity, T(r,t)→ T 0 for r→∞, implies that c 1 = 0. The second boundary condition, originating from Fourier's law and the even distribution of the heating power, where κ is the thermal conductivity of the semi-infinite material, and P 0 is the deposited harmonic power amplitude. This equation is used to determine c 2 , Then, the spatial-dependent part of the temperature becomes, and the temperature profile is given by, The power amplitude of the second harmonic defined by Joule heating is P 0 =I rms 2 R c , where I rms is the root mean square current passing through contact with resistance R c at the distance b from the contact center. This power is deposited at an angular frequency ω 0 = 2ω, where ω is the angular frequency of the current. The temperature profile is then, Finally, we make use of the limit b→0 since we are only considering point heat sources at the current pins, and the expression becomes,

Appendix II
This appendix shows the procedure used in this study to accurately calculate the contact resistances (R c,l ) of each pin. The contact resistances were obtained by fitting load resistance measurements of different configurations, first, on a calibration sample (a nickel thin film, the probe material), and then, on the measured samples (bismuth telluride and skutterudite).
The load resistance between pins 1 and 2, R load,1,2 , is given by, where R lead,l is the lead resistance (resistance from the electronics until the contacts) of pin l, and R sample,1,2 is the resistance of the sample between pin 1 and pin 2. For the calibration measurements, a perfect contact between the probe and the sample is assumed (R c,l = 0 between nickel and nickel), and the sample resistance is estimated as [24], where R S is the sheet resistance, r 1,2 is the distance between pins 1 and 2, and a contact radius b = 129 nm is assumed to be constant for all pins. Hence, the measurements at different configurations define a system of equations where the only unknowns are the lead resistances and the sheet resistance.
Once the lead resistances are determined, the load resistance measurements on the samples under evaluation can be used to determine the contact resistances by defining a similar system of equations (now, the only unknowns are the contact resistances). Notice that, in this case, the sample resistance term is not present since we are measuring 3D samples where the load resistance is adequately describe by the superposition of spreading resistance from the two current pins.