Phonon Drag Contribution to Thermopower for a Heated Metal Nanoisland on a Semiconductor Substrate

Abstract

The possible contribution of phonon drag effect to the thermoelectrically sustained potential of a heated nanoisland on a semiconductor surface was estimated in a first principal consideration. We regarded electrons and phonons as interacting particles, and the interaction cross-section was derived from the basic theory of semiconductors. The solution of the equation of motion for average electrons under the simultaneous action of phonon drag and electric field gave the distributions of phonon flux, density of charge carriers and electric potential. Dimensional suppression of thermal conductance and electron-phonon interaction were accounted for but found to be less effective than expected. The developed model predicts the formation of a layer with a high density of charge carriers that is practically independent of the concentration of dopant ions. This layer can effectively intercept the phonon flow propagating from the heated nanoisland. The resulting thermoEMF can have sufficient magnitudes to explain the low-voltage electron emission capability of nanoisland films of metals and sp²-bonded carbon, previously studied by our group. The phenomenon predicted by the model can be used in thermoelectric converters with untypical parameters or in systems for local cooling.

1. Introduction

Field electron emitters based on nanostructured materials are considered promising for the use in vacuum microelectronics devices and for several other applications [1,2]. However, some of them also pose a challenge for physicists, since the physical mechanism that allows electrons to be emitted at room temperature in a relatively weak electric field (of the order of a few V/μm [3,4,5,6]) remains unclear [7]. Nanoisland films of sp² carbon or metal on crystalline Si substrates represent a “model” type of such emitters—the one that is most convenient for experimental study and theoretical consideration. Low-field emission from carbon island films was described in [8,9]. These studies have shown that the films possessed no sharp protrusions, as well as no low-workfunction spots or large workfunction jumps with the inherent “patch fields”, any of which could explain the observed emission capability. Later [10], it was found that films of refractory metals on silicon substrates can also emit electrons in an electric field of approximately the same magnitude. The insensitivity of the effect on the film material indicates that the processes that determine the emission mechanism take place not in the islands themselves, but in adjacent regions of the substrate. In [7,10,11,12], it was suggested that the key process may be associated with a thermoelectric effect, when the thermal energy released in emission centers produces electric fields that facilitate electron emission but disappear thereafter and therefore cannot be detected in subsequent studies of the emitter surface. In [11], a quantitative estimate of the expected magnitude of the thermoelectric effect was given. Yet this estimate was not sufficiently substantiated, since it did not account for the specifics of the thermoelectric effect at the nanoscale. It is well known that a description of kinetic phenomena becomes more complicated when the size of the system approaches the mean free paths of charge carriers and phonons. The phonon mean free path L_ph in silicon at room temperature can exceed 100 nm, whereas the emitting coatings in [8,9,10] contained nanoislands with a size d of the order of 10 nm, and the nanoscale effects should certainly be accounted for.

The field electron emission is usually regarded as “cold”, and nevertheless, nanometer-sized field emission centers can operate in harsh conditions, including relatively high temperature. In our experiments [8,9,10], we often observed optical radiation from emitting spots. In low-current regimes, it had a blue color, typical for luminescence. However, when the extracted current approached the threshold of destruction of an emission center, its color changed to red and then yellow, witnessing the thermal nature of the radiation. This means that in the near-threshold regimes, the centers were heated to 1000 K or an even higher temperature. This observation agrees with literary data on energy distributions of electrons emitted from nanocarbon field-emission structures—the spectra included substantially broadened peaks [13,14]. The corresponding electron temperature was estimated in [15] as ca. 2000 K. Even while the electron temperature in these conditions can be different from the lattice temperature, and the observed broadening can partially originate from other phenomena (e.g., field penetration in the emitting structure), it is clear that the temperature of the centers of “cold” electron emission can amount to at least several hundred Kelvins (albeit, not high enough for notable thermionic emission). For nanosized centers, this also means very high temperature gradients, possibly resulting in high gradients of the thermoelectric potential, which can play a vital role in the emission mechanism.

In this paper, we consider a simple model for maybe a more relevant assessment (in comparison with the one given in our previous work [11]) of the thermoelectric potential generated near an emission center of a nanoisland film or near another heated nanometer-size spot on the surface of a semiconductor.

2. Theory

2.1. Background

In bulk semiconductors at temperatures of 300 K and above, the main contribution to the thermopower originates from the temperature dependence of the mobility and concentration of charge carriers. Another part is associated with the phonon drag—a phenomenon of transfer of mechanical momentum from atomic vibrations (phonons) to mobile charge carriers. At the macroscale, the drag component dominates only at cryogenic temperatures, when the phonon mean free path is comparable to the size of macroscopic samples [16,17]. At room temperature, and if the doping level is not very high, the probability for a phonon to scatter on other phonons, impurity atoms or defects is much greater than it is for it to scatter on electrons, which reduces the phonon drag. For point contacts of micron- and submicron-scale sizes, the existing theory predicts and experiments confirm further reduction of the drag contribution to thermopower [18,19], as the boundary scattering of phonons “at the contact aperture” is added to the other processes competing with electron-phonon interaction for the mechanic momentum carried by lattice vibrations. However, for yet smaller sizes and in situations characterized by very high temperature gradients, the tendency may change [20,21].

It has been shown [22,23] that the thermal and thermoelectric properties of thin films, nanoparticles and molecules can be controlled by adjusting the position of the Fermi level (by doping or electric field) and hence the charge carrier concentration. In this paper, we will consider the possibility of changing the carrier concentration in a semiconductor near the point contact due to the action of the thermoelectric fields themselves. If under certain conditions the concentration of free carriers can increase significantly in a region smaller than the mean free path of phonons with respect to other scattering channels, this can lead to a notable increase in the phonon drag contribution and thermopower as a whole. The size-related damping of phonon drag will also be taken into account.

We expect the considered effect to be fundamentally nonlinear and gradient-dependent. This excludes the use of the best-developed theoretical models (such as BTE [24,25]) assuming near-equilibrium local phonon and electron distributions and takes into account only the linear terms with respect to deviation from equilibrium—these theories are known to fail when relative gradients of variables are large [26,27,28]. The very general Landauer formalism [20,29] employs potential diagrams, whereas phonon drag forces are non-potential (the effect is also called “the phonon wind”): the integral of the drag force over a closed trajectory can be non-zero. Hence, the following consideration will be based on the most basic physical principles, such as energy and momentum conservation. Our approach also differs from the method used in other studies, where expected manifestations of nonlocality were estimated by evaluating higher-order derivatives for the coordinate dependence of temperature [26].

