Numerical Simulation of Electro-Thermal Properties in FDSOI MOSFETs Down to Deep Cryogenic Temperatures

: An original reformulation of the thermopower S , heat conductivity K and heat capacitance C in bulk silicon for electrons and phonons is first proposed. Closed-form analytical approximations for these coefficients as a function of Fermi level, temperature and/or carrier concentration are developed for implementation in TCAD simulation. These analytical expressions for S , K and C have been employed to simulate the electro-thermal properties of a FDSOI (Fully depleted silicon on insulator) MOSFET versus front gate voltage from room down to very low temperature. The obtained results allow discriminating the electron and phonon contributions to the whole properties. These analyses could be very useful to further performing TCAD simulations of FDSOI MOSFETs down to very low temperature for full assessment of electro-thermal performances.


Introduction
The cryogenic microelectronics is still a crucial research subject as allowing circuit performance enhancements in terms of operation speed, turn-on behavior, thermal noise reduction, punch-through decrease [1][2][3][4][5][6][7][8][9]. It can be used in many fields of application such as high speed computing, detection and sensing, spatial electronics and lately in readout CMOS electronics and quantum-bit MOS devices for quantum computing [10,11]. In this context, there are many challenges in characterization and modelling of MOS devices down to deep cryogenic temperatures. Recently, efforts have been made for achieving TCAD device simulations down to deep cryogenic temperatures [12][13][14]. However, electro-thermal numerical simulations are still missing and are very challenging to be performed down to very low temperatures. They could bring paramount information about the electro-thermal properties such as electronic conductivity, thermopower, Peltier coefficient and thermal conductivity of MOS devices operated down to sub Kelvin temperature and further used for instance in self-heating simulation.
Therefore, in this work, we first aim at developing original reformulations of thermopower, heat conductivity and heat capacitance for electrons and phonons in bulk silicon as a function of Fermi level position, temperature and/or carrier concentration, as well as closed-form analytical approximations useful for TCAD simulation. Then, these approximations are used for the assessment by numerical simulations of electro-thermal performances of FDSOI (Fully depleted silicon on insulator) MOSFETs operated down to very low temperatures (100mK). These results will constitute a first step to further 2D/3D TCAD simulations of electrothermal properties (transport, self-heating…) in FDSOI MOSFETs at deep cryogenic temperatures.

Theoretical Background
In this section, we first reformulate in an original way the electronic parameters such as the electronic thermopower , e S the electronic heat conductivity e K and the electronic heat capacitance e C needed for the simulation of the channel properties in a FDSOI MOSFET. We also recall the phonon-related parameters such as the phonon heat capacitance , ph C the phonon heat conductivity th K and the phonon drag thermopower ph S necessary for a complete simulation of charge and heat transport. To this end, we also establish closed-form analytical approximations of the thermal and transport coefficients useful for TCAD simulation.

Electronic Parameters
The electron concentration in a bulk semiconductor is obtained from the Fermi-Dirac statistics as [15], where ( ) C N T is the effective density of states ( and In the case of degenerate (metallic) statistics (deg.), i.e. 0,  up to , f u such that the carrier concentration, the electronic thermopower, the reduced electronic heat conductivity and the reduced heat capacitance can be obtained using the Sommerfeld expansion of the integrals (8), (9) and (10) in the form [15], Note that Eq. (17) translates the so-called Wiedemann-Franz law, which gives the electronic heat conductivity in metals.
In Figure 1 are reported the typical dependence of the carrier concentration with reduced Fermi energy as given by Fermi-Dirac statistics (1) as well as its Maxwell-Boltzmann limit (11) and power law degenerate metallic approximation (15).  Figure 2 shows the corresponding variations of the electronic thermopower, the reduced electronic heat conductivity and the reduced heat capacitance with reduced Fermi energy as obtained from Fermi-Dirac statistics (8), (9) and (10), as well as their Maxwell-Boltzmann (12), (13) and (14) and degenerate (16), (17) and (18) (12) and (14) as well as of Eqs (16) and (18) indicate that   / / e C n k is equal to 3 / (2 / ) e S k q as confirmed by the plot of Figure 2. Actually, this original feature, which is not evident to find by direct inspection of Eqs (8) and ( (8), (9) and (10), as well as from Maxwell-Boltzmann (12), (13) and (14) and degenerate (16), (17) and (18)  T and / e C n versus f u or n over the whole carrier concentration range from semiconductor to metallic regimes. Based on the Maxwell-Boltzmann and degenerate asymptotic limits presented above and on usual interpolation smoothing method already employed for n, for example in cryogenic TCAD simulations [14], it is easy to construct the following analytical expressions for , n , It should be noted that these formulas well reduce to their Maxwell-Boltzmann and degenerate asymptotic limits for 0  f u or 0  f u , respectively (see Eqs (12), (13) and (14) or Eqs (16), (17) and (18)). The parameters in the analytical formulas have been optimized in order to minimize the error with respect to the exact results (see below). Figure 3a shows the variation of the carrier concentration with reduced Fermi energy as obtained from the rigorous Fermi-Dirac statistics (1) and from the analytical formula (15), both in log and linear scales. As can be seen from the figure, Eq. (19) provides a very good approximation to the exact solution with an error inferior to  5% over the whole reduced Fermi energy range varying from -20 to +20 (see Figure 3b), which is quite sufficient for TCAD simulation of real devices [14].  In Figure 4 are compared the variations of the electronic thermopower e S and of the reduced electronic heat conductivity /  e K T with reduced Fermi energy f u as calculated from the exact integrals (8) and (9) to those obtained with the analytical formulas (20) and 21). As can be seen, Eqs (20) and (21) offer very good approximations to the exact results with an overall error below  8% and  12%, respectively, over the whole Fermi energy range (see Figure 4b) for both the thermopower e S and reduced electronic heat conductivity T with a better accuracy for e S .
In Figure 5 are reported the variations of the reduced electronic heat capacitance / e C n with reduced Fermi energy f u as computed from the exact integral (10) along with those obtained with the analytical formula (22). As can also be seen, Eq. (22) gives a very good approximation to the exact result with an overall error below 8% over the whole Fermi energy range as for ( ) It should also be mentioned that the electronic heat capacitance of bulk silicon has been extracted from low temperature measurements of thermal capacitance by Kobayashi et al [20] for various doping levels and can be well fitted by Eqs (10) or (22) without any parameter adjustment as can be seen from Figure 6, thereby validating our reformulation of e C . (19), (20), (21) and (22) do constitute simple and analytical closed-form expressions for calculating the carrier concentration, the electronic thermopower, the reduced electronic heat conductivity and the reduced electronic heat capacitance in silicon as a function of reduced Fermi level from semiconductor regime (Maxwell-Boltzmann statistics) up to metallic regime (degenerate statistics) with a good accuracy for practical use in device numerical simulation.

Phononic-Related Parameters
As for e C , the phonon heat capacitance in a crystal can be evaluated from / ,    ph ph C U T ph U being the lattice vibrations (phonon) energy and can be equated to as [19], is the Bose-Einstein distribution function and 2 ( )  g E E the phonon density of states at low energies (i.e. for the acoustic branches). This readily yields the so-called Debye formula [19], where N is the atom density (=5×10 22 /cm 3 for silicon) and d T is the Debye temperature ( 645K  d T for silicon).
At high temperature (  d T T ), the phonon heat capacitance tends to a constant 3 , corresponds to the Dulong and Petit heat capacitance limit. At low temperature (  d T T ), the integral in Eq. (24) tends to 12.π 4 /45, such that the phonon heat capacitance ph C reduces to the usual 3 T power law [19], 3 4 12 As it was done for the electronic heat capacitance, a closed-form analytical expression can be developed for ( ) ph C T based on the two asymptotic limits at low and high temperature, yielding the simple expression, (26) provides a very good approximation to the exact formula with an error less than  12% over the whole temperature range up d T .
Figure 7b compares the variation with temperature of the phonon heat capacitance as obtained from the Debye model (Eqs 24 or 26) to bulk silicon experimental data from Ref. [21]. It indicates that the Debye model provides a sufficiently good description of the phonon heat capacitance over a wide range of temperature needed for device simulation. The lattice thermal conduction is governed by the heat transport by phonon diffusion due to a temperature gradient such as the heat flux  [22,23], ph s v is the phonon relaxation time, which is in general a function of temperature and phonon energy. If we consider, in first approximation, that  ph is independent of energy, we can factorize the product 2  s ph v outside the integral, such that the phonon heat conductivity can be related to the phonon heat capacitance according to the kinetic theory as [24], Figure 8 shows how Eq. (28) can be used to fit reasonably well the bulk silicon data from Ref. [23] for ( ), ph K T provided that the phonon relaxation time is made explicitly temperature dependent as However, in section 3, a more accurate empirical combined power law vs T will be used to fit the ( ) ph K T silicon data for device simulation purpose (see Eq (42)). Another important phenomenon appearing at low temperature is the so-called phonon drag effect, which enhances the electronic thermopower due to the interaction between diffusing phonons and electrons. The phonon drag effect has been addressed both experimentally and theoretically in the 50's by Geballe and Hull [25] and Herring [26], whereas further developments have been proposed later by Cantrell and Butcher [27] and Mahan [28]. In first instance, the Herring formula should be used to depict the phonon drag thermopower ph S and can be formulated using the Matthiessen rule as [26], where / ( )     eff q n is the effective mobility and ( )  effph T is the phonon-limited electron mobility component described by 1.5 ( ) 1400 ( / 300)     effph T T for bulk silicon [29]. As can be seen from Figure 9, Eq. (29) well reproduces the modelling results of Mahan [28], which have been calibrated on Geballe and Hull silicon experimental data [25] after having adjusted the temperature dependence of the phonon drag relaxation time such as Note that this phonon relaxation time is 2 orders of magnitude lower than the  ph entering the phonon heat conductivity formula (28). This discrepancy has been discussed in the literature and could be attributed to the difference of nature between pure phonon transport and electron-phonon interaction processes [26][27][28]. However, it should be mentioned that it is close to the one used in Mahan et al [28] and has the same trend versus T , indicating that Eq. (29) is physically sound. In any case, this is not so important since it is calibrated on both experimental and theoretical silicon data and therefore can be used in first approximation to evaluate the phonon drag thermopower in bulk silicon. Nevertheless, Eq. (29) giving the phonon drag thermopower can be considered as too simplistic to be applied at low temperature where degeneracy occurs. This is why further improvement to (29) has been proposed by Cantrell and Butcher [27] to account for the electron energy dispersion in 2D inversion layers [30]. Applying their derivation to 3D electron gas leads to the generalized expression, where the energy mobility functions ( , )  E T and ( , )  ph E T are given by, In Eq. (31), ( , )  ph E T is the phonon limited mobility function given by [ Note that the parameters used in Eq. (32) have been properly calibrated such that the effective mobility / ( )     eff q n calculated with Eqs (1), (3), (31) and (32) fits reasonably well the standard electron mobility data versus doping level at room temperature for silicon [31] as shown in Figure 10. Besides, it is worthwhile to examine which are the asymptotic limits of the Cantrell and Butcher equation (30) for Maxwell-Boltzmann and degenerate metallic statistics. In the case of MB statistics, Eq (30) can be well approximated by, which is nothing else than Herring's formula (29).
In the case of degenerate statistics, Eq. (30) reduces to, indicating that ph S now depends directly on the phonon limited mobility function at Fermi energy ( , )  ph f E T . Figure 11a shows the variations of the phonon drag thermopower ph S with temperature as obtained from Cantrell and Butcher model [27] of Eq. (30) and corresponding MB and degenerate approximations given by Eqs (33) and (34) These results also reveal that, for low doped silicon used in FDSOI MOS channel, the Herring formula of Eq. (29) calibrated on Mahan modelling results and Geballe and Hull experimental silicon data should constitute a first approximation for the phonon drag thermopower evaluation. Nevertheless, it should be mentioned that it is not able to describe the temperature and 2D carrier density dependence obtained at very low temperature in MOS inversion layers by Gallagher et al [32]. This issue goes beyond the scope of the present theoretical analysis and should be further discussed.  (33) and (34), and, (b) as obtained from Eq. (30) (solid line) and extrapolated Mahan's data [28] (symbols) for low doped silicon (Parameters: vs = 9×10 5 cm/s, τph0 = 18ns and T0= 58K).

Results and Discussion
The electro-thermal properties of a FDSOI MOS structure have been simulated after solving the Poisson equation in two dimensions using the finite element partial differential equation software FlexPDE [33]. The simulated SOI structure is illustrated in Figure 12 and features a 2nm EOT (equivalent oxide thickness) top oxide thickness, a 10nm undoped silicon film and a 20nm buried oxide (BOX). No source and drain electrodes were considered here since limiting our analysis to 1D profiles, sufficient to evaluate in a long channel the sheet quantities integrated over the silicon thickness. For simplicity, the front (resp. back) gate voltage 1 was directly applied to the top (resp. bottom) oxide external boundary. In this case, the electrical potential V across the structure follows the Poisson equation as, where  ox and  si being the oxide and silicon permittivities, respectively, and 0 V (0.55V) the mid-gap potential.
The carrier concentrations in Eq. (35a) were computed using the analytical approximation of Eq. (19). Quantum confinement effects were taken into account using Hansch's quantum correction at each Si-SiO2 interface with a quantum length of 1nm as in [14]. Once the electron density ( ) n x was known from Poisson's equation solution, the electronic thermopower e S , the reduced electronic heat conductivity / ( )  e K T and the reduced electronic heat capacitance / e C n across the silicon film were then calculated using Eqs (20), (21) and (22).
In Figure 13   The sheet carrier density s N was then calculated by integrating the electron concentration ( ) n x across the silicon channel as, 0 ( )   si t s N n x dx (36) Typical variations of the inversion charge density s N with front gate voltage 1 g V are shown in Figure 15 for various temperatures both in log and linear scales. One should clearly notice the usual strong increase of the subthreshold slope as the temperature is lowered due to the onset of Maxwell-Boltzmann statistics [1,2]. In contrast, at strong inversion above threshold, as is well known [1,2], the sheet carrier density s N becomes nearly temperature independent. The local electrical conductivity ( )  x was then evaluated by considering a standard mobility law both depending on temperature and local vertical electric field x F as [14,34], where c F is a critical field (1MV/cm) allowing to emulate the mobility degradation at high gate voltage due to the surface roughness scattering. The low field empirical mobility law 0 ( )  T accounts for the mobility increase with temperature reduction governed by phonon scattering and for mobility saturation at low temperature due to prevailing neutral defect scattering [2,35]. Knowing the local electrical conductivity ( )  x , it is then possible to calculate the sheet electrical conductivity  s by integrating the parallel contribution to charge transport of each channel layer as, (38) In Figure 16 are reported typical evolutions of the sheet electrical conductivity  s with front gate voltage 1 g V for various temperatures from 300K down to 1K, as obtained from Eqs (37) and (38). As for the 1 ( ) s g N V characteristics, one should also note the strong increase of the subthreshold slope with temperature lowering when Maxwell-Boltzmann's statistics prevails. However, in this case, the sheet electrical conductivity  s is significantly improved above threshold at low temperature due to the low field mobility enhancement. The onset of a zero temperature coefficient point (ZTC) in the 1 ( )  s g V is noticeable and is the consequence of the 1/T mobility dependence. It should be pointed out that this is in full agreement with usual cryogenic transfer characteristics data obtained in long channel Si MOSFETs [2,35]. It should also be mentioned that the effective mobility deduced from the sheet conductivity / ( )    eff s s qN exhibits variation with temperature from 400 cm 2 /Vs at 300K up to 2500 cm 2 /Vs at 10K (not shown here), which are in good agreement with the experimental data obtained on FDSOI MOSFETs [35,36]. Once having evaluated the local and sheet electrical conductivity, it is now possible to calculate the sheet electronic thermopower es S , the sheet electronic heat conductivity es K and the sheet electronic heat capacitance es C by integrating the parallel contribution of each channel layer to electronic and heat transport as [18,37], It should be mentioned that the thermal leakage in the front and BOX oxides surrounding the channel have been overlooked since their heat conductivity is at least two orders of magnitude lower than those of bulk silicon.
In Figure 18a are compared the variations with temperature of the total sheet thermal conductivity tots K with the electronic ( es K ) and phonon ( phs K ) components as obtained using Eqs (40) [39]. Figure 19 shows the variations with temperature of the total channel thermal response time  tot along with its phonon 2 ( / (6 ))   . As can be seen from the figure, the thermal response time is dominated by the phonon contribution above 10K, while it is mainly controlled by the electronic contribution below 1K. Note that for such a gate length, the channel thermal response time lies in the picosecond range below 10K.  The variations of the Peltier coefficient Pel and thermoelectric figure-of-merit factor zT with front gate voltage 1 g V can be evaluated as shown in Figure 20. As can be seen, the Peltier coefficient variations with gate voltage mainly mimic those of the thermopower above threshold, whereas there is a relative temperature independency below threshold. Similarly, the thermoelectric figure-of-merit factor variations with gate voltage are mainly controlled by those of the electrical conductivity below threshold, whereas they are mostly reflecting those of 2   s es S in strong inversion, since ths K is constant with 1 g V as being mostly dominated by the sheet phonon heat conductivity term (see Figure 18). In any case, both the Peltier coefficient and the thermoelectric figure-of-merit factor take values well below one, indicating a poor thermoelectric efficiency for such a FDSOI MOSFET, especially at cryogenic temperatures, indicating that such structures are not optimized for pure thermoelectric applications. Finally, in Figure 21 are compared, for completeness, the variations with temperature of the sheet electronic thermopower es S for two front gate voltages 1 g V and of the phonon drag thermopower ph S as given by the extrapolated Mahan model [28] of Eq. (29) calibrated on Geballe and Hull silicon data [25]. It appears from the figure that the phonon drag thermopower component should dominate the electronic thermopower component, especially below T50K (see Figure 11). However, these results have to be taken with caution since the simple model of Eq. (29) does not account for the 1 / n and 3 T dependence of the thermopower observed at very low temperature e.g. by Gallagher et al in Si MOS inversion layers [32]. Therefore, a more accurate modelling of the phonon drag thermopower component should be developed, especially at very low temperatures, as a function of carrier concentration (i.e. also Fermi level) and temperature to be further included in FDSOI MOSFET simulation, bringing more reliable electro-thermal simulation results. But this issue goes beyond the scope of this paper. Variations with temperature T of the sheet electronic thermopower Ses for two front gate voltages Vg1 (red solid lines) and of the phonon drag thermopower Sph as given by Mahan model [28] of Eq. (29) (blue dashed line) as obtained in a FDSOI structure (tox1 = 2nm, tox2 = 20nm, tsi = 10nm, Vg2=0).

Summary and Conclusions
In this work, we have proposed an original reformulation of the thermopower S, heat conductivity K and heat capacitance C coefficients for electron in bulk silicon as a function of Fermi level position and temperature applicable from Maxwell-Boltzmann to degenerate statistics. Similarly, the phonon related coefficients have also been reformulated in this context. Besides, we have also developed closed-form analytical approximations for these coefficients versus Fermi level, temperature and/or carrier concentration, and which are useful for implementation in TCAD simulation. Then, these analytical expressions for S, K and C have been further employed to simulate the electro-thermal properties of a FDSOI MOSFET versus front gate voltage from room down to very low temperature (0.1K). The obtained results enable the discrimination between the electron and phonon contributions to the whole properties. These analyses pave the way to further 2D/3D fully coupled electro-thermal numerical simulations of FDSOI MOSFETs down to very deep cryogenic temperatures for full assessment of their electro-thermal performances (transport, self-heating…).