*Article* **Scanning Thermal Microscopy of Ultrathin Films: Numerical Studies Regarding Cantilever Displacement, Thermal Contact Areas, Heat Fluxes, and Heat Distribution**

**Christoph Metzke 1,2 , Fabian Kühnel <sup>1</sup> , Jonas Weber 1,3,4 and Günther Benstetter 1,\***


**Abstract:** New micro- and nanoscale devices require electrically isolating materials with specific thermal properties. One option to characterize these thermal properties is the atomic force microscopy (AFM)-based scanning thermal microscopy (SThM) technique. It enables qualitative mapping of local thermal conductivities of ultrathin films. To fully understand and correctly interpret the results of practical SThM measurements, it is essential to have detailed knowledge about the heat transfer process between the probe and the sample. However, little can be found in the literature so far. Therefore, this work focuses on theoretical SThM studies of ultrathin films with anisotropic thermal properties such as hexagonal boron nitride (h-BN) and compares the results with a bulk silicon (Si) sample. Energy fluxes from the probe to the sample between 0.6 µW and 126.8 µW are found for different cases with a tip radius of approximately 300 nm. A present thermal interface resistance (TIR) between bulk Si and ultrathin h-BN on top can fully suppress a further heat penetration. The time until heat propagation within the sample is stationary is found to be below 1 µs, which may justify higher tip velocities in practical SThM investigations of up to 20 µms−<sup>1</sup> . It is also demonstrated that there is almost no influence of convection and radiation, whereas a possible TIR between probe and sample must be considered.

**Keywords:** scanning thermal microscopy (SThM); numerical study; finite element analysis (FEA); boron nitride; h-BN; ultrathin films; heat transfer; thermal contact; penetration depth; stationary time

#### **1. Introduction**

Since the thermal properties of thin films vary significantly from those of the corresponding bulk materials [1–4], new promising materials for micro and nanoscale devices require a detailed investigation with advanced techniques. Moreover, due to several atomic and molecular effects, such as grain boundaries, or the transition to ballistic heat transport, the thermal characterization becomes increasingly challenging. Methods for thin-film thermal characterization are also limited by various factors such as film thickness or anisotropic material properties. One possible method is SThM, which is specifically designed to characterize the local thermal properties of thin films. SThM is applied in an AFM system together with additional measurement equipment. It is a method in which a cantilever is in direct physical contact with a sample. The sample is scanned in a special pattern to obtain local thermal properties. SThM thermal images are likely to be influenced by the sample's topography, which has been explained in literature in recent years [5–12]. To ensure a correct interpretation of the recorded measurement results, a deep understanding of heat

**Citation:** Metzke, C.; Kühnel, F.; Weber, J.; Benstetter, G. Scanning Thermal Microscopy of Ultrathin Films: Numerical Studies Regarding Cantilever Displacement, Thermal Contact Areas, Heat Fluxes, and Heat Distribution. *Nanomaterials* **2021**, *11*, 491. https://doi.org/10.3390/ nano11020491

Academic Editor: Orion Ciftja

Received: 20 January 2021 Accepted: 15 February 2021 Published: 16 February 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

transfer during SThM measurements is essential to comprehend the origin and impact of those and further effects. This work aims to provide these insights.

In this work, the goal is to accurately predict the heat transfer process during realistic SThM measurements of ultrathin films such as h-BN and a bulk Si sample. Therefore, theoretical calculations and finite element analysis (FEA) are performed. FEA is a versatile tool to simulate, e.g., heat transfer or mechanical problems that can be described by mathematical equations. Similar SThM measurements with h-BN have been performed in recent works [12]. We tried to create theoretical measurement setups that comply with real scenarios so that the theoretical results may be adopted by other researchers to compare them with practical measurements or to explain certain effects of SThM applied to ultrathin films.

The literature in recent years especially focused on qualitative local thermal properties. Leitgeb, Fladischer et al. investigated 500 nm tungsten films employing SThM based on the 3ω method and compared the results to the time domain thermoreflectance technique. Said tungsten films would have a thermal conductivity between 151.4 Wm−1K <sup>−</sup><sup>1</sup> and 156.0 Wm−1K <sup>−</sup><sup>1</sup> depending on the heat treatment [13,14]. Chen et al. estimated the thermal conductivity of a single SiO<sup>2</sup> nanoparticle at 300 K to 0.95 <sup>±</sup> 0.08 Wm−1<sup>K</sup> <sup>−</sup><sup>1</sup> using SThM. They also found out that the TIR between the probe and the nanoparticle accounts for 70% of the total thermal resistance. Therefore, TIRs would have to be considered in thermal conductivity measurements [15]. Existing TIRs from metal/non-metal and nonmetal/non-metal contacts within a sample were studied by Park et al. using ultrahigh vacuum SThM. They suggest the presented method to be used actively for nanoscale TIR measurements and show the significant contribution of TIRs to the entire heat dissipation at the nanoscale [16]. Chirtoc et al. conclude that the heat management of nanofabricated thermal probes might be optimized by decreasing the TIRs between tip and sample [17]. To obtain local thermal properties using SThM, a tip calibration is necessary. A new calibration method for local temperature measurements is introduced by Nguyen et al., whose active thermal microchips could be used for different SThM probes. In addition, they estimate the TIRs between tip and sample, which would have a great impact on the results [18]. Recent literature demonstrates the necessity of detailed insights into the heat transfer process of SThM measurements, but also the great possibilities of SThM. This work aims to provide values regarding energy fluxes, heat distribution and the influence of possible TIRs between tip and sample. Said values can hardly be found in the literature, especially for h-BN or Si samples, which are interesting in our field of research. We want to support other researchers to fully understand and interpret upcoming effects in their practical SThM investigations.

The measurement and simulation procedure applied in this research work is a sequence of interdependent steps. First, the geometry of a used SThM probe was investigated via SEM (Section 4.1). The information of the SEM images provided the fundamentals to model the probe in SolidWorks (SW). The geometry was then imported into COMSOL Multiphysics (CM). Here, the cantilever displacement under specific forces (Section 4.2) was simulated. After a theoretical calculation of the thermal contact area (TCA) between tip and sample (Section 4.3), we were able to model realistic measurement scenarios using a sample with an ultrathin top surface, which exhibits anisotropic thermal conductivity (Section 4.4). Parametric sweep studies enabled the simulation of different probe temperatures and various anisotropic material properties. Subsequently, these simulations were compared to those with a standard bulk Si sample with isotropic thermal properties (Section 4.5). Both simulation setups deliver particular results concerning the heat distribution and the heat fluxes from tip to sample. Henceforth, the time until the heat transfer process becomes stationary was investigated (Section 4.6). The influence of convection and radiation on the present simulations was also taken into account (Section 4.7), as well as the influence of a possible TIR between tip and sample (Section 4.8). Finally, the applied meshing strategy is verified to demonstrate the reliability of our results (Section 4.9).

#### **2. Materials and Methods**

AFM is a method to characterize surfaces according to topographical, mechanical, and electrical properties on a nanometer scale. One of the key parts is a microfabricated probe with an ultrasharp tip to contact the sample surface. Tip radii are in the small nm range (depending on the desired application), e.g., Bruker VITA-DM-NANOTA-200 with a tip radius of up to 30 nm and a tip height of 3–6 µm for a new tip [19,20]. The sample is scanned with a predefined number of lines and readings per line, which results in a line-by-line image. SThM is a method to qualitatively map local thermal conductivities and temperatures as a subcategory of AFM that requires additional measuring equipment. In standard applications, SThM measurements are performed in contact mode. Here, the tip is in direct physical contact with the surface under investigation. The force between tip and sample is held constant in most cases and can be defined by the user during the measurement. Typically, the force of contact mode measurements is within the range of 10 nN to 100 nN [21]. Moreover, a thermal signal is measured and assigned simultaneously to the corresponding topography. Thereunto, a thermal resistive probe is heated with a specific heating power. The temperature of the probe depends on the heat exchange between tip and sample and, therefore, on the sample's local thermal conductivity. If the local thermal conductivity of the sample is high, more heat can spread into the sample, causing a temperature decrease of the probe. The AFM relies on a feedback algorithm and a Wheatstone bridge to evaluate the local thermal conductivities. As a result, SThM measurements create two images simultaneously, one topography and one thermal image of the same position. Such images can be found in [12]. For more details regarding the SThM measurement process, please refer to [22].

In our field of research, special focus is on the thermal characterization of ultrathin films such as h-BN with thicknesses of approximately 10 nm, which are supposed to have anisotropic thermal properties. Recent works demonstrated the possibility of such SThM measurements [12]. This work continues said research and compares it to "standard" SThM measurements of bulk Si with isotropic thermal properties. We deploy the following SThM probe, soft- and hardware:


#### **3. Theoretical Background**

#### *3.1. Cantilever Displacement*

The displacement of a beam w(l), which is fixed at one end and stressed by a single force F at the other end, can be calculated using Equation (1), where l is the cantilever length and E Young's modulus. The moment of inertia I of a cantilever with a rectangular cross-section can be calculated with I = b · h 3 /12, in which b is the width and h is the height of the cantilever [23].

$$\mathbf{w}(\mathbf{l}) = \frac{1}{\mathbf{3}} \cdot \frac{\mathbf{F} \cdot \mathbf{l}^3}{\mathbf{E} \cdot \mathbf{I}} \text{ in [m]} \tag{1}$$

The bending stress σ<sup>b</sup> can be calculated according to Equation (2) by the division of the bending moment M<sup>b</sup> = F · l and the moment of resistance W = b · h 2 /6 (in case of a rectangular beam) [23].

$$
\sigma\_{\rm b} = \frac{\mathbf{M\_{b}}}{W} \text{ in } \left[ \text{Nm}^{-2} \right] \tag{2}
$$

In Section 4.2, Equations (1) and (2) are used to calculate the cantilever displacement and the bending stress in SThM measurements.

#### *3.2. Hertzian Surface Pressure*

Hertzian surface pressure occurs when two rigid bodies touch each other. In the special case of a sphere with radius R, Poisson's ratio ν<sup>1</sup> and Young's modulus E<sup>1</sup> touching a flat surface with Poisson's ratio ν<sup>2</sup> and Young's modulus E<sup>2</sup> under a force F, the indentation depth d can be calculated using Equation (3) [24].

$$\mathbf{d} = \sqrt[3]{\left(\frac{3\mathbf{F}}{4\,\mathrm{E}'\cdot\sqrt{\mathrm{R}}}\right)^2} \,\mathrm{in}\,\mathrm{[m]}\tag{3}$$

E' is the combined Young's modulus and can be calculated according to Equation (4) [24].

$$\mathbf{E}' = \frac{\mathbf{E}\_1 \cdot \mathbf{E}\_2}{\left(1 - \mathbf{v}\_1^2\right) \cdot \mathbf{E}\_2 + \left(1 - \mathbf{v}\_2^2\right) \cdot \mathbf{E}\_1} \text{ in } \left[\text{Nm}^{-2}\right] \tag{4}$$

The touching radius r can then be calculated with r = √ R · d [24]. It must be stated that using the touching radius r does not lead to the mathematical exact contact area as the contact area is not a flat circle but a small section of a sphere. However, the calculated indentation depths are by far smaller than the touching radius (d(r) = r <sup>2</sup>/R; e.g., d = 0.1 nm and r = 5.7 nm with rtip = 300 nm and tip force = 100 nN) and can thus be neglected. We, therefore, use the touching radius r for the calculation of the TCAs in SThM measurements (Section 4.3) as an adequate approximation for the simulations in Section 4.4 to Section 4.9.

#### *3.3. Heat Transfer*

Equation (5) represents the general heat conduction equation in the three-dimensional case without inner heat sources, in which the quotient k/(ρ · c) is expressed by the thermal diffusivity a [25].

$$\frac{\text{dT}}{\text{dt}} = \frac{\text{k}}{\rho \cdot \text{c}} \cdot \left( \frac{\partial^2 \text{T}}{\partial \text{x}^2} + \frac{\partial^2 \text{T}}{\partial \text{y}^2} + \frac{\partial^2 \text{T}}{\partial \text{z}^2} \right) = \text{ a} \cdot \Delta \text{T} \tag{5}$$

If there are inner heat sources h(t) (e.g., joule heating or chemical reactions), the inhomogeneous heat conduction Equation (6) follows (5) [25]. In this work, thermal conduction is applied to the entire geometry of the simulations in Section 4.4 to Section 4.9. CM solves said equations using numerical methods.

$$\frac{\text{dT}}{\text{dt}} = \left| \mathbf{a} \cdot \Delta \mathbf{T} + \mathbf{h}(\mathbf{t}) \right| \tag{6}$$

Convection is a mass bound energy transport caused by macroscopic particle movement. Considering free convection, this flow mainly results from the temperature-dependent particle movement and the thermal buoyancy. The reason for this is that fluids normally expand with increasing temperatures and therefore exhibit lower densities. Forced convection occurs if, e.g., a ventilator causes the flow field. Forced convection then overlaps with free convection. The influence of free convection is studied in Section 4.7 using Equation (7), where α is the heat transfer coefficient, which depends on the materials, surfaces and ambient conditions, Aconv the convective area and T<sup>2</sup> and T<sup>1</sup> the temperatures of the material and the surrounding fluid (e.g., air), respectively [25]. In this work, the influence of convection can be neglected. This is further discussed in Section 4.7.

$$\dot{\mathbf{Q}}\_{\rm conv} = \boldsymbol{\alpha} \cdot \mathbf{A}\_{\rm conv} \cdot \left(\mathbf{T}\_2 - \mathbf{T}\_1\right) \text{in } [\mathbf{W}] \tag{7}$$

Thermal radiation is not mass bound. Therefore, it occurs even under high vacuum conditions. Electromagnetic waves run from one body surface to another. Each body absorbs and emits radiation. In the special case of a small body with an area Arad and temperature T<sup>2</sup> surrounded by a much greater emissive area with a temperature T1, the net radiative heat flow of the body can be calculated with Equation (8), using the Stefan-Boltzmann law [25].

$$\dot{\mathbf{Q}}\_{\rm rad} = -\boldsymbol{\varepsilon}\_{\rm rad} \cdot \boldsymbol{\sigma} \cdot \mathbf{A}\_{\rm rad} \cdot \left(\mathbf{T}\_2^4 - \mathbf{T}\_1^4\right) \text{ in } [\mathcal{W}] \tag{8}$$

The negative radiation constant −rad depends on the considered body, and σ is the Stefan-Boltzmann constant. When T<sup>2</sup> is greater than T<sup>1</sup> the sign of . Qrad is negative which means, that the considered body loses energy. In this work, the influence of radiation of the heated SThM cantilever can be neglected. This is further discussed in Section 4.7.

#### **4. Results and Discussion**

#### *4.1. SEM Investigation of a Used SThM Tip*

To obtain vital information about the geometry of the cantilever and tip, we performed an SEM investigation of the thermal SThM probe Bruker VITA-DM-NANOTA-200 [20], which was already in use. This facilitated modeling the probe in detail in SW and paved the way for the simulations. SEM investigations require a conductive connection between the sample and the holder to avoid electrical overcharging caused by trapped electrons from the electron beam. Therefore, the probe was fixed on a sample holder using conductive silver and gold in a sputtering process. Figure 1a–c shows the probe before and after the preparation process, as well as inside the sample chamber of the SEM. To collect precise geometry data, the probe was investigated in two positions, one side view and one top view. The results are depicted in Figure 1d–i. Therefore, the tip radius was estimated to be approximately 300 nm. The geometry data for the simulations were created upon these images and are depicted in Figure 1j.

We assume a tip radius of 300 nm to represent an "average" used probe. However, an "average" used probe can hardly be defined as it depends on various factors such as the number of measurements, tip velocities and forces, sample materials, and many more. Mostly, SThM probes are used as long as they deliver reliable results, and from our experience, such tips exhibit tip radii within the region of 300 nm. Furthermore, new tips may degrade much faster than used tips, which is the reason why larger tip radii of around 300 nm may occur more often in practical SThM measurements. As the main residue is not located directly at the TCA (see Figure 1f), we assume it to have no influence. However, residues may increase the TIR between tip and sample. This effect is further discussed in Section 4.8.

**Figure 1.** (**a**) Photography image of the scanning thermal microscopy (SThM) probe Bruker VITA-DM-NANOTA-200 before preparation for SEM. (**b**) SEM image of a fully prepared tip. The sputtering area, wires for the heating current, and the cantilever are visible. (**c**) Image of the sample fixation inside the SEM. (**d**) Side-view close-up of the tip in (**e**) and (**f**). (**g**) Top-view close-up of the tip in (**h**) and (**i**). (**j**) Dimensions of the probe used for the simulations and isometric view. This SEM investigation made it possible to estimate the geometry of the cantilever and the radius of the used tip, which is approximately 300 nm. **Figure 1.** (**a**) Photography image of the scanning thermal microscopy (SThM) probe Bruker VITA-DM-NANOTA-200 before preparation for SEM. (**b**) SEM image of a fully prepared tip. The sputtering area, wires for the heating current, and the cantilever are visible. (**c**) Image of the sample fixation inside the SEM. (**d**) Side-view close-up of the tip in (**e**,**f**). (**g**) Top-view close-up of the tip in (**h**,**i**). (**j**) Dimensions of the probe used for the simulations and isometric view. This SEM investigation made it possible to estimate the geometry of the cantilever and the radius of the used tip, which is approximately 300 nm.

#### *4.2. Cantilever Displacement and Von Mises Stress 4.2. Cantilever Displacement and Von Mises Stress*

To study the displacement and the von Mises stress of the cantilever, a 3D-simulation was set up. The cantilever (geometry see Figure 1j) is fixed by the two connectors on one end. The simulated load is applied directly to the cantilever tip on the other end. All other edges and areas are set to be free. Moreover, gravity is set to a value of 9*.*81 ms−2 in the zdirection. For the cantilever, the material bulk Si was used with Young's modulus of 131 GPa and a Poisson's ratio of 0.221 [26]. The mesh was built with the option *Physics-*To study the displacement and the von Mises stress of the cantilever, a 3D-simulation was set up. The cantilever (geometry see Figure 1j) is fixed by the two connectors on one end. The simulated load is applied directly to the cantilever tip on the other end. All other edges and areas are set to be free. Moreover, gravity is set to a value of 9.81 ms−<sup>2</sup> in the z-direction. For the cantilever, the material bulk Si was used with Young's modulus of 131 GPa and a Poisson's ratio of 0.221 [26]. The mesh was built with the option *Physics-*

*controlled extremely fine*, which provides the finest automatically built mesh in CM. The study was performed in the stationary regime. *controlled extremely fine*, which provides the finest automatically built mesh in CM. The study was performed in the stationary regime.

Figure 2a illustrates the total displacement in µm, (b) represents the von Mises stress in Nm−<sup>2</sup> under a tip force of 100 nN, which acts directed towards the tip in the negative z-direction. Logically, the total displacement reaches its maximum value at the largest x-coordinate, whereas the von Mises stress achieves its maximum directly at the clamping. Here are some important simulation results for a tip force of 100 nN. The simulation results for 10 nN and 1000 nN can be derived by the multiplication of the following values with 0.1 and 10, respectively: Figure 2a illustrates the total displacement in µm, (b) represents the von Mises stress in Nm−2 under a tip force of 100 nN, which acts directed towards the tip in the negative zdirection. Logically, the total displacement reaches its maximum value at the largest xcoordinate, whereas the von Mises stress achieves its maximum directly at the clamping. Here are some important simulation results for a tip force of 100 nN. The simulation results for 10 nN and 1000 nN can be derived by the multiplication of the following values with 0.1 and 10, respectively:


**Figure 2.** Simulation results under a tip force of 100 nN. (**a**) Total displacement in µm. (**b**) Von Mises stress in Nm<sup>−</sup>2. **Figure 2.** Simulation results under a tip force of 100 nN. (**a**) Total displacement in µm. (**b**) Von Mises stress in Nm−<sup>2</sup> .

These values were also calculated manually for a simplified cantilever with a rectangular cross-section (b = 24.5 µm, h = 2 µm and l = 194.5 µm) fixed on only one end using the Equations (1) and (2). The results of the calculations are w (l = 194.5 µm) = 0.115 µm and σb = 1.19 × 106 Nm−2. Compared to the simulation results above (0.114 µm and 1*.*26 × 106 Nm−2), it is obvious that the exact geometry around the tip has just little influence on the displacement and the maximum von Mises stress. The breaking strength of Si is assumed to be in the range from 5 × 107 Nm−2 to 20 × 107 Nm−2 [27]. Compared to the maximum von Mises stress at the fixed end of the cantilever (1*.*26 × 106 Nm−2), there is a safety factor of approximately 40 to 160 before the cantilever breaks. Based on this calculation, the maximum tip forces that can be applied to the tip without breaking the cantilever are in the range from 4 µN to 16 µN. However, such values will not be reached in reasonable These values were also calculated manually for a simplified cantilever with a rectangular cross-section (b = 24.5 µm, h = 2 µm and l = 194.5 µm) fixed on only one end using the Equations (1) and (2). The results of the calculations are w (l = 194.5 µm) = 0.115 <sup>µ</sup>m and <sup>σ</sup><sup>b</sup> = 1.19 <sup>×</sup> <sup>10</sup><sup>6</sup> Nm−<sup>2</sup> . Compared to the simulation results above (0.114 µm and 1.26 <sup>×</sup> <sup>10</sup><sup>6</sup> Nm−<sup>2</sup> ), it is obvious that the exact geometry around the tip has just little influence on the displacement and the maximum von Mises stress. The breaking strength of Si is assumed to be in the range from 5 <sup>×</sup> <sup>10</sup><sup>7</sup> Nm−<sup>2</sup> to 20 <sup>×</sup> <sup>10</sup><sup>7</sup> Nm−<sup>2</sup> [27]. Compared to the maximum von Mises stress at the fixed end of the cantilever (1.26 <sup>×</sup> <sup>10</sup><sup>6</sup> Nm−<sup>2</sup> ), there is a safety factor of approximately 40 to 160 before the cantilever breaks. Based on this calculation, the maximum tip forces that can be applied to the tip without breaking the cantilever are in the range from 4 µN to 16 µN. However, such values will not be reached in reasonable SThM measurements.

#### SThM measurements. *4.3. Calculation of TCAs in Realistic SThM Investigations*

lated thermal contact areas are depicted in Figure 3.

*4.3. Calculation of TCAs in Realistic SThM Investigations*  To estimate TCAs between an SThM probe and different samples, we performed a mathematical evaluation regarding Equations (3) and (4). We calculated six different cases: a new Si probe with a tip radius of 30 nm and an "average" used to probe with a tip radius of 300 nm each combined with the Si sample, SiO2 sample and h-BN sample. The tip radius (300 nm) for an "average" used probe was determined by SEM investigations (Section 4.1). For the calculations, we used the material properties in Table 1. The calcu-To estimate TCAs between an SThM probe and different samples, we performed a mathematical evaluation regarding Equations (3) and (4). We calculated six different cases: a new Si probe with a tip radius of 30 nm and an "average" used to probe with a tip radius of 300 nm each combined with the Si sample, SiO<sup>2</sup> sample and h-BN sample. The tip radius (300 nm) for an "average" used probe was determined by SEM investigations (Section 4.1). For the calculations, we used the material properties in Table 1. The calculated thermal contact areas are depicted in Figure 3.

**Table 1.** Material properties are used for the theoretical calculation of the thermal contact areas (TCAs) in Figure 3. Values for SiO<sup>2</sup> are approximated from the predefined material in COMSOL Multiphysics (CM). **Table 1.** Material properties are used for the theoretical calculation of the thermal contact areas (TCAs) in Figure 3. Values for SiO2 are approximated from the predefined material in COMSOL


*Nanomaterials* **2021**, *11*, 491 8 of 21

**Figure 3.** Calculated TCAs in realistic SThM measurements for different probe-sample combinations. The red curves represent a new probe with a tip radius of 30 nm, and the blue curves represent an "average" used probe with a tip radius of 300 nm. A similar investigation can be found in [30]. **Figure 3.** Calculated TCAs in realistic SThM measurements for different probe-sample combinations. The red curves represent a new probe with a tip radius of 30 nm, and the blue curves represent an "average" used probe with a tip radius of 300 nm. A similar investigation can be found in [30].

We chose the above described tip-sample combinations for reasons of reusability as these configurations are comparatively common and might also occur in similar SThM investigations of other researchers. We also used the calculated TCAs of the h-BN curves We chose the above described tip-sample combinations for reasons of reusability as these configurations are comparatively common and might also occur in similar SThM investigations of other researchers. We also used the calculated TCAs of the h-BN curves to calculate the thermal contact radius (TCR) for the simulations in Section 4.4 to Section 4.9.

to calculate the thermal contact radius (TCR) for the simulations in sections 4.4 to 4.9. Appropriate material properties for h-BN are not easy to estimate, as in ultrathin (2D) materials, they are also strongly dependent on factors such as material orientation and exact film thickness. However, in our calculations concerning the h-BN sample, the Poisson's ratio has just little influence on the calculated contact area. The main factor is Young's modulus, as it is assumed to be extremely high in comparison to Si and SiO2. Consequently, if other researchers might find different values of Young's modulus of h-BN appropriate for their individual case, the contact areas can be approximated to lie in between the h-BN and SiO2 curves in Figure 3. Prerequisites are a similar measurement setup and Young's modulus between 70 and 850 GPa. Nevertheless, it must be stated that Appropriate material properties for h-BN are not easy to estimate, as in ultrathin (2D) materials, they are also strongly dependent on factors such as material orientation and exact film thickness. However, in our calculations concerning the h-BN sample, the Poisson's ratio has just little influence on the calculated contact area. The main factor is Young's modulus, as it is assumed to be extremely high in comparison to Si and SiO2. Consequently, if other researchers might find different values of Young's modulus of h-BN appropriate for their individual case, the contact areas can be approximated to lie in between the h-BN and SiO<sup>2</sup> curves in Figure 3. Prerequisites are a similar measurement setup and Young's modulus between 70 and 850 GPa. Nevertheless, it must be stated that the lower Young's modulus, the higher becomes the influence of different values for the Poisson's ratio.

the lower Young's modulus, the higher becomes the influence of different values for the Poisson's ratio. For our simulations, we assume the calculated ideal contact areas in this section to be For our simulations, we assume the calculated ideal contact areas in this section to be realistic, as the impact of surface roughness decreases with lower surface roughnesses. We also assume roughness to be neglectable as we work with samples exhibiting surface

realistic, as the impact of surface roughness decreases with lower surface roughnesses.

roughnesses in the sub-nm area. In the following sections, we use the calculated TCAs to simulate the heat transfer in SThM measurements. roughnesses in the sub-nm area. In the following sections, we use the calculated TCAs to simulate the heat transfer in SThM measurements.

We also assume roughness to be neglectable as we work with samples exhibiting surface

#### *4.4. Heat Spreading in SThM Measurements Regarding Ultrathin h-BN Film 4.4. Heat Spreading in SThM Measurements Regarding Ultrathin h-BN Film*

*Nanomaterials* **2021**, *11*, 491 9 of 21

To study the heat spreading in SThM measurements on a sample with an ultrathin h-BN film on top, this simulation was set up. To study the heat spreading in SThM measurements on a sample with an ultrathin h-BN film on top, this simulation was set up.

**Simulation setup:** For the cantilever (dimensions see Figure 1j), the CM-predefined material bulk Si was set. The rectangular sample (footprint 100 µm × 100 µm) consists of the CM-predefined material h-BN with a thickness of 10 nm on top of 10 µm thick bulk Si. Regarding the anisotropic thermal properties of h-BN, we assume a cross-plane thermal conductivity of h-BN of 1 Wm−1K−<sup>1</sup> . The in-plane thermal conductivity of h-BN was defined using the parameter *k*, which varies depending on the simulations [31]. The thermal contact resistance between h-BN and bulk Si is assumed to be 4 <sup>×</sup> <sup>10</sup>−<sup>8</sup> Km2W−<sup>1</sup> . However, this value can only be estimated and will vary in practical measurements depending on material quality, ambient conditions and manufacturing processes. We chose this value inspired by investigations on similar material stacks such as graphene/Si [32] since the atomic structure of graphene is similar to h-BN. Villaroman et al. estimate the thermal contact resistance between graphene and Si to 3.1–4.9 <sup>×</sup> <sup>10</sup>−<sup>8</sup> Km2W−<sup>1</sup> [32]. Values in the same order of magnitude can be found in [33,34]. **Simulation setup:** For the cantilever (dimensions see Figure 1j), the CM-predefined material bulk Si was set. The rectangular sample (footprint 100 µm × 100 µm) consists of the CM-predefined material h-BN with a thickness of 10 nm on top of 10 µm thick bulk Si. Regarding the anisotropic thermal properties of h-BN, we assume a cross-plane thermal conductivity of h-BN of 1 Wm−1K−1. The in-plane thermal conductivity of h-BN was defined using the parameter , which varies depending on the simulations [31]. The thermal contact resistance between h-BN and bulk Si is assumed to be 4 × 10−8 Km2W−1. However, this value can only be estimated and will vary in practical measurements depending on material quality, ambient conditions and manufacturing processes. We chose this value inspired by investigations on similar material stacks such as graphene/Si [32] since the atomic structure of graphene is similar to h-BN. Villaroman et al. estimate the thermal contact resistance between graphene and Si to 3.1–4.9 × 10−8 Km2W−1 [32]. Values in the same order of magnitude can be found in [33,34].

The initial value for the temperature was set to 293.15 K for all boundaries. Between tip and sample, an ideal thermal contact was defined. We calculated the TCR of the TCA as follows: With Figure 3 (blue squared h-BN curve), we assumed a realistic TCA of ~103.4 nm<sup>2</sup> for rtip = 300 nm and a tip force of 100 nN, which leads to a TCR of ~5.7 nm (TCA = TCR<sup>2</sup> · Π). This appears to be a realistic value for an ideal thermal contact between a used SThM probe and a sample with an ultrathin h-BN film on top. The temperature of the topside of the cantilever is defined using the parameter *temp*, which varies between 50 ◦C and 200 ◦C depending on the simulations. These are appropriate cantilever temperatures for practical SThM investigations [12]. The remaining outer areas of the cantilever were defined as thermal isolating as well as the top surface of the sample, while the remaining outer areas of the sample take over the ambient temperature of exactly 293.15 K. An overview of the simulation setup and boundary conditions is given in Figure 4. The center point of the TCA defines the origin of the coordinate system. The initial value for the temperature was set to 293*.*15 K for all boundaries. Between tip and sample, an ideal thermal contact was defined. We calculated the TCR of the TCA as follows: With Figure 3 (blue squared h-BN curve), we assumed a realistic TCA of ~103.4 nm2 for rtip = 300 nm and a tip force of 100 nN, which leads to a TCR of ~5.7 nm (TCA = TCR<sup>ଶ</sup> ⋅ Π). This appears to be a realistic value for an ideal thermal contact between a used SThM probe and a sample with an ultrathin h-BN film on top. The temperature of the topside of the cantilever is defined using the parameter *temp*, which varies between 50 °C and 200 °C depending on the simulations. These are appropriate cantilever temperatures for practical SThM investigations [12]. The remaining outer areas of the cantilever were defined as thermal isolating as well as the top surface of the sample, while the remaining outer areas of the sample take over the ambient temperature of exactly 293*.*15 K. An overview of the simulation setup and boundary conditions is given in Figure 4. The center point of the TCA defines the origin of the coordinate system.

**Figure 4.** Simulation setup and boundary conditions. Left: The general setup consists of the heated SThM probe and the sample. Right: zoom into the TCA with ideal thermal contact between tip and sample. **Figure 4.** Simulation setup and boundary conditions. Left: The general setup consists of the heated SThM probe and the sample. Right: zoom into the TCA with ideal thermal contact between tip and sample.

**Meshing strategy:** To simulate the thermal contact with sufficient resolution, a usercontrolled mesh was created. We defined the maximum dimension of a triangular mesh element for the smallest area of the entire simulation, the TCA of the cantilever, to 0.2 nm. Moreover, a circle with a radius of 500 nm around the TCA was defined, in which the **Meshing strategy:** To simulate the thermal contact with sufficient resolution, a usercontrolled mesh was created. We defined the maximum dimension of a triangular mesh element for the smallest area of the entire simulation, the TCA of the cantilever, to 0.2 nm. Moreover, a circle with a radius of 500 nm around the TCA was defined, in which the maximum dimension of a mesh element is 5 nm. For the remaining geometry of the

cantilever and sample, we used the predefined meshing strategy with element size normal. In the present simulations, this meshing strategy represents a good tradeoff between accuracy and simulation time and delivers trustworthy results. The meshing strategy is evaluated in Section 4.9, in which the reliability of our results is demonstrated. Figure 5 illustrates the applied meshing strategy in more detail. maximum dimension of a mesh element is 5 nm. For the remaining geometry of the cantilever and sample, we used the predefined meshing strategy with element size normal. In the present simulations, this meshing strategy represents a good tradeoff between accuracy and simulation time and delivers trustworthy results. The meshing strategy is evaluated in Section 4.9, in which the reliability of our results is demonstrated. Figure 5 illustrates the applied meshing strategy in more detail.

**Figure 5.** Meshing strategy zooming from the geometry overview to the small TCA. **Figure 5.** Meshing strategy zooming from the geometry overview to the small TCA.

**Parametric studies:** To investigate the cross- and in-plane heat distribution for different cases, a parametric sweep study in the stationary regime was performed. Every single combination of the parameters *k* and *temp* was calculated, which are defined as **Parametric studies:** To investigate the cross- and in-plane heat distribution for different cases, a parametric sweep study in the stationary regime was performed. Every single combination of the parameters *k* and *temp* was calculated, which are defined as follows:


are specified in K. **Simplifications:** Simulations can never represent the real world as the results are only as good as the chosen simulation model. We tried to create the simulation models as realistic as possible. However, with regards to simulation time and the reusability of our results, we also had to introduce simplifications. The chosen simplifications do not decisively influence our results and are therefore justified. The main simplifications are listed **Simplifications:** Simulations can never represent the real world as the results are only as good as the chosen simulation model. We tried to create the simulation models as realistic as possible. However, with regards to simulation time and the reusability of our results, we also had to introduce simplifications. The chosen simplifications do not decisively influence our results and are therefore justified. The main simplifications are listed below:

below: • We assume an ideal thermal contact between tip and sample in Sections 4.4 to 4.7. TIRs in practical SThM measurements are quite hard to estimate as they depend on numerous factors such as surface roughness, material combination, vertical steps, tip material and geometry, contaminations, tip force, and some more. Due to these numerous influences, TIR will also vary greatly during a single SThM measurement. In literature, values in the range of 10−8 Km2W−1 to 10−10 Km2W−1 for different probes, samples, and measurement scenarios can be found [35–37]. Indeed, we assume the TIR to vary in a broader range due to the great number of possible measurement scenarios. To realize a comparison between the simulations in sections 4.4 and 4.5, we assume an ideal thermal contact in said sections. The influence of a possible TIR is further studied in Section 4.8. It must also be stated that the present investigation in Section 4.4 focuses on ultrathin h-BN samples, which in fact are super flat as they are built of a stack of single h-BN layers. Roughnesses of high-quality h-BN are assumed to be less than 0.4 nm [38], which is significantly lower compared to a tip • We assume an ideal thermal contact between tip and sample in Section 4.4 to Section 4.7. TIRs in practical SThM measurements are quite hard to estimate as they depend on numerous factors such as surface roughness, material combination, vertical steps, tip material and geometry, contaminations, tip force, and some more. Due to these numerous influences, TIR will also vary greatly during a single SThM measurement. In literature, values in the range of 10−<sup>8</sup> Km2W−<sup>1</sup> to 10−<sup>10</sup> Km2W−<sup>1</sup> for different probes, samples, and measurement scenarios can be found [35–37]. Indeed, we assume the TIR to vary in a broader range due to the great number of possible measurement scenarios. To realize a comparison between the simulations in Sections 4.4 and 4.5, we assume an ideal thermal contact in said sections. The influence of a possible TIR is further studied in Section 4.8. It must also be stated that the present investigation in Section 4.4 focuses on ultrathin h-BN samples, which in fact are super flat as they are built of a stack of single h-BN layers. Roughnesses of high-quality h-BN are assumed to be less than 0.4 nm [38], which is significantly lower compared to a tip radius of 300 nm. As surface roughness has a great impact on TIR, values for h-BN samples should be lower compared to samples with higher roughnesses;

• We neglect a possible water meniscus around the tip. Other researchers also propose that a water meniscus can often be neglected in SThM measurements [22]. On one hand, a present water drop could cause heat conduction and may increase the amount of heat flux between tip and sample slightly. On the other hand, SThM measurements can also be performed under high vacuum conditions, where the appearance of a water meniscus can be excluded. We, therefore, focus on simulations without a possible water meniscus; • We neglect a possible water meniscus around the tip. Other researchers also propose that a water meniscus can often be neglected in SThM measurements [22]. On one hand, a present water drop could cause heat conduction and may increase the amount of heat flux between tip and sample slightly. On the other hand, SThM measurements can also be performed under high vacuum conditions, where the appearance of a water meniscus can be excluded. We, therefore, focus on simulations without a possible water meniscus;

radius of 300 nm. As surface roughness has a great impact on TIR, values for h-BN

samples should be lower compared to samples with higher roughnesses;

*Nanomaterials* **2021**, *11*, 491 11 of 21


**Results:** To compare single measurements, it is necessary to define useful measurable parameters. Hence, we define the thermal active radius (TAR) located on the top surface of the sample starting at the geometric center point of the tip. The temperature until the TAR is greater than 294.15 K, which represents a temperature rise of 1 K in comparison to ambient conditions. In Figure 6b, the heat flow in the z-direction is investigated. It is obvious that with increasing *temp*, the temperature at TCA also increases. However, an increasing *k* results in a greater heat spreading in x and y-direction, which leads to a stronger temperature decrease within the ultrathin h-BN film and the tip (compare curves of the same color each in Figure 6b. It can also be observed that the heat penetration in z-direction ends at the z-coordinate 10 nm, which is exactly where the h-BN ends. The thermal contact resistance between h-BN and Si causes this rapid temperature decrease. Thus, no heat is transferred into the bulk Si. The green curve in Figure 6a illustrates the isothermal line of the TAR. With increasing *k* and *temp,* this line moves to greater values of x, which can also be seen in Figure 7b–d. **Results:** To compare single measurements, it is necessary to define useful measurable parameters. Hence, we define the thermal active radius (TAR) located on the top surface of the sample starting at the geometric center point of the tip. The temperature until the TAR is greater than 294.15 K, which represents a temperature rise of 1 K in comparison to ambient conditions. In Figure 6b, the heat flow in the z-direction is investigated. It is obvious that with increasing *temp,* the temperature at TCA also increases. However, an increasing *k* results in a greater heat spreading in x and y-direction, which leads to a stronger temperature decrease within the ultrathin h-BN film and the tip (compare curves of the same color each in Figure 6b. It can also be observed that the heat penetration in z-direction ends at the z-coordinate 10 nm, which is exactly where the h-BN ends. The thermal contact resistance between h-BN and Si causes this rapid temperature decrease. Thus, no heat is transferred into the bulk Si. The green curve in Figure 6a illustrates the isothermal line of the TAR. With increasing *k* and *temp,* this line moves to greater values of x, which can also be seen in Figure 7b–d.

**Figure 6.** Simulation regarding ultrathin (10 nm) h-BN film on top of Si. (**a**) Stationary heat distribution for *k* = 1 and *temp* = 100 °C on the x/z-plane through the center of TCA (y = 0). (**b**) Temperature curves along the z-coordinate through the center point of TCA for various simulated cases. **Figure 6.** Simulation regarding ultrathin (10 nm) h-BN film on top of Si. (**a**) Stationary heat distribution for *k* = 1 and *temp* = 100 ◦C on the x/z-plane through the center of TCA (y = 0). (**b**) Temperature curves along the z-coordinate through the center point of TCA for various simulated cases.

Figure 7a illustrates the heat distribution along the x-direction (y = z = 0). It can be seen that higher values of *k* result in a greater TAR. This can also be seen in Figure 7b–d for the special case *temp* = 100 °C. With an increasing *k,* the green isothermal circle line of the TAR moves to greater values (64 nm @ *k* = 1 and 232 nm @ *k* = 100), indicating that heat spreading on the x/y-plane of h-BN increases. Furthermore, we can see that the heat spreading effect between *k* = 1 and *k* = 10 is greater than between *k* = 10 and *k* = 100 (compare curves of the same color in Figure 7a. There seems to be a kind of saturation effect for an increasing *k*. Surely, the temperature of TCA increases with increasing *temp* (Figure 7a). A higher *k* leads to an effective greater thermal conductivity of h-BN. The small TCA Figure 7a illustrates the heat distribution along the x-direction (y = z = 0). It can be seen that higher values of *k* result in a greater TAR. This can also be seen in Figure 7b–d for the special case *temp* = 100 ◦C. With an increasing *k,* the green isothermal circle line of the TAR moves to greater values (64 nm @ *k* = 1 and 232 nm @ *k* = 100), indicating that heat spreading on the x/y-plane of h-BN increases. Furthermore, we can see that the heat spreading effect between *k* = 1 and *k* = 10 is greater than between *k* = 10 and *k* = 100 (compare curves of the same color in Figure 7a. There seems to be a kind of saturation effect for an increasing *k*. Surely, the temperature of TCA increases with increasing *temp* (Figure 7a). A higher *k* leads to an effective greater thermal conductivity of h-BN. The small TCA then has an increased proportion of the entire thermal resistance, which is the reason for lower temperatures at the TCA with *temp* = const. and *k* increasing.

*Nanomaterials* **2021**, *11*, 491 13 of 21

**Figure 7.** Simulation regarding ultrathin (10 nm) hexagonal boron nitride (h-BN) film on top of Si. (**a**) Temperature curves along the x-coordinate on the top surface (y = z = 0) for different simulated cases. (**b**) to (**d**) Heat distribution and thermal active radius (TAR) on x/y-plane for *temp* = 100 °C and *k* = 1 in (**b**), *k* = 10 in (**c**) and *k* = 100 in (**d**). (**e**) TAR in dependency of *temp* and *k* for all simulated cases. (**f**) Interpolated color plot of the simulated TARs in (**e**). (**g**) normal total energy flux (NTEF) through the TCA in dependency of *temp* and *k* for all simulated cases. (**h**) Interpolated color plot of the simulated NTEFs in (**g**). **Figure 7.** Simulation regarding ultrathin (10 nm) hexagonal boron nitride (h-BN) film on top of Si. (**a**) Temperature curves along the x-coordinate on the top surface (y = z = 0) for different simulated cases. (**b**) to (**d**) Heat distribution and thermal active radius (TAR) on x/y-plane for *temp* = 100 ◦C and *k* = 1 in (**b**), *k* = 10 in (**c**) and *k* = 100 in (**d**). (**e**) TAR in dependency of *temp* and *k* for all simulated cases. (**f**) Interpolated color plot of the simulated TARs in (**e**). (**g**) normal total energy flux (NTEF) through the TCA in dependency of *temp* and *k* for all simulated cases. (**h**) Interpolated color plot of the simulated NTEFs in (**g**).

A comparison of the NTEFs to measured values of other researchers is quite hard. There is an almost infinite number of possible measurement setups, and our specific setup could not be found in literature so far, to the best of our knowledge. However, the following references used different setups but may verify our simulated range of the NTEFs. Figure 7e shows the simulated values for the TAR as a graphical representation. The TAR increases with an increasing *temp,* but there seems to be a kind of saturation effect. The TAR logically also increases with higher *k* as heat spreading is greater there. Values for

Hwang et al. measured heat fluxes during null-point SThM in the range of approximately

the TAR lie in between 46 nm (@ *temp* = 50 ◦C and *k* = 1) and 357 nm (@ *temp* = 100 ◦C and *k* = 100). Figure 7f delivers a color plot of the TARs in dependency of *k* (x-axes) and *temp* (*y*-axes). This color plot was created using the simulated values in Figure 7e). Intermediate points are interpolated. Depending on the case, other researchers can estimate the TAR for their research using Figure 7a–f. The thermal penetration depths for all cases is exactly 10 nm.

We also simulated the normal total energy flux (NTEF) in the stationary state, which flows through the small TCA for every case. Figure 7g illustrates the simulated cases. The NTEF seems to increase nearly linearly for an increasing *temp*. The NTEF also increases with greater *k*. Values for the NTEF vary from 0.6 µW (@ *temp* = 50 ◦C and *k* = 1) to 32.9 µW (@ *temp* = 100 ◦C and *k* = 100). Figure 7h shows a color plot of the NTEF in dependency of *k* (x-axes) and *temp* (y-axes). This color plot was created using the simulated values in Figure 7g. Intermediate points are interpolated. It, therefore, allows a rough estimation for other research for different cases.

A comparison of the NTEFs to measured values of other researchers is quite hard. There is an almost infinite number of possible measurement setups, and our specific setup could not be found in literature so far, to the best of our knowledge. However, the following references used different setups but may verify our simulated range of the NTEFs. Hwang et al. measured heat fluxes during null-point SThM in the range of approximately 1 µW (Teflon-coated surface) and 6 µW (SiO<sup>2</sup> surface) using a thermocouple probe [39]. Assy and Gomès used a Wollaston wire probe at 140 ◦C and a Kelvin nanotechnology probe at 65 ◦C on germanium and Si samples. They calculated the probe Joule power relative difference ∆P/P = (P<sup>c</sup> − Pa)/Pc, where P<sup>c</sup> and P<sup>a</sup> represent the probe Joule power in contact and out of contact, respectively. Unfortunately, only relative values of ∆P/P ranging from 0.003 to 0.058 are presented [40].

#### *4.5. Heat Spreading in SThM Measurements Regarding a Bulk Si Sample*

To compare the results of Section 4.4 to a realistic SThM measurement with a bulk Si sample, this simulation was set up. The only difference to the simulation in Section 4.4 is that the 10 nm thick h-BN layer is removed. The main results are presented in Figure 8.

Figure 8a shows the heat distribution for *temp* = 100 ◦C on the x/z-plane, whereas (c) shows the same heat distribution on the x/y-plane. Since the Si sample features isotropic thermal conductivity, the heat spreading is more or less circular around the tip. The green line is the TAR-line, where the temperature rise regarding the ambient temperature of 293.15 K is 1 K. As a reason of the mesh density, this line is not exactly a semicircle. The ideal semicircle is represented by the red dotted line, which overlaps the green semicircle. It can be seen that there is only a little deviation from the ideal semicircle caused by the mesh density. Compared to the h-BN sample (Section 4.4), we can see that the penetration depth is not limited by a TCR at 10 nm, and therefore, the penetration depth is much larger than for the h-BN sample. Figure 8b illustrates the temperature drop in z-direction exactly through the center of the TCA for all simulated cases. It is obvious that with increasing *temp,* the temperature at the TCA and in general also increases. Temperatures at z = 0 nm, which is directly at the physical center point of the TCA, are in between 301.8 K (@ *temp* = 50 ◦C) and 336.1 K (@ *temp* = 200 ◦C). Compared to the h-BN sample, we can see that the temperature drop within the cantilever is much higher. This is because Si has a higher overall thermal conductivity than h-BN, whereby the thermal resistance of the small TCA creates a greater proportion of the entire thermal resistance. Figure 8d shows the temperature curves along the x-direction (y = z = 0). As expected, with an increasing *temp* the curves also rise.

Figure 8e illustrates the TAR and the NTEF in dependency of *temp* for all simulated cases of the Si sample and compares them to the simulations of the isotropic case (*k* = 1) of the h-BN sample in Section 4.4. It is obvious that the NTEFs are significantly greater for the Si sample and increase with increasing *temp*. In general, the NTEFs for the Si sample are in between 28.3 µW (@ *temp* = 50 ◦C) and 126.8 µW (@ *temp* = 200 ◦C), the TARs are ranging

from 29 nm (@ *temp* = 50 ◦C) to 132 nm (@ *temp* = 200 ◦C). The TARs of the Si sample are greater compared to the h-BN sample from *temp* ≈ 90 ◦C upwards. Values for the NTEF are also significantly greater compared to the maximum curves of the h-BN sample with *k* = 100 in Figure 7g. The root cause for this is the greater overall thermal conductivity of Si. *4.5. Heat Spreading in SThM Measurements Regarding a Bulk Si Sample*  To compare the results of Section 4.4 to a realistic SThM measurement with a bulk Si sample, this simulation was set up. The only difference to the simulation in Section 4.4 is that the 10 nm thick h-BN layer is removed. The main results are presented in Figure 8.

1 µW (Teflon-coated surface) and 6 µW (SiO2 surface) using a thermocouple probe [39]. Assy and Gomès used a Wollaston wire probe at 140 °C and a Kelvin nanotechnology probe at 65 °C on germanium and Si samples. They calculated the probe Joule power relative difference ∆P/P = (Pc − Pa)/Pc, where Pc and Pa represent the probe Joule power in contact and out of contact, respectively. Unfortunately, only relative values of ∆P/P rang-

*Nanomaterials* **2021**, *11*, 491 14 of 21

ing from 0.003 to 0.058 are presented [40].

**Figure 8.** Simulation regarding bulk Si sample. (**a**) Stationary heat distribution for *temp* = 100 °C on the x/z-plane through the center of the TCA (y = 0). (**b**) Temperature curves along the z-coordinate through the center point of the TCA for different simulated cases. (**c**) Heat distribution and TAR on x/y-plane for *temp* = 100 °C. (**d**) Temperature curves along the x-coordinate on the top surface (y = z = 0) for different simulated cases. (**e**) TAR and NTEF in dependency of *temp* for all simulated cases of the Si sample and comparison to simulations with the h-BN sample (isotropic case with *k* = 1). **Figure 8.** Simulation regarding bulk Si sample. (**a**) Stationary heat distribution for *temp* = 100 ◦C on the x/z-plane through the center of the TCA (y = 0). (**b**) Temperature curves along the z-coordinate through the center point of the TCA for different simulated cases. (**c**) Heat distribution and TAR on x/y-plane for *temp* = 100 ◦C. (**d**) Temperature curves along the x-coordinate on the top surface (y = z = 0) for different simulated cases. (**e**) TAR and NTEF in dependency of *temp* for all simulated cases of the Si sample and comparison to simulations with the h-BN sample (isotropic case with *k* = 1).

#### Figure 8a shows the heat distribution for *temp* = 100 °C on the x/z-plane, whereas (c) shows the same heat distribution on the x/y-plane. Since the Si sample features isotropic *4.6. Stationary Time of Heat Distribution*

To study the time until heat dissipation is stationary (tstat), we performed the subsequent simulations. The simulations are based on Sections 4.4 and 4.5, with the difference being time-dependent instead of stationary. We consider the minimum and maximum case regarding the NTEF through the TCA for both the h-BN and the Si sample. We consider the time-dependent temperature curve for the geometric center point of the TCA. Those four curves are depicted in Figure 9.

curves also rise.

conductivity of Si.

*4.6. Stationary Time of Heat Distribution* 

Those four curves are depicted in Figure 9.

**Figure 9.** Time-dependent simulation of the temperature rise at the center point of the TCA for the minimum and maximum cases of the h-BN and Si sample, respectively. **Figure 9.** Time-dependent simulation of the temperature rise at the center point of the TCA for the minimum and maximum cases of the h-BN and Si sample, respectively.

thermal conductivity, the heat spreading is more or less circular around the tip. The green line is the TAR-line, where the temperature rise regarding the ambient temperature of 293.15 K is 1 K. As a reason of the mesh density, this line is not exactly a semicircle. The ideal semicircle is represented by the red dotted line, which overlaps the green semicircle. It can be seen that there is only a little deviation from the ideal semicircle caused by the mesh density. Compared to the h-BN sample (Section 4.4), we can see that the penetration depth is not limited by a TCR at 10 nm, and therefore, the penetration depth is much larger than for the h-BN sample. Figure 8b illustrates the temperature drop in z-direction exactly through the center of the TCA for all simulated cases. It is obvious that with increasing *temp,* the temperature at the TCA and in general also increases. Temperatures at z = 0 nm, which is directly at the physical center point of the TCA, are in between 301.8 K (@ *temp* = 50 °C) and 336.1 K (@ *temp* = 200 °C). Compared to the h-BN sample, we can see that the temperature drop within the cantilever is much higher. This is because Si has a higher overall thermal conductivity than h-BN, whereby the thermal resistance of the small TCA creates a greater proportion of the entire thermal resistance. Figure 8d shows the temperature curves along the x-direction (y = z = 0). As expected, with an increasing *temp* the

Figure 8e illustrates the TAR and the NTEF in dependency of *temp* for all simulated cases of the Si sample and compares them to the simulations of the isotropic case (*k* = 1) of the h-BN sample in Section 4.4. It is obvious that the NTEFs are significantly greater for the Si sample and increase with increasing *temp*. In general, the NTEFs for the Si sample are in between 28.3 µW (@ *temp* = 50 °C) and 126.8 µW (@ *temp* = 200 °C), the TARs are ranging from 29 nm (@ *temp* = 50 °C) to 132 nm (@ *temp* = 200 °C). The TARs of the Si sample are greater compared to the h-BN sample from *temp* ≈ 90 °C upwards. Values for the NTEF are also significantly greater compared to the maximum curves of the h-BN sample with *k* = 100 in Figure 7g. The root cause for this is the greater overall thermal

To study the time until heat dissipation is stationary (tstat), we performed the subsequent simulations. The simulations are based on Section 4.4 and 4.5, with the difference being time-dependent instead of stationary. We consider the minimum and maximum case regarding the NTEF through the TCA for both the h-BN and the Si sample. We consider the time-dependent temperature curve for the geometric center point of the TCA.

As a result, we can say that the tstat for every simulated case is below 1 µs. In general, it seems that higher cantilever temperatures lead to a slightly greater tstat. In practical SThM measurements, the cantilever is usually moving over the sample with a specific velocity. This is not representable in simulations. The tip velocity in practical SThM measurements is normally below 20 µms−<sup>1</sup> [12], which means that in 1 µs, the tip moves less than 20 pm. This is by far smaller than the atomic radius of Si or the lattice constant of h-BN. Therefore, simulations in the stationary regime in Sections 4.4 and 4.5, which assume a motionless cantilever, are justified.

#### *4.7. Influence of Convection and Radiation*

To study the influence of radiation, we set up the same simulation as in Section 4.5, with the only difference of surface radiation being enabled for all areas which were isolated in Section 4.5. Surface emissivity rad of Si was set to 0.67 [41]. We also consider the minimum and maximum cases regarding the NTEF. By comparing the NTEF and the TAR for the minimum and maximum cases, we obtain the same results whether radiation being enabled or disabled (28.3 µW and 29 nm @ *temp* = 50 ◦C; 126.8 µW and 132 nm @ *temp* = 200 ◦C; compare also Figure 8e). Therefore, we can say that the influence of radiation can be neglected in our simulations. To roughly estimate the amount of radiative heat losses, a radiative area Arad with a radius of 100 nm and a constant surface temperature T<sup>2</sup> of 393.15 K with an ambient temperature T<sup>1</sup> of 293.15 K and an rad of 0.67 is considered. Using Equation (8), we calculate radiative heat losses of ~19.7 pW.

To roughly estimate the influence of convection a manual calculation with a convective heat transfer coefficient <sup>α</sup> of <sup>−</sup>5 Wm−<sup>2</sup> <sup>K</sup>−<sup>1</sup> [42], a convective area Aconv with a radius of 100 nm, a constant surface temperature T<sup>2</sup> of 393.15 K and an ambient temperature T<sup>1</sup> of 293.15 K is considered. Using Equation (7), we reach a loss of power caused by free convection of ~15.7 pW.

Finally, we may state that estimated convective and radiative heat losses in this work are extremely low compared to the simulated NTEF in Figures 7g and 8e. A comparative simulation also shows no differences regarding the NTEF and the TAR with radiation being enabled or disabled. Therefore, we neglect the influence of convection and radiation in our simulations. Furthermore, in the literature, the influence of radiation and convection is neglected in special cases of SThM measurements [22].

#### *4.8. Influence of the TIR at the TCA regarding the TAR and the NTEF*

To study the influence of a possible TIR between probe and sample directly at the TCA, these simulations were set up. TIRs in practical SThM measurements are quite hard to estimate as they depend on numerous factors such as surface roughness, material combination, vertical steps, tip material and geometry, contaminations, residue, tip force, and some more. Due to these numerous influences, the TIR will also vary significantly during a single SThM measurement. In literature, values in the range of 10−<sup>8</sup> Km2W−<sup>1</sup> to 10−<sup>10</sup> Km2W−<sup>1</sup> for different probes, samples, and measurement cases can be found [35–37]. Indeed, we assume the TIR to vary in a broader range due to the great number of possible measurement scenarios. To enable a comparison of the simulations in Sections 4.4 and 4.5, we assume an ideal thermal contact in said sections. In contrast, here, the influence of a varying TIR is investigated. Figure 10 illustrates the simulation results. We created the same simulations as in Section 4.4 (with *k* = 1) and Section 4.5 for *temp* = 50 ◦C and *temp* = 200 ◦C with the difference of varying the TIR from 5 <sup>×</sup> <sup>10</sup>−<sup>7</sup> Km2W−<sup>1</sup> to 1 <sup>×</sup> <sup>10</sup>−<sup>13</sup> Km2W−<sup>1</sup> . *Nanomaterials* **2021**, *11*, 491 17 of 21

**Figure 10.** TAR and NTEF in dependency of the TIR at the TCA. (**a**) Sample: h-BN (compare Section 4.4). (**b**) Sample: Si (compare Section 4.5). **Figure 10.** TAR and NTEF in dependency of the TIR at the TCA. (**a**) Sample: h-BN (compare Section 4.4). (**b**) Sample: Si(compare Section 4.5).

The case with an ultrathin h-BN film on the top surface is demonstrated in Figure 10a, representing ultrathin films with low thermal conductivities. It can be deduced that both, the TAR and the NTEF increase with decreasing TIR and converge exactly at the same values for the corresponding case in Section 4.4 (46 nm and 0.6 µW @ *temp* = 50 °C; 79 nm and 3.8 µW @ *temp* = 200 °C). This convergence starts approximately from 1 × 10−<sup>10</sup> Km2W−1 downwards. Figure 10b represents the case of the bulk Si sample as a sample with high thermal conductivity. Values for the TAR and the NTEF also converge at similar values compared to the corresponding cases in Section 4.5 (29 nm and 28.3 µW @ *temp* = 50 °C; 132 nm and 126.8 µW @ *temp* = 200 °C). Compared to the h-BN sample, this convergence effect starts with significantly lower TIR values at approximately 1×10−12 Km2W−1. However, it seems that in samples with a higher thermal conductivity, the influence of the TIR also increases. Nevertheless, it must be stated that SThM measurements with different material combinations cannot be compared directly, as the TIR depends on numerous influence factors as listed above. As a result, we may conclude that a TIR at the TCA has some influence and will reduce the ideal simulated values of the TAR and the NTEF in sections 4.4 and 4.5. This effect should be considered by other researchers assuming specific values for the TIR. The case with an ultrathin h-BN film on the top surface is demonstrated in Figure 10a, representing ultrathin films with low thermal conductivities. It can be deduced that both, the TAR and the NTEF increase with decreasing TIR and converge exactly at the same values for the corresponding case in Section 4.4 (46 nm and 0.6 µW @ *temp* = 50 ◦C; 79 nm and 3.8 <sup>µ</sup>W @ *temp* = 200 ◦C). This convergence starts approximately from <sup>1</sup> <sup>×</sup> <sup>10</sup>−<sup>10</sup> Km2W−<sup>1</sup> downwards. Figure 10b represents the case of the bulk Si sample as a sample with high thermal conductivity. Values for the TAR and the NTEF also converge at similar values compared to the corresponding cases in Section 4.5 (29 nm and 28.3 µW @ *temp* = 50 ◦C; 132 nm and 126.8 µW @ *temp* = 200 ◦C). Compared to the h-BN sample, this convergence effect starts with significantly lower TIR values at approximately 1 <sup>×</sup> <sup>10</sup>−<sup>12</sup> Km2W−<sup>1</sup> . However, it seems that in samples with a higher thermal conductivity, the influence of the TIR also increases. Nevertheless, it must be stated that SThM measurements with different material combinations cannot be compared directly, as the TIR depends on numerous influence factors as listed above. As a result, we may conclude that a TIR at the TCA has some influence and will reduce the ideal simulated values of the TAR and the NTEF in Sections 4.4 and 4.5. This effect should be considered by other researchers assuming specific values for the TIR.

#### *4.9. Mesh Verification*

findings.

*4.9. Mesh Verification*  To demonstrate the reliability of the results, we set up two parametric simulations. They are based on the simulations in Section 4.4 (h-BN sample) and 4.5 (Si sample). For the h-BN sample, the minimum case regarding the NTEF with *temp* = 50 °C and *k* = 1 was To demonstrate the reliability of the results, we set up two parametric simulations. They are based on the simulations in Section 4.4 (h-BN sample) and Section 4.5 (Si sample). For the h-BN sample, the minimum case regarding the NTEF with *temp* = 50 ◦C and *k* = 1

simulated, whereas for the Si sample, the maximum case regarding the NTEF with *temp* =

between 0.09 nm and 10 nm to obtain information about the dependency of the NTEFs and the TARs on the mesh density and to confront them with the simulation time. Additionally, the maximum size in the circle around the TCA (see Figure 5 center) is defined to be 25 times larger than the maximum size within the TCA. Figure 11 illustrates the

was simulated, whereas for the Si sample, the maximum case regarding the NTEF with *temp* = 200 ◦C was performed. The sweep parameter in these simulations is the maximum size of a triangular mesh element within the TCA (see Figure 5 right). This parameter was swept between 0.09 nm and 10 nm to obtain information about the dependency of the NTEFs and the TARs on the mesh density and to confront them with the simulation time. Additionally, the maximum size in the circle around the TCA (see Figure 5 center) is defined to be 25 times larger than the maximum size within the TCA. Figure 11 illustrates the findings. *Nanomaterials* **2021**, *11*, 491 18 of 21

**Figure 11.** Dependency of the NTEFs and the TARs on the mesh density and confrontation with the simulation time. (**a**) Sample: h-BN (compare Section 4.4). (**b**) Sample: Si (compare Section 4.5). **Figure 11.** Dependency of the NTEFs and the TARs on the mesh density and confrontation with the simulation time. (**a**) Sample: h-BN (compare Section 4.4). (**b**) Sample: Si (compare Section 4.5).

It is obvious that with a finer mesh (left part of the *x*-axis), the values for the NTEFs and the TARs show a kind of saturation effect. For the h-BN sample, there is only little change for x values smaller than 0.2 nm, which was the chosen mesh density in all previous simulations. For the Si sample, this saturation effect starts below 0.4 nm on the *x*-axis. On the other hand, the simulation times increase enormously with x values below 0.2 nm (h-BN sample) and 0.1 nm (Si sample). The final results of the present work were obtained by performing 60 single simulations, excluding the significantly larger number of "presimulations". Therefore, the chosen meshing strategy with a maximum size of a triangular mesh element at the TCA of 0.2 nm provides a good compromise between accuracy and simulation time. It is obvious that with a finer mesh (left part of the *x*-axis), the values for the NTEFs and the TARs show a kind of saturation effect. For the h-BN sample, there is only little change for x values smaller than 0.2 nm, which was the chosen mesh density in all previous simulations. For the Si sample, this saturation effect starts below 0.4 nm on the *x*-axis. On the other hand, the simulation times increase enormously with x values below 0.2 nm (h-BN sample) and 0.1 nm (Si sample). The final results of the present work were obtained by performing 60 single simulations, excluding the significantly larger number of "presimulations". Therefore, the chosen meshing strategy with a maximum size of a triangular mesh element at the TCA of 0.2 nm provides a good compromise between accuracy and simulation time.

#### **5. Conclusions 5. Conclusions**

This work provides detailed insights into the heat transfer process during realistic SThM measurements. An SEM investigation of a used SThM probe made it possible to model an "average" used probe and calculate realistic values for the TCAs. Realistic values for thermal penetration depths, TARs, and NTEFs through the tip in contact with an h-BN or a Si sample, respectively, are provided. This allows other researchers to estimate said values for their special measurement setup and may help to interpret practical SThM measurements correctly or to explain occurring effects such as topography influences. In addition, the presented values for the TARs may help to evaluate the lateral resolution of SThM measurements as the TARs help to interpret the effect of adjacent layers. Similar to the proposal of other researchers [22], it could be shown that the influence of convection This work provides detailed insights into the heat transfer process during realistic SThM measurements. An SEM investigation of a used SThM probe made it possible to model an "average" used probe and calculate realistic values for the TCAs. Realistic values for thermal penetration depths, TARs, and NTEFs through the tip in contact with an h-BN or a Si sample, respectively, are provided. This allows other researchers to estimate said values for their special measurement setup and may help to interpret practical SThM measurements correctly or to explain occurring effects such as topography influences. In addition, the presented values for the TARs may help to evaluate the lateral resolution of SThM measurements as the TARs help to interpret the effect of adjacent layers. Similar to the proposal of other researchers [22], it could be shown that the influence of convection and radiation may be neglected in such studies.

and radiation may be neglected in such studies. From the authors' point of view, one of the most interesting findings of this study is the great impact of possible TIRs, which may not be neglected. This work may also justify higher tip velocities in practical SThM measurements as tstat is estimated below 1 µs. As a single SThM measurement can take more than 1 h, a possible increase of the tip velocities From the authors' point of view, one of the most interesting findings of this study is the great impact of possible TIRs, which may not be neglected. This work may also justify higher tip velocities in practical SThM measurements as tstat is estimated below 1 µs. As a single SThM measurement can take more than 1 h, a possible increase of the tip velocities may accelerate practical measurements without a loss of thermal accuracy.

may accelerate practical measurements without a loss of thermal accuracy. However, tstat only represents the stationary time of the heat propagation of the tip-sample contact. The

to a tip velocity of 1 µms−1 without significant differences. Said practical observation, therefore, agrees with the theoretical findings in this work and may justify tip velocities of up to 20 µms−1. Nevertheless, the values obtained in this work are only theoretical results, which could hardly be verified as almost no comparative results can be found in the

However, tstat only represents the stationary time of the heat propagation of the tip-sample contact. The time constant of the entire probe may slow down the sensing mechanism. In recent practical measurements, we tried to use higher tip velocities and compared the thermal images to a tip velocity of 1 µms−<sup>1</sup> without significant differences. Said practical observation, therefore, agrees with the theoretical findings in this work and may justify tip velocities of up to 20 µms−<sup>1</sup> . Nevertheless, the values obtained in this work are only theoretical results, which could hardly be verified as almost no comparative results can be found in the literature so far. In the future, similar practical measurements need to be performed to verify the presented values.

**Author Contributions:** Conceptualization, C.M.; methodology, C.M.; software, C.M.; validation, C.M., F.K., J.W. and G.B.; formal analysis, C.M.; investigation, C.M.; resources, G.B.; data curation, C.M.; writing—original draft preparation, C.M.; writing—review and editing, C.M., F.K., J.W. and G.B.; visualization, C.M.; supervision, G.B.; project administration, G.B.; funding acquisition, G.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by "Bayerisches Staatsministerium für Wirtschaft, Landesentwicklung und Energie StMWi", grant number DIE-2003-0004//DIE0111/02. The APC was funded by Technische Hochschule Deggendorf.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available on request from the corresponding author. The data are not publicly available due to the large amount of simulation data and the difficulty to present native modeling and simulation files in a reasonable way.

**Acknowledgments:** Many thanks to our laboratory engineer Edgar Lodermeier for the support during the SEM investigations. This work contains little parts of the Master Thesis of Christoph Metzke: Scanning Thermal Microscopy—Measurements and Finite Element Method Simulations, Deggendorf, 2019.

**Conflicts of Interest:** The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

#### **References**