2.2. Problem Statement and Geometry

We will consider the heat flow propagating from a heated metal island into a thick semiconductor plate. The size d of the thermal contact between the media has the order of several nanometers or tens of nanometers [8,9,10]. In the semiconductor, the heat is transferred by phonons. Their mean free path L_ph satisfies the condition L_ph >> d (e.g., refs. [30,31] give the value L_ph = 210 nm for bulk Si at 300 K). Hence, their propagation is geometrically similar to radial expansion from a point. We will consider only a region of the semiconductor adjacent to the perpendicular to the interface drawn through the center of the contact area, inside a solid angle of 1 steradian. We will assume that conditions are met that allow us to regard the problem in this region as one-dimensional: the motion of charge carriers and phonons occurs in the direction from a single center, and all variables depend on a single coordinate r. The contact corresponds to the coordinate r = d (see Figure 1); it has the temperature T₁ and the electric potential U. For simplicity, the contact/semiconductor boundary is ascribed the shape of a spherical segment, with an area equal to d². The remote part of the semiconductor plate corresponds to r→∞, temperature T₂ < T₁ and zero electric potential.

The metal is ideal, with infinite thermal and electrical conductivity. The heat in the modelled region of the semiconductor is transferred by lattice vibrations only. At temperatures from 300 K and higher, more than 90% of the heat flux is transferred by phonons with wavelengths of the order of a few nanometers [25]. Thus, in a problem with the characteristic scale of the order of 10 nm, the phonons can be regarded as localized particles with certain positions (r) in space. In the case of d << L_ph, it is necessary to take into account the ballistic (at least partially) nature of heat transfer in the modeled region, which implies that we cannot use the concept of temperature, even a local one [32]. The only exceptions are boundaries: the contact itself (with temperature T₁) and the remote part of semiconductor bulk (T₂). We will consider the phonon population as having two components. One component has an equilibrium distribution corresponding to the bulk temperature T₂. Its parameters are not explicitly included in the equations, but they affect other parameters of the problem: mobility of charge carriers μ, the density of the heat flux through the contact q₀, and the mean free path for the second phonon population component L_ph. This second component is a ballistic unidirectional flow of phonons propagating from the contact into the semiconductor. At the contact, this component carries heat energy flux (or power) P₀ ≡ P(d) = q₀d², and the power decreases with the coordinate as P(r) = q(r)r² due to scattering processes. Since the region under consideration is small, we will assume that with each scattering event, the interacting “ballistic” phonon “dies” (i.e., is excluded from consideration), and its energy is subtracted from the flux.

For certainty, we take an n-type semiconductor with only one species of mobile carriers—conduction electrons with concentration n(r) and with immobile donor impurity ions with a constant concentration n_D. Despite the small size of the modeled region, we consider these quantities to be continuously distributed in space. At the equilibrium (for T₁ = T₂), we have n(r) = n_D for all positions, and no electric field exists. However, the redistribution of electron density caused by the electron-phonon interaction leads to the appearance of an electric field E(r), which produces thermoelectric current I and/or voltage U. In the proposed model, the presence of a velocity distribution of electrons will be ignored—the equations will only include the velocity v_e(r) of the “average” electron, with v_e(r) = 0 in the equilibrium. When deviating from equilibrium, the local velocity will be determined by electron mobility μ, by the resultant force (the sum of the electrical force −eE(r) and the force of phonon drag) and by diffusion of electrons against the gradient of n(r).

2.3. Boundary Condition

In [10,11], the heat power generated in the emission center was estimated as 1–2 µW. However, in other situations, it may be necessary to determine the heat flux from the applied temperature difference. Fourier’s law is not applicable at distances r < L_ph where phonon distribution has two components. For such cases, a formula was proposed that is an analogue of the formula for electromagnetic blackbody radiation (see details in [24,31,33]). For distances << L_ph (to which the temperature difference is applied), it gives the values of heat flux density independent of the distance and much lower than the values predicted by Fourier’s law; even for zero distance, the heat flux density is finite. The “surface brightness” (i.e., the heat flux density per unit solid angle) of a “black” (acoustically matching the medium) metal body with temperature T is expressed as follows:

$$I_{ph}(T)={\displaystyle \frac{{\sigma }_{ph}(T)T^4}{\pi }}.$$

(1)

Here, σ_ph(T) is the phonon analogue of the Stefan–Boltzmann constant that is temperature-independent only at T << θ, where θ is the Debye temperature for the medium. In general [24,31,33], its value is determined as follows:

$${\sigma }_{ph}(T)={\displaystyle \frac{3k_B^4}{8{\pi }^2{\hslash }^3{v}_a^2}}{\displaystyle \underset{0}{\overset{\theta /T}{\int }}\frac{x^{3}dx}{e^x-1}}.$$

(2)

Here, k_B is the Boltzmann constant, v_a is the speed of sound or, more precisely, the mode-averaged phonon velocity. For the density of the heat flux from a unit area of a metal body at temperature T₁ into a unit solid angle of a semiconducting medium at temperature T₂ (see Figure 1), we can write the following:

$$q_0={\displaystyle \frac{1}{\pi }}({\sigma }_{ph}(T_1)T_1^4-{\sigma }_{ph}(T_2)T_2^4).$$

(3)

Figure 2 shows a plot of this dependence for the Si medium (v_a = 6400 m/s [30,31]), which is very close to linear. Since the calculation was performed for a “black body”, the plotted values are probably overestimated. On the other hand, they are in sensible agreement with the macroscopic formulas. Taking the bulk thermal conductivity of silicon 149 W/m·K and a temperature difference of 100 K applied to a layer of 210 nm (the mean free path of phonons in Si [31], i.e., the minimum distance for which the macroscopic formulas are approximately correct), we obtain a heat flux density of approximately 7 × 10¹⁰ W/m² or 7 × 10⁻⁸ W/nm²—which matches the value given by the graph in Figure 2 for T₁ − T₂ = 100 K. Hence, we can use formula (3) in our calculations. However, when describing the propagation of the “radiated” phonons, we must account for nanoscale effects by the introduction of additional phonon scattering factors.

2.4. Electron-Phonon Interaction

The probability of phonon scattering on electrons is determined by the concentration of the latter n(r), which can differ from the concentration of the ionized dopant n_D. Let us represent the contribution of electrons to the inverse mean free path of phonons as a product σn(r). The coefficient σ has the dimension of area and therefore can be considered as a mean cross-section of the electron–phonon interaction. Classical theory of electron–phonon interaction (see, for instance, in [34]), can be used to estimate its value σ₀ in the limit r→∞, i.e., in equilibrium conditions, in the absence of size effects. In our previous work [11], a convenient intermediate formula can be found that refers to a partial flow of phonons with energy ћω and velocity v_a, carrying a heat flux P_ћω through a spot of area A. The probability of interaction for a single electron per unit time can be expressed as follows:

$$d{N}_{\hslash \omega -e}/dt={\displaystyle \frac{16\pi {\hslash }^3}{M{m}^{2}a{v_a}^3}}\cdot {\displaystyle \frac{P_{\hslash \omega }}{\hslash \omega }}\cdot {\displaystyle \frac{1}{A}}.$$

(4)

(Compared to [11], some notations have been changed.) Here M, a and m are the atomic mass, lattice parameter and effective mass of the electron, respectively. This relationship can be rewritten as follows:

$$d{N}_{\hslash \omega -e}/dt={\displaystyle \frac{P_{\hslash \omega }}{\hslash \omega }}\cdot {\displaystyle \frac{{\sigma }_0}{A}},$$

(5)

where

$${\sigma }_0={\displaystyle \frac{16\pi {\hslash }^3}{M{m}^{2}a{v_a}^3}}.$$

(6)

This combination of material characteristics and fundamental constants has the dimension of area and can be interpreted as the cross-section of the electron–phonon interaction. Indeed, the first factor of the right-hand side of the formula (5) is the number of phonons crossing the spot per unit of time. Multiplying it by the ratio of the cross-section σ₀ to the spot area A, we really obtain the frequency of interaction events per electron. It is noteworthy that in the Debye approximation used, the cross-section σ₀ does not depend on phonon energy and is determined only by material properties.

However, the efficiency of electron–phonon interaction is known to decrease when it is restricted in space [18,19]. In the contemplated model, this property can be described by a coordinate dependence of the parameter σ in the following form:

σ(r) = σ₀ s(r),

(7)

where s(r) is a factor accounting for the confinement suppression of electron–phonon interaction. It is natural to assume that s(∞) = 1. In the literature [18], we find that the effective degree of electron–phonon interaction for contacts of small size l decreases as l/L_e, where L_e is the mean free path of charge carriers. For our region of interest, which has the conical shape with transverse dimension l equal to the coordinate r, we can describe this factor by the simplest function with suitable asymptotes for r→∞ and r/L_e→0:

s(r) = 1 − exp(−r/L_e).

(8)

However, in the region with strong electric fields (unlike the case considered in [18]), electron motion is also restricted by the distribution of electric potential. Hence, the restricted electron free path can be re-assessed as (L_e⁻¹ + r⁻¹)⁻¹. This gives us an alternative expression for s(r):

s(r) = 1 − exp[−r(L_e⁻¹ + r⁻¹)] = 1 − exp(−1 − r/L_e).

(9)

(We do not put a factor in front of the exponential to avoid confusing the base of the natural logarithm with the electron charge.)

When choosing between possible representations (8) and (9) of the function s(r), one can also take into account the fact that in the general case, the question of the relative strength of the electron–phonon interaction in local regions cannot be considered finally resolved—opinions have been expressed in favor of both its dimensional suppression [35,36,37,38] and its enhancement [39,40].

2.5. Attenuation of Phonon Flow

We use the “narrow beam approximation” for phonons: all scattered phonons of the ballistic flux are excluded from further consideration, because the scattered phonons most probably acquire sufficient transverse (other than r) velocity components to leave the modelled region without the second scattering; thus they carry out their remaining energy and momentum. As the model includes several phonon scattering channels, Matthiessen’s rule requires summation of their contributions to the inverse travel time. However, in the Debye approximation, phonons are characterized by a single velocity v_a, which allows summing the contributions to the inverse path length (−P′(r)/P(r))—the prime here and everywhere below means the derivative with respect to the coordinate.

One of the scattering channels is the electron–phonon interaction considered above. Another one is the scattering of phonons of the ballistic flow by phonons of the equilibrium population of the medium. Only U-type (Umklapp) events should be taken into account, since normal (N-) scattering does not change the total momentum carried by the heat flow. The contribution of this channel, together with scattering by defects and impurities, can be accounted for by the value of the phonon mean free path in bulk material L_ph. Assessment of the role of U-type interactions between different phonons of ballistic flow is more complicated. Its contribution obviously depends on the intensity of the flow and therefore is a function of the coordinate. Its upper estimate can be obtained if we take the contact temperature T₁ when specifying L_ph, instead of the bulk temperature T₂. However, this method would most likely significantly overrate this channel’s role. As we know, the phonon drag is predominantly determined by the long-wave part of the acoustic phonon population, and U-processes for these phonons are less probable or completely impossible. Thus, U-interactions within the propagating flow will not be included in the model at this stage, which is best justified for temperatures << θ or for small temperature differences T₁ − T₂.

In nanoscale systems, the phonon mean free path is additionally reduced by boundary scattering. Its contribution for nanowires is estimated by the inverse wire diameter. For the considered conical region (Figure 1), the diameter grows with coordinates and is equal to r, so the boundary scattering contribution to the inverse phonon mean free path, by analogy, will be taken as 1/r. This can also be interpreted as a consequence of the mismatch between the spherical wave front shape and the concept of phonons as quanta of plane waves of atomic displacements.

The resulting law of attenuation of the ballistic phonon flux propagating into the semiconductor can be given by the following formula:

$$P^{\prime }(r)=-P(r)\left(L_{ph}^{-1}+r^{-1}+\sigma (r)n(r)\right)=-P(r)\left(\gamma (r)+\sigma (r)n(r)\right).$$

(10)

Here, γ(r) = (L_ph⁻¹ + r⁻¹) is the attenuation coefficient for all scattering processes except the electron–phonon interaction.

2.6. Phonon Drag

Being “heavy” quasiparticles, phonons carry more momentum per unit energy than electrons. According to the theory (for instance, in [41]), the average momentum received by an electron is determined only by the first stage of the phonon scattering process—the absorption of the phonon by the electron (later, the phonon is re-emitted in a random direction; hence, the average value of the momentum it carries is zero). The momentum transmitted per unit of time equals the non-potential force, which eventually produces the thermoelectric potentials and/or current. A formula for the drag force per electron immersed in a unidirectional phonon flow of P(r)/r² flux density has been derived in [11] on the basis of the standard phonon drag theory [34]:

$$F_{therm}\approx {\displaystyle \frac{16\pi {\hslash }^3}{M{m}^{2}a{v_a}^4}}{\displaystyle \frac{P(r)}{r^2}}={\displaystyle \frac{{\sigma }_0}{v_a}}{\displaystyle \frac{P(r)}{r^2}}.$$

(11)

In the Debye approximation (v_a = const(ω)), it does not depend on spectral composition of the phonon flow. In our case, we must account for possible localization effects of substituting σ₀ with σ(r) (see (6)–(9)). In addition, relation (11) was obtained for a stationary average electron, i.e., for zero macroscopic current. If the average electron has a non-zero velocity v_el, then the additional Doppler factor (1 − v_el/v_a) should be introduced in the formula—so that the average velocity of electrons accelerated by the wave of atomic displacements does not exceed the velocity of this wave:

$$F_{therm}\approx {\displaystyle \frac{\sigma (r)}{v_a}}\left(1-{\displaystyle \frac{v_{el}}{v_a}}\right){\displaystyle \frac{P(r)}{r^2}}.$$

(12)

2.7. Electric Field and Electric Current

Current I must be independent of the coordinate r and at position is proportional to the local collective velocity of conduction electrons v_el and to their concentration:

$$I=en(r)v_{el}(r)r^2.$$

(13)

Factor r² is the cross-section area of the modelled region at position r. The current in our model is mainly driven by the phonon drag force, electric field E(r) and diffusion of charge carriers:

$$I=\mu F_{therm}(r)n(r)r^2-\mu eE(r)n(r)r^2-eD{n}^{\prime }(r)r^2,$$

(14)

where μ and D are the electron mobility and diffusion coefficient, respectively.

We assume that the workfunction values for the metal and the semiconductor are equal, which means the absence of intrinsic contact potentials, charges and polarization associated with chemical heterogeneity. Potential distribution in the system, as it can be demonstrated a-posteriory, is very unbeneficial for the appearance of holes, and the holes would also be pushed away by the phonon drag. Thus, the electric field E(r) originates only from the space charge of electrons with concentration n(r), with dopant ions with uniformly distributed density n_D and with “mirror reflection” (induced charge) of these charges in the conductive contact. Here we assume that the one-dimensional radial character of potential distribution in the considered region is somehow ensured—by engineering of the material parameter outside the region and/or by selecting appropriate shapes of surrounding electrodes, as is done, for instance, in Pierce-type electron guns. In this case, the relationship between distributions of electric field and electron concentration can be obtained using the Gauss theorem in spherical coordinates with the obvious boundary condition E(∞) = 0:

$$E(r)={\displaystyle \frac{e}{\epsilon {\epsilon }_0{r}^2}}{\displaystyle \underset{r}{\overset{\infty }{\int }}(n(\xi )-n_D){\xi }^{2}d\xi }.$$

(15)

According to this formula, E(r) is determined by the total electric charge located farther from the coordinate center, at ξ > r. It would seem more natural to use the charge “inside” (at ξ < r), but this also includes the induced charge in the contact (island), which is more difficult to determine, while the condition E(∞) = 0 means that the absolute values of charges at ξ < r and at ξ > r are equal. Knowing E(r), we can calculate the thermoelectric voltage between the contact and the bulk:

$$U={\displaystyle \underset{d}{\overset{\infty }{\int }}E(r)dr}.$$

(16)

Substituting (12), (13) and (15) into (14), we come to the following relation:

$$I={\displaystyle \frac{\mu \sigma (r)}{v_a}}\left(1-{\displaystyle \frac{I}{e{v}_{a}n(r)r^2}}\right)n(r)P(r)-{\displaystyle \frac{\mu e^2}{\epsilon {\epsilon }_0}}{\displaystyle \underset{r}{\overset{\infty }{\int }}(n(\xi )-n_D){\xi }^{2}d\xi }\cdot n(r)-eD{r}^2{n}^{\prime }(r).$$

(17)

This integro-differential equation defines the relationship between the sought functions P(r) and n(r). The function σ(r) that is also present here must be explicitly pre-determined by (8) or (9), or in some other form. Current I serves as a free parameter, but it can be related to voltage U through the properties of the external circuit, for example, its electrical resistance R = U/I. In this case, the current I can also be expressed via n(r)—by means of relation (16).

Equation (10) also connects the functions n(r) and P(r). In combination with (17) and with other relations mentioned in the previous paragraph, it gives us a closed system of equations. Its solution, in principle, provides a complete (albeit significantly simplified) quantitative description of phonon drag, a nanocontact, and its contribution to the thermopower. However, the general solution of this system cannot be obtained in analytical form.

2.8. Open-Circuit Regime

Analysis of the equations obtained above is greatly simplified for the case of an open external circuit (I = 0), when U is equal to the thermoEMF. Further simplification is achieved by neglecting the diffusion term in (17)—the post-analysis has shown that its effect on the final result is insubstantial. Then, after reductions, we come to the following:

$$\sigma (r)P(r)={\displaystyle \frac{v_a{e}^2}{\epsilon {\epsilon }_0}}{\displaystyle \underset{r}{\overset{\infty }{\int }}(n(\xi )-n_D){\xi }^{2}d\xi }.$$

(18)

To get this, we divided both sides by n(r), assuming that n(r) ≠ 0. However, it should be noted here that n(r) = 0 also represents a solution, even if a trivial one.

For further numerical calculations, it is convenient to introduce dimensionless variables. Dimensionless power is defined as W(r) = P(r)/β₀, where β₀ is a combination of material parameters and fundamental constants:

$${\beta }_0={\displaystyle \frac{v_a{e}^2}{{\sigma }_0\epsilon {\epsilon }_0}}.$$

(19)

For this variable, (18) transforms to the following:

$$s(r)W(r)={\displaystyle \underset{r}{\overset{\infty }{\int }}(n(\xi )-n_D){\xi }^{2}d\xi }.$$

(20)

After derivation, it is as follows:

$$s(r)W^{\prime }(r)+s^{\prime }(r)W(r)=-(n(r)-n_D)r^2,$$

(21)

or

$$W^{\prime }(r)=-{\displaystyle \frac{(n(r)-n_D)r^2}{s(r)}}-{\displaystyle \frac{s^{\prime }(r)}{s(r)}}W(r).$$

(22)

On the other hand, formula (10) for the dimensionless power is as follows:

$$W^{\prime }(r)=-W(r)\cdot \left(\gamma (r)+\sigma (r)n(r)\right).$$

(23)

By equating the right-hand sides of (22) and (23) and expressing n(r), we obtain an algebraic relation between n(r) and W(r):

$$n(r)={\displaystyle \frac{n_D+\frac{W(r)}{r^2}\left(\gamma (r)s(r)-s^{\prime }(r)\right)}{1-\frac{{\sigma }_{0}W(r)}{r^2}{\left(s(r)\right)}^2}},$$

(24)

where γ(r) and s(r) are explicitly defined known functions. Substituting this into (23), we arrive at the differential equation for W(r):

$$W^{\prime }(r)=-W(r){\displaystyle \frac{\gamma (r)+{\sigma }_{0}s(r)\left(n_D-\frac{W(r)}{r^2}{s}^{\prime }(r)\right)}{1-\frac{{\sigma }_{0}W(r)}{r^2}{\left(s(r)\right)}^2}}$$

(25)

At the first glance, the function W(r) can be determined by solution of this equation with a boundary condition for the heat flux through the contact W(d) ≡ W₀ = q₀ d²/β₀. However, this is correct only for low q₀, below the critical value

$$q_c={\displaystyle \frac{{\beta }_0}{{\sigma }_0{\left(s(d)\right)}^2}}.$$

(26)

Otherwise (i.e., for q₀ > q_c), equations (24) and (25) cannot have physically reasonable solutions for the coordinate range d < r < r_c, where r_c is the positive root of the equation:

$$W(r_c)={\displaystyle \frac{r_c^2}{{\sigma }_0{\left(s(r_c)\right)}^2}}.$$

(27)

It is easy to show that r_c > d when q₀ > q_c; equity q₀ = q_c corresponds to r_c = d when the “critical” region d < r < r_c narrows to the contact interface.

For r > r_c, the denominator in (25) is positive, so the derivative W′(r) is negative. (For the range of interest, the numerators in (24) and (25) are positive. Otherwise, the situation would be even more complex and interesting than that described below.) Consequently, the function W(r) representing the residual power of the ballistic phonon flux monotonically decreases with the coordinate. This natural behavior is violated at r < r_c. The denominator in (25) is negative here, W′(r) is positive and the heat power function W(r) increases with distance from the contact. Meanwhile, the electron concentration takes negative values (as noted earlier, this does not mean the presence of holes)—n(r) function is described by (24), where the denominator of the right-hand side is the same as in (25). This result should be considered physically invalid. It only means that at r < r_c, the balance of forces for electrons cannot be achieved—the drag force here always exceeds the electrostatic force. At the same time, the open-circuit condition I = 0 can still be satisfied—by the solution n(r) = 0, discarded when the formula (18) was derived. It means the absence of mobile charge carriers between d and r_c. Thus, the full solution of the considered mathematical problem for electron concentration in the case of q₀ > q_c can be written as follows:

$$\left\{\begin{array}{l}n(r)=0,ford<r\le r_c\\ n(r)={\displaystyle \frac{n_D+\frac{W(r)}{r^2}\left(\gamma (r)s(r)-s^{\prime }(r)\right)}{1-\frac{{\sigma }_{0}W(r)}{r^2}{\left(s(r)\right)}^2}},forr>r_c\end{array}\right..$$

(28)

Accordingly, the solution for the residual heat flux function W(r) given by (25) is correct everywhere for q₀ ≤ q_c, but for q₀ > q_c, it is valid only at r > r_c. In the carrier-depleted region d < r < r_c, the function W(r) is described by a simpler equation, which follows from (23):

$$W^{\prime }(r)=-W(r)\gamma (r)\equiv -W(r)\cdot \left(L_{ph}^{-1}+r^{-1}\right).$$

(29)

With the boundary condition W(d) ≡ W₀ = q₀ d²/β₀, its solution is as follows:

$$W(r)={\displaystyle \frac{q_0{d}^3}{{\beta }_{o}r}}\exp \left(-{\displaystyle \frac{r-d}{L_{ph}}}\right).$$

(30)

Representing the residual heat flux, function W(r) must be continuous. Therefore, relations (27) and (30) are valid simultaneously at r = r_c:

$$W(r_c)\equiv W_c={\displaystyle \frac{q_0{d}^3}{{\beta }_o{r}_c}}\exp \left(-{\displaystyle \frac{r_c-d}{L_{ph}}}\right)={\displaystyle \frac{r_c^2}{{\sigma }_0{\left(s(r_c)\right)}^2}},$$

(31)

or

$$r_c^3=d^3{\displaystyle \frac{q_0{\sigma }_0{\left(s(r_c)\right)}^2}{{\beta }_o}}\exp \left(-{\displaystyle \frac{r_c-d}{L_{ph}}}\right).$$

(32)

This algebraic equation allows us to calculate r_c. This parameter (the outer radius of the depleted region) and the value W(r_c) ≡ W_c calculated from (31) can further be used to determine the boundary condition necessary to solve Equation (25) at the interval r ≥ r_c in the case of q₀ > q_c.

The procedure described above allows us to obtain a numerical solution for the residual heat flux function W(r). Within the framework of the model under consideration, this is sufficient to assess the phonon drag contribution to the thermoelectric voltage.

2.9. Thermoelectric Potential for Open Circuit Regime

The formulas for the potential U depend on whether the density of the incident heat flux q₀ is less or greater than the critical value q_c. In the case of q₀ ≤ q_c, the condition (20) (expressing the balance between the electric and phonon drag forces) is valid over the entire coordinate range. Comparing (20) with (15), we arrive at the following:

$$E(r)=es(r)W(r)/(\epsilon {\epsilon }_0{r}^2).$$

(33)

Therefore (see also (16)), the thermoelectric potential can be obtained by integrating the function W(r), which is the solution of Equation (25), with the boundary condition W(d) = P₀/β₀:

$$U={\displaystyle \frac{e}{\epsilon {\epsilon }_0}}{\displaystyle \underset{d}{\overset{\infty }{\int }}(s(r)W(r)/r^2)dr}.$$

(34)

In the case of q₀ > q_c, the formulas for E(r) and U are different. In the depleted layer (at r < r_c), the electron concentration is zero. Hence, for this coordinate range,

$$E(r)={\displaystyle \frac{e}{\epsilon {\epsilon }_0{r}^2}}\left[{\displaystyle \underset{r}{\overset{r_c}{\int }}(-n_D){\xi }^{2}d\xi }+{\displaystyle \underset{r_c}{\overset{\infty }{\int }}(n(\xi )-n_D){\xi }^{2}d\xi }\right]={\displaystyle \frac{e}{\epsilon {\epsilon }_0{r}^2}}\left[{\displaystyle \frac{-n_D}{3}}\left(r_c^3-r^3\right)+s(r_c)W(r_c)\right].$$

(35)

Taking W(r_c) from (31), we can re-write this in another form:

$$E(r)={\displaystyle \frac{e}{\epsilon {\epsilon }_0{r}^2}}\left({\displaystyle \frac{r_c^2}{{\sigma }_{0}s(r_c)}}-{\displaystyle \frac{n_D}{3}}\left(r_c^3-r^3\right)\right).$$

(36)

Farther from the contact (at r > r_c), the condition (20) is valid, and the electric field should be determined from (33).

Voltage U is calculated by integrating E(r) over the entire coordinate range. Hence, for the case of high heat flux intensity (q₀ > q_c), we have the following:

$$U={\displaystyle \underset{d}{\overset{\infty }{\int }}E(r)dr}={\displaystyle \underset{r_c}{\overset{\infty }{\int }}E(r)dr}+{\displaystyle \underset{d}{\overset{r_c}{\int }}E(r)dr}={\displaystyle \frac{e}{\epsilon {\epsilon }_0}}\left({\displaystyle \underset{r_c}{\overset{\infty }{\int }}\frac{s(r)W(r)}{r^2}dr}+{\displaystyle \underset{d}{\overset{r_c}{\int }}\left[\frac{r_c^2}{{\sigma }_{0}s(r_c)r^2}-\frac{n_D}{3}\left(\frac{r_c^3}{r^2}-r\right)\right]dr}\right),$$

(37)

or:

$$U={\displaystyle \frac{e}{\epsilon {\epsilon }_0}}\left({\displaystyle \underset{r_c}{\overset{\infty }{\int }}\frac{s(r)W(r)}{r^2}dr}+{\displaystyle \frac{r_c^2}{{\sigma }_{0}s(r_c)}}\left({\displaystyle \frac{1}{d}}-{\displaystyle \frac{1}{r_c}}\right)+n_D\left({\displaystyle \frac{r_c^2}{2}}-{\displaystyle \frac{d^2}{6}}-{\displaystyle \frac{r_c^3}{3d}}\right)\right).$$

(38)

In this formula, the function W(r) is the solution of differential Equation (25) with the boundary condition (27).

In the case of q₀ = q_c, the depleted layer is narrowed to zero width at r_c = d, and formulas (34) and (38) coincide.

3. Results of Numerical Simulations

The behavior of the one-dimensional model for the open-circuit regime described above was tested by numerical simulations, which were performed with the use of Wolfram Mathematica software 14.1. For these early tests, we took material parameters for Si at 300 K [11,31], ignoring their temperature dependence: L_ph = 210 nm, v_a = 6400 m/s, ε = 12, a = 0.543 nm, M = 28 a.m.u., m = 1.08 m_e. This gave as the following the parameters of the model calculated in accordance with formulas (6), (19) and (26): σ₀ = 8.8 nm², β₀ = 1.75 × 10⁻⁷ W, q_c = 2.0 × 10⁻⁸ W/nm². Other model parameters were varied. This included varying the value of dopant concentration n_D, but the results for U differed significantly (by more than a few percent) only at extremely high heat flux densities. Below, we show only the data for the lowest doping level n_D = 10¹⁵ cm⁻³ = 10⁻⁶ nm⁻³.

Figure 3 shows the results of simulations for two values of contact diameter d and three values of the incident flux density q₀. The latter is expressed through its ratio to the “critical” flux density q_c. According to the graph in Figure 1, in the “blackbody” approximation of the value q_c = 2.0 × 10⁻⁸ W/nm² corresponds to a temperature difference of ΔT ≈ 30 K; for real interfaces, this estimate will be greater. Parts (a) and (b) of the figure present the functions of electron density n(r). Curves 1 and 2 show that at small and moderate values of q₀, the excess electron density is concentrated near the contact. At q₀ > q_c (curve 3), a density maximum with infinite peak value appears, separated from the contact by an electron-depleted region. Over a ca. 10 nm-wide range near the maximum, the electron concentration is greater than 10⁻³ nm⁻³, which means that it exceeds the bulk value n = n_D = 10⁻⁶ nm⁻³ by three orders of magnitude or more.

Parts (c) and (d) in Figure 3 present distributions of the residual power of the ballistic phonon flux W(r) normalized by its boundary value W₀. The decrease in W(r) for all the plots in Figure 3c,d occurs at distances significantly shorter than L_ph = 210 nm used in the calculation, i.e., it is mainly determined by scattering on electrons and/or the contribution of the “geometric” factor 1/r. Hence, the choice of the L_ph value, within reasonable limits, would not have a substantial effect on the calculation results. Curve 1 corresponds to the minimum value q₀ = 0.1q_c, so the excess electron density (n(r) − n_D) is relatively low—its contribution to phonon flux attenuation is much less than the contribution of boundary scattering. Notably, at distances below r_c, this curve goes very close to curve 3 representing the data for the highest intensity of phonon flow q₀ = 10q_c. In the latter case, electron density at r < r_c is zero, i.e., not far from the low dopant concentration n_D = 10⁻⁶ nm⁻³. In general, the difference between curve 1 and the other curves in Figure 3c,d reflects the share of the flux absorbed by the excess electrons in each case. For d = 10 nm (part (c)), this share is less than for d = 30 nm (part (d)), due to the stronger effect of the boundary phonon scattering.

The dependences of the voltage U on the heat flux density q₀ (that is, in fact, on the contact temperature difference) are shown in Figure 4a. The calculations were performed for three contact sizes d, without taking into account the suppression of electron–phonon interaction due to confinement (we took s(r) = 1). Yet, the other dimensional effect—boundary (geometrical) scattering of the phonon flux—caused noticeable reduction in the thermoEMF. This conclusion follows from the observed downward shift of the curves with decreasing d. As expected, the plots are nonlinear—with complex behavior near q₀ = q_c, when the depletion region appears. The growth of U with increasing q₀ is generally not very fast—slower than linear. However, the voltage values predicted by the model are even higher than could be expected. This result will be discussed in the following section.

Figure 5 compares the calculation results for different forms of the function s(r) describing the confinement effect on electron–phonon interaction—see (8) and (9). The smallest contact size d = 10 nm is chosen, since for large sizes, the effect is less pronounced. The upper plot (curve 1) corresponds to L_e = 0, when the effect is totally disregarded (s(r) = 1). For the largest value of this parameter L_e = 100 nm (curve 3), we see reduction of U at low q₀ by approximately an order of magnitude. However, when we also took into account the effect of localization on the electrons themselves according to formula (9), the result (curve 4) was close to the result for L_e = 0.

In Figure 4b, the data of Figure 4a are re-plotted to show the thermoelectric voltage dependence, not on the incident flux density q₀ but on the total power at the contact P₀ = q₀d². In this form, the result can be used when considering the problem of the potential of low-field electron emission centers [8,9,10,11], where the boundary condition is set not for the emitting island temperature but the total power released there. The graph in Figure 4b shows that the highest thermoelectric voltage is achieved for the smallest island size. In more detail, the dependence of U against the island size is disclosed by the graph in Figure 6. The shown plot corresponds to P₀ = 1 μW—the total power estimates given in articles [10,11]. We expected to find a maximum in the dependency U(d) corresponding to the highest emission capacity of the islands, but the plot in Figure 6 is monotonous. Analysis shows that with the material parameters and formulas used, the maximum is located at yet lower values of d, where the model is apparently irrelevant.

4. Discussion

4.1. Physical Mechanism

The results of 1d numerical simulations revealed a physical mechanism that can affect the phonon drag contribution to the thermopower in a metal-semiconductor nanocontact. With a sufficiently high temperature difference and a large heat flux through the nanocontact, this phenomenon can prevail over the factors that decrease the drag contribution: namely, the boundary scattering of phonons and the suppression of electron–phonon interaction by confinement. This can be achieved due to a different distribution of electric charge established with the application of a large temperature difference to a nanocontact—in comparison both with a lower temperature difference and with the macroscopic case.

In the macroscopic case, the temperature gradient and the phonon drag force it produces (as well as the action of other mechanisms for creating thermopower) are distributed in the sample approximately uniformly, while the excess charge of free carriers is formed at the sample boundaries or in contact electrodes. Thus, the phonon drag acts on mobile carriers whose concentration is equal to the equilibrium value, i.e., the concentration of dopant ions. If this value is low and the medium temperature is high (non-cryogenic), the scattering of phonons on phonons and impurities is a much more probable process than their scattering on electrons, which makes the phonon drag inefficient.

The performed study showed that the situation changes in the case of nanocontact. The excess concentration of free carriers is formed at a small distance from the nanocontact (see Figure 3a,b), which can be much smaller than the phonon mean free path with respect to scattering on phonons, impurity atoms, and defects. If the temperature difference has the proper sign and sufficient magnitude, the electron concentration in this region can significantly exceed the bulk value. Then the ratio of the probabilities of phonon scattering on phonons and on charge carriers shifts towards the latter, even taking into account additional scattering at the contact boundaries. This can lead to the increase in the phonon drag contribution to thermopower and/or thermoelectric current. Within the framework of the model, at a heat flux density of the order of q_c (see formula (26)), a concentration of free carriers really achieves the value at which scattering by electrons becomes the dominant scattering process for phonons. Behavior of the U(q₀) near q₀ = q_c is complex because of the formation of a layer depleted of the carriers. The influence of this process in the thermoelectric current in regimes with finite external load can be even more dramatic.

Discussing the phonon drag near the nanocontact, we can consider possible analogies with other physical phenomena. One of them is the pressure of acoustic waves, which can be described as a result of the coherent addition of phonons emitted by an oscillating boundary. With increasing distance, the coherence is lost, and the effect transforms into ordinary thermodynamic pressure as a result of the incoherent action of phonons. The electromagnetic analogue is the ponderomotive force acting on charged particles near a metal tip to which an rf potential is applied [42]. This force is determined by the near field, which at short distances has a magnitude significantly greater than the magnitude of the radiated electromagnetic wave. It is quadratic with respect to the field magnitude and pushes charges of any sign out of the region of its concentration—just like phonon drag, which acts on electrons and holes in the same direction. Yet another analogy can be drawn to the thermionic electron emission, with different regimes of current limitation by the space charge and by the temperature.

4.2. Relevance of the Model

Numerical modeling performed on the basis of the proposed one-dimensional theoretical model gave us peculiarly high values of thermoelectric voltage. They can be compared with the predictions of classical macroscopic theory as follows. According to previous estimates, the temperature difference ΔT = 30 K between the contact and the semiconductor bulk in the blackbody approximation corresponds to a heat flux density through the contact ca. 2 × 10⁻⁸ W/nm², i.e., equal to q_c in the model under discussion. For this value of q₀, the model predicts (see Figure 4a) a phonon drag thermoEMF in the range of U = 0.5–3.9 V, which corresponds to effective Seebeck coefficients U/ΔT as high as 17–130 mV/K. This greatly, by one to two orders of magnitude, exceeds the macroscopic values of the Seebeck coefficient even for lightly doped silicon (ca. 0.4 mV/K), as well as the maximum values given by the Mott formula (up to 1.5 mV/K). However, it should be noted that at this stage, the constructed model is based on numerous simplifications, which in most cases lead to an overestimation of the final values. Thus, the presented results of numerical modeling should be considered as upper estimates, which can be useful for assessing the prospects for the implementation of the predicted phenomenon and the need for further research in this direction and improvement of the model.

One of the major simplifications consisted of ignoring the actual distribution of electron velocities and writing the equations for the average electron. In effect, this means adopting a zero value for the Debye radius, a parameter used to account for the thermal velocity spread in electrostatic problems of semiconductor and plasma physics. The n(r) distributions shown in Figure 3a,b are the narrowest of those obtained in the simulations. They correspond to the bulk concentration of ionized dopant atoms n_D = 10¹⁵ cm⁻³ = 10⁻⁶ nm⁻³. The Debye radius for this concentration and temperature of 300 K is approximately 130 nm, and widening the distribution of n(r) by this value would obviously change its shape very significantly. However, it should be noted that near the maximum, the carrier concentration here greatly exceeds the value of 10¹⁸ cm⁻³ = 10⁻³ nm⁻³. If this value is used to calculate the Debye radius, it will be only ~4 nm, and such a broadening of the n(r) profile will be less substantial. Due to this, neglecting the thermal spread of electron velocities seems justified.

The constructed model also neglects the discreteness of the electric charge, although an electron concentration of 10¹⁸ cm⁻³ corresponds to only one particle in a 10 nm cube. However, this approach is generally accepted: in quantum mechanical considerations, electrostatic potentials are often associated with a spatially continuous probability density function. A certain error may be introduced by using the blackbody phonon emission formula to describe the properties of the metal/semiconductor interface, but we see no alternatives here.

Some of the sources of errors are obvious and can be taken into account in the next stages of model development. For instance, this concerns the use of material parameters for 300 K in all calculations. The largest value of q₀/q_c = 10 used in the simulations corresponds to ΔT = 300 K or more, and material parameters (such as L_ph) should be taken as temperature- or position-dependent. In addition, the largest used contact size d = 100 nm hardly satisfies the condition d << L_ph, which must also be accounted for. Speaking more generally, the condition for applicability of the 1d model should be formulated more strictly: as d << L_Σ, where L_Σ is the mean free path of phonons with respect to all types of interactions, including the electron–phonon interaction. Obviously, this quantity is position-dependent and can be obtained only by solution of the equations for which it serves as a relevance criterion. However, the situation of the necessity for a-posteriori relevance checking is not unusual in physics. An even more serious problem concerns the use of parameters defined for equilibrium conditions in a purely non-equilibrium description.

A significant error in numeric results can be associated with the use of the approximate Formula (6) for definition of the cross-section σ₀, since the predicted effects are roughly proportional to σ₀². The formula is based on very bold consideration [11,34] where electron–phonon interaction strength was assessed, simply speaking, from upper estimates of interatomic forces. Hence, the accuracy of this assessment may be rather low. In principle, it could be improved via comparison with the existing data on phonon drag in bulk silicon [17].

The use of a one-dimensional model for description of the process is another serious simplification. In a more rigorous approach, the phonon flux should be considered anisotropic. It is natural to assume that the flux density is maximum near the normal to the contact spot and close to zero near the boundary of the semiconductor plate. Then the largest thermoelectric voltage generated in the direction of the maximum phonon flux density will be shunted by other parts of the contact region with lower flux density. This will lead to the formation of stationary eddy currents [27] with a subsequent decrease in the observed thermoelectric voltage—a similar effect was discussed in [19]. Therefore, a two- or three-dimensional model will most probably show lower values of thermoEMF and current, especially in the presence of surface electrical conductivity along the semiconductor boundary. The thermoEMF predicted by this early one-dimensional model can be realized only with the intentional or accidental fulfillment of special conditions that ensure the elimination of the shunting effect of the boundary.

Summarizing all these apparent objections, we might conclude that the model in its present state can provide only very rough estimates of quantities and maybe some insights into the physics of thermoelectricity on the nanoscale. The above-mentioned sources of error should be given special attention later when developing a more detailed model.

Yet on the other hand, the thermopower coefficient values exceeding the macroscopic ones by several orders of magnitude are in good agreement with the benchmarks given by very different theories for the “ideal” thermoelectric [43,44], for “supernodes” [45], for molecular junctions [46,47,48] and for nanostructures employing quantum interference effects [21,44]. Being much simpler, the suggested model may describe the same “ideal” conditions. In addition, the obtained results agree with the results of works [8,9,10], where high electron emission capability of discontinuous films (presumably enhanced by the thermoelectric effect) was correlated with the presence of islands with a size of ca. 10 nm. The estimates of the thermoelectric voltage obtained in this work are of the order of several volts, which may be sufficient to activate the emission by the “patch field” mechanism. In the case of island films, favorable conditions for producing high thermoelectric potentials could have been realized as a result of a rare combination of random factors, which explains the small number of emission centers observed in [8,9,10].

4.3. Prospects of Further Study and Implementation of the Effect

It should be noted that although the proposed model predicts high values of thermoelectric voltage at nanocontacts at large temperature differences, direct experimental observation of the effect may be difficult. The maximum current generated by the mechanism under consideration is small due to the low speed of the produced motion charge carriers, which is limited by the Doppler effect at a level no higher than v_a. Another limitation is related to the depletion layer appearing at high values of heat flux density. Even if the total depletion described by the present model is mitigated in regimes with non-zero current, the charge carrier concentration in this layer is still expected to be low, implying a high internal electrical resistance of the thermoelectric voltage source. Its effect may be particularly strong due to the possible influence of surface conductivity and eddy currents mentioned above. Experimental study of the effect can probably be carried out using an atomic force microscope. Comprehensive description of scanning thermal microscopy method can be found in the classic study [49], and a more recent short review [50] is devoted to local measurements of the thermoelectric coefficient. A good example of the method was given by the work [51], where the absence of size-related suppression of thermopower for point contacts down to 4 nm in size was experimentally demonstrated. Layout and preliminary feasibility analysis for our own planned experiments were disclosed in our previous paper [52].

The next step in our work with the suggested model will consist of numeric simulations of current-voltage characteristics described by Equations (14), (16) and (17) for finite resistance of the external circuit, i.e., for non-zero current regimes. The results of the present work show a drastic increase in mobile carriers’ concentration near the contact for some regimes and depletion for others, which must have an even stronger effect on electric conductivity and current than that on potentials. In further stages of our studies, the basic limitations of the presently over-simplified model will be addressed. We expect that proper regard of two-dimensional or three-dimensional character of the processes (e.g., eddy currents) can substantially reduce the voltage estimates, and the same applies to the more correct definition of the interaction cross-section σ. Both these problems can be solved by two-dimensional simulations using a finite-element software—simulations for a macroscale planar system can help in definition of σ. The equations derived in the present work can be relatively easily adopted to the 2D geometry. However, the statement of the relevant conditions at the semiconductor boundary can represent a serious problem and can strongly influence the results of the simulations.

High internal resistance and low currents of structures based on nanocontacts determine the spheres of their possible implementation. They may be similar to the ones suggested for molecular junctions or supernodes [45,46,47,48], but only if large temperature differences are available. It seems likely that in the near future, they can hardly compete with bulk thermoelectric materials for energy applications. However, there remains the possibility of its use in niche areas where unique (atypical) parameters will be demanded, such as local cooling [45,53].

5. Conclusions

In this paper, we considered the possible contribution of phonon drag to thermopower for a special situation of a metal/semiconductor nanocontact with high temperature difference applied. The study was performed to find a plausible explanation for the phenomenon of low-field cold emission of electrons by nanoisland films of metals or carbon on silicon wafers, previously investigated by the authors.

Assuming that in the region adjacent to the heated nanocontact, the phonon distribution is far from equilibrium, we used a very straightforward approach for description of electron–phonon interaction representing both phonons and mobile charge carriers as classical particles. A system of one-dimensional equations for function of phonon flux power and density of conduction-band electrons was suggested and analyzed in detail for the case of an open-circuit regime. Numeric simulations performed on the basis of the developed equations have shown that in these conditions, the electron density near the contact can drastically increase, which results in enhanced contribution of phonon drag to thermopower even at elevated (>300 K) temperatures. The resulting effective thermoelectric coefficients predicted by the suggested model are one to two orders of magnitude greater than its value for macroscopic samples.

In discussion, we agree that the suggested model in its present state may be too over-simplified to provide reliable numeric values of quantities. However, it may be useful for understanding the distinguishing features of thermoelectric phenomena at the nanoscale and can be used as a staring framework for future more advanced studies.