Dr Michael Cassé

Bruna Cardozo Paz

هللا بسم والصالة والسالم أ على المرسلين كرم

I would like to thank the following people, without whom I would not have been able to complete this research, and without whom I would not have made it through. Firstly, I would like to express my deepest appreciation to my jury committee, for having the kindness of reviewing our work and attend my defense. Your

Table of contents

Conventionally, quantum computers consist of two parts: a quantum processor that comprises a set of qubits, and a classical electronic interface part required to perform the control and readout of quantum states. Factually, in order to preserve the compactness and reliability of the quantum computer and to ensure a reduced signal latency, the qubits were implemented in the cryogenic chamber along with the classical electronic interface so as to ensure the higher accuracy and the lower noise of the signals provided by these control/readout electronics. Owing to the back biasing feature, FDSOI transistors can regulate the threshold voltage, require low voltage supply, making them consequently able to exhibit low power consumption, which make them a good candidate for cryogenic operation.

The designing of reliable and optimized circuits for deep cryogenic operation requires suitable compact models to be embedded in the process design kit packages. Several efforts including this thesis are going on currently in order to provide such models. Nonetheless, the realization of standard compact models that are mature enough still requires longer time and further endeavors.

To tackle such task two approaches have been adopted by the research community. The first approach consists of adapting existing standard models for deep cryogenic operation. Such adaptation is made using empirical solutions in order to improve their accuracy and predictability in cryogenic operation. Such approach comes along with a lot of advantages, since these standard models contain already all the additional effects and features and are numerically robust. However, these advantages come along with the limitations that these standard models are not adapted for cryogenic operation, namely because they consider a 3D gas of electrons and use Maxwell-Boltzmann (MB) statistics to describe their distribution. The second approach consists of building physics-based models aimed for cryogenic operation from scratch. Such choice allows to overpass the limitations of the first approach, for instance, 2-D electron gas along with the use of Fermi-Dirac (FD) statistics can be considered initially. Nevertheless, these models are not mature enough to be implemented in Process Design Kits (PDK).

Expressly, the key stone into developing mature and suitable compact models is to understand the underlying physics that rules the MOS transistor behavior at cryogenic temperatures. For this end, in the frame of the present thesis we explore the different physical effects that manifest in transistors operating at cryogenic temperature such as the statistics that describe the electron distribution at these conditions, the subthreshold slope saturation at low temperature, and the mobility law evolution.

In our study, three levels of modeling are performed. The first one consists of the self-consistent solution of Poisson and Schrödinger (PS) equations. Such PS simulations are a useful tool for understanding the physics that governs FDSOI transistors down to deep cryogenic temperatures. The second level of modeling is presented by the numerical model which meant to be a solid background for the development of the analytical model. The third level of modeling is of an analytical nature. Expressly, the PS simulations certify the charge and current solutions predicted by the numerical model, which in turn certifies the same predicted solutions by the analytical compact model.

Correspondingly, in the introduction chapter, i.e. Chapter one, we will be exposing the frame from which emerges the challenge that we are trying to solve throughout this thesis. Accordingly, the components of a quantum computer are exposed along with the reasons the electronic interface needs to be located in the cryogenic chamber. Then, we will be discussing the justifications for the choice of CMOS technology and precisely the FDSOI transistors as the ultimate candidate to be implemented in the classical interface for deep cryogenic temperature operation. In this context, the urgent demand for compact models describing the operation of FDSOI devices at cryogenic temperatures emerges as it is a crucial element for the process-design kits to assemble reliable and optimized circuits. To produce such compact models one can seek two approaches, the first one starts from existing standard compact models, which are originally built for room temperature operation and uses empirical formulas in order to adapt them to cryogenic temperature operation. The second approach considers building fully-physics based compact models that are dedicated for cryogenic device operation. Indeed, in the frame of the present thesis we pursue the second approach.

In the course of the Chapter that uncovers the underlying physics that reigns the device behavior at deep cryogenic operation i.e. Chapter two of the thesis, we pave the way to justify some principal assumptions and choices that have to be made subsequently. Chiefly, the Maxwell-Boltzmann (MB) approximation validity down to cryogenic temperatures is discussed by the exposition of the reasons for which it does not hold in our present study and the reasons why it can be retained in some specific situations. Such discussion will be followed by the demonstration of the exponential band tails exhibited by the two-dimensional density of states. Ensuing, we expose two approaches to compute the subthreshold slope saturation manifested at deep cryogenic temperatures. In the same frame, we propose a description of the conductivity function through the use of the Kubo-Greenwood integral along with the diffusivity function in the degenerate statistics regime. Such exposition will be accompanied by the exhibition of the bell-shape mobility law that rules the electronic mobility at cryogenic temperatures. At the end of this chapter, we present an appealing effect that has been observed experimentally for backbiased FDSOI transistors and propose a legitimate explanation.

In the third chapter, dedicated to the exposition of the performed Poisson-Schrödinger simulations, we display an in-situ analysis of the band diagrams along with the subbands population mechanism appropriate for FDSOI structure operating at deep cryogenic temperatures. Such simulations were curiously performed as well at the 𝑇 → 0𝐾 limit. Correspondingly, we expose the electrostatic parameter curves namely the gate-to-channel capacitance curves that exhibit an appealing two-plateaus behavior for back-biased FDSOI structures. At the end of this chapter, the simulated gate-to-channel capacitance curves are confronted to collected C-V measurements in order to be validated.

The fourth chapter is dedicated to the development of a numerical model that is aimed to the keystone for the compact model development. Such numerical model would be based on a system of two coupled charge equations. The establishment of such system would be accompanied by certain choices concerning the charge coupling term and the quantum shift function. In this frame, we propose an extended form of the quantum shift function that is suitable for different geometrical configurations and for back-biased structures. At the end of this chapter, mathematical expressions that describe the numerical integrals involved for the computation of the drain current in the case of Fermi-Dirac statistics are presented.

In the course of the Fifth chapter, dedicated to the development of the required compact model, we range from the presented system of coupled equations in order to get an exact analytical solution for the surface potentials. Such solution is derived through a step-by-step technique that considers the application of a number of error correction steps. Similarly, closed-form analytical expressions were demonstrated for the diffusion current computation directly, however, concerning the drift current computation we presented a two-slope inversion charge linearization technique applicable for back-bias structures that considers the computation of the respective slopes at the source and saturation ends. Finally, a few short channel effects were implemented to the core model such as the velocity saturation, the DIBL and charge sharing phenomena, and the parasitic resistance effects.

Chapter 1: Introduction

In his document published in 1982 [1], Richard Feynman proposed the idea of a computer that "will do exactly the same as nature", i.e. a new kind of computer that imitates the physical laws of quantum mechanics like superposition or entanglement, i.e. a "quantum computer". Quantum computers can solve real world NP-complete problems 1 proficiently, such as efficient search in extremely large datasets, factorization of large integers in their prime factors, simulations of quantum systems for the optimization of drug synthesis, materials and industrial chemical processes [2].

In a quantum computer, standard logic bits are replaced by quantum bits (qubits), whose states can be represented as a point on the surface of a three-dimensional sphere, the so-called Bloch sphere. In this concept, standard logic '1' and '0' are replaced by quantum states |0> and |1>, and are manipulated, so as to exploit the fundamental phenomena of quantum mechanics for computation. Qubits can exist in a superposition of both states |0> and |1> simultaneously, which results in a computing power that doubles with every additional qubit, thus resulting in a massive speedup with respect to traditional computers. For example, it has been estimated that the state of a 50-qubit system, which corresponds to about one petabits, cannot be stored in the memory of the world's most powerful computers today [2].

In its fundamental core, a quantum computer comprises:

• A quantum processor, which consists of a set of qubits.

• A classical electronic interface required to perform the control and readout of quantum states.

Given that each qubit technology/implementation has its own strengths and weaknesses, the discussion is still ongoing about the choice of a standard qubit to build a large-scale quantum computer. Due to their similar nanofabrication techniques to those used in the microelectronics industry, silicon spin qubits and superconducting qubits or the so-called "Solid-state qubits" are probably the best choices for scalability (see Figure 1). The advantages of the superconducting qubits are their simple fabrication and easy control; but size-wise they are very large compared to silicon qubits, the formers are of the order of micrometers whereas the latter of the order of nanometers; they need lower operating temperature, 10mK for superconducting qubits compared to 100mK for silicon qubits. On the other hand, silicon qubits have finer compatibility to microelectronics industry than superconducting qubits, but they are more susceptible to defects, disorders or strains, invoking a long-lasting endeavor to determine the right operating voltages [3].

Figure 1. Solid-state qubits: to the left a spin qubit, to the right a superconductive qubit [4].

In the framework of the MOS-Quito project, CEA-Leti has developed its own qubit; a siliconon-insulator nanowire field-effect transistor with two gates G1 and G2 (see Figure 2). The superposition and entanglement of the electron spins under the gates set the quantum information, spin-resonance techniques are used to control or "rotate" the spins by sending nanoseconds voltage pulses and GHz microwave bursts to the gates G1 and G2. RF reflectometry is commonly used to read-out the spin states of the qubit by connecting an LCnetwork to the gate [5].

Quantum states can be very easily disrupted by the heat generated vibrations, which must be eliminated for the quantum states to manifest. Thus, the quantum processor must be cooled down to deep cryogenic temperatures, typically between 10 -100 mK.

The lifetime of a quantum state, so-called the coherence time, only lasts for a very short time, normally in the order of nanoseconds or microseconds in best cases [6]. Such periods are not long enough to execute any practical computation. Hence, they need to be maintained longer, a task supervised by the classical electronic interface. Considering the qubit sensitivity, the classical electronic interface needs to grant the following extremely challenging specifications [2]:

• Accuracy: to ensure the optimum operation of the quantum processors, the electronic signals feeding the qubits must be highly accurate in terms of amplitude, timing, frequency and phase. • Noise: the electronic signals feeding the qubits must transmit very limited noise, to guarantee the non-alteration of the quantum states. • Bandwidth: Controlling solid-state qubits is done through the generation of microwave bursts ranging from few GHz to tens GHz, in addition to the current and voltage pulses of tens MHz.

Moreover, to transcend the short coherence time, information redundancy has been proposed through the error correction techniques, by implementing the quantum information in a large number of qubits, i.e. trading off the simplicity of the system by its reliability, making the threshold of number of qubits required to perform practical computations even higher [2].

In essence, the classical interface must take care of two functions, the execution of the quantum algorithm and the fidelity of the computation beyond coherence time. High fidelity conveys that the quantum-controller needs to bring back the qubits to its initial state with a probability of 99.99% [4]. In addition, the latency of the error-correction loop must be lower than the qubit coherence time.

A generic control and readout platform is composed of several subsystems, a multiplexing and demultiplexing amplification matrix or the switch matrix, an I-to-V converter, a low pass filter, analog to digital conversion ADC components and digital to analog conversion DAC components [2], [8]. The switch matrix is a key component that directs a particular waveform to a particular qubit based on its digital address, it is placed close to the qubits, to avoid latency and synchronization problems if it propagates on a path comparable to its wavelength [6].

At first, the qubits were implemented in a Helium dilution chamber, while the control/readout electronics were kept at a distance at room temperature. The connection between the two is realized through a number of physical wires, and each qubit is controlled solely. This might be practical for simpler prototype systems where the number of qubits is below 100 [4], and compactness is not required. The limitations of this implementation are basically due to the thermal load of the large number of cables, and the latency of the error-correction loop. The result would be a big, expensive, unpractical quantum computer with low reliability [2]. Thus, such an approach cannot be maintained for a large-scale quantum computer with big density of qubits.

Scaling up the number of qubits requires a new approach, one that can ensure a compact design, while satisfying all thermal requirements, and improving reliability and debuggability of the overall system [4]. Such needs are attained if the classical logic interfaces are located in the cryogenic chamber, which would be accompanied with the advantages of an enhanced clock speed, an improved noise performance, a reduced signal latency/timing errors, and larger bandwidth [6]. In such an installation (see Figure 3), the long wiring can be removed from the system and replaced by interconnects, and the control and readout of multi-qubits will be performed simultaneously in proximity [8], resulting in a more compact and more reliable system [2]. Nonetheless, this approach will add another constraint to the classical interface performance, the power dissipation constraint. The thermal budget due to the power dissipation has to be within the thermal absorption limits of the refrigeration system. Taking advantage of the different distribution of thermal budget across the system, the control architecture can be placed at a certain stage of temperature [4]. A feasible control platform power dissipation of 30 mW, making it operating at liquid helium temperature (LHT) i.e. 4.2K [6].

Based on what has been said previously, the choice of technology for constructing this layer of classical control is largely dictated by its functionality at very low temperatures, its high speed and compactness. The electronic devices that have shown functionality at cryogenic temperatures are JFETs, HEMTs, superconducting devices based on Josephson junction, compound semiconductors, and CMOS transistors. Regarding the low power consumption, integration of billions of transistors on a single chip, and the very mature industrial technology, CMOS is the chosen candidate [2], [4], [6], [8].

The quest for understanding the underlying device's physics at cryogenic temperatures has started decades ago [9]. Some physical phenomena that are dominant at room temperature become non-dominant at cryogenic temperature, and equivalently, some physical phenomena that were of minor importance at room temperature become dominant at cryogenic temperatures, those physical phenomena will be discussed in detail in Chapter 3.

Bulk CMOS operating at cryogenic temperatures is characterized by an increase in the threshold voltage due to the carrier freeze-out, and by some kink effects that are caused by the dopant freeze-out. Nonetheless, those kink effects are either reduced or absent in a more mature and refined CMOS technologies such as the FDSOI technology. In comparison to bulk CMOS, FDSOI transistors can perform an in-situ control of the threshold voltage due to its back biasing, a low power consumption due to its low voltage supply, and a low variability due to the undoped silicon channel [10], which could offer the optimum cryogenic device performance [11].

Compact models describing the operation of MOS devices at cryogenic temperatures are crucial for the designing of reliable and optimized circuits. Nowadays, process-design kits lack models for all devices at deep cryogenic temperatures, whether it is a MOSFET, a qubit, or a passive device. The task of building compact models for cryogenic operation needs to be tackled urgently, mainly because today's circuit simulations fail to function as expected down to cryogenic temperatures. Consequently, excessive simulations are needed to be done to account for the changes in different parameters.

There are two approaches to deal with this task, the first one is to take existing standard compact models, which are originally built for room temperature operation and try to adapt them to cryogenic temperature operation through empirical formulas. This approach does not require a long time but as the underlying physics at cryogenic temperatures is not the same as at room temperature, the result would be non-physical compact models that would be very limited. Such approach can be adopted when the resources (whether human or budget or time) are limited. It should be noted that for this approach the circuit performance is not guaranteed, and they are designed with a non-negligible degree of uncertainty.

The second approach, which is a research approach and the one we choose in the frame of this thesis, is to build fully-physics based compact models that are dedicated for cryogenic device operation. Such an approach is more time consuming and the process of developing appropriate compact models is progressive, starts with a core model that predicts the long channel transistor behavior, followed by continuous add-ons/improvements for different effects such as small transistors effects, access resistances, self-heating effects, to name a few. Compared to the first approach, this approach should be less risky and the output design-wise is much more appropriate and precise.

State-of-the-art of standard compact models at cryogenic temperatures:

Current standard compact models can scale down to liquid Nitrogen temperature 77K, but at liquid Helium temperature 4.2K and below some discontinuities start appearing in the moderate inversion region. This is due to two things basically, the lacking/incorrect modeling of the device physics at this range of temperatures (as it was not initially planned for this aim), and the numerical issues generated by the very large arguments of the exponential functions used in the Maxwell-Boltzmann statistics, leading to their explosion and the model's crashing.

Moreover, at cryogenic temperatures, the Maxwell-Boltzmann approximation is not valid anymore; the intrinsic carrier concentration 𝑛 𝑖 becomes extremely small (they are found to be lying outside the range of the IEEE double precision 10 -308 → 10 308 ; for example at 4.2K, 𝑛 𝑖 = 10 -678 𝑐𝑚 -3 ), resulting in enormous arithmetic underflows in the implemented analytical expressions [12].

Efforts to models the operation of CMOS transistors started prior to the need for quantum cryogenic controllers [13], [14]. Figure 4 for example exposes the simulations obtained using the BSIM model versus the experimental data measured at 100K of an NMOSFET. This deficit of modelling the transistor's behavior at this temperature is due to the inappropriate projection of the built-in temperature dependences of the BSIM model. Such projections are not valid down to this range of temperatures [14]. Subsequent to the emergence of a quantum controller need, modelling efforts were held predominately in the four axes; each one will be discussed and highlighted hereafter.

❑ The EKV MOSFET model:

The EKV MOSFET is a simplified charge-based model that can be used for advanced CMOS technologies, and require four parameters to fit the 𝐼 𝑑 -𝑉 𝑔 transfer characteristics: the slope factor 𝑛, the specific current-per-square parameter 𝐼 𝑠𝑝𝑒𝑐 , the threshold voltage parameter 𝑉 𝑇0 , and the velocity saturation parameter 𝐿 𝑠𝑎𝑡 [15]. In such a model, the concept of inversion coefficient IC, a parameter that quantifies the channel's inversion, is introduced to replace the overdrive voltage. Both the normalized transconductance efficiency 𝑔 𝑚 𝐼 𝑑 ⁄ and the normalized output conductance 𝑔 𝑑𝑠 𝐼 𝑑 ⁄ are found to be dependent of the inversion coefficient [15]. This IC-based methodology has been proved valid for back-biased FDSOI transistors at room temperature [16].

At cryogenic temperatures on the other hand, the major increase of the slope factor parameter n at 4.2K is attributed to the interface-trapping phenomenon, followed by a small adjustment to introduce the change in the subthreshold slope induced by the back-gate biasing. The threshold voltage shift due to the incomplete ionization, a phenomenon that characterize doped CMOS technologies at low temperatures (this effect will be discussed in Chapter III), as well as the threshold voltage shift due to the back-gate biasing and the Fermi-Dirac shifting are captured in the 𝑉 𝑇0 model parameter adjustments, and the 𝐿 𝑠𝑎𝑡 parameter decreases at 4.2K because of the reduced phonon scattering. The normalized transconductance efficiency designmethodology is proved still valid for FDSOI transistors at 4.2K [17].

The published modeling works of the 28 nm FDSOI transistor based on the EKV MOSFET model [17], [18] do not account in an intrinsic manner for Fermi-Dirac statistics, do not demonstrate the C-V characteristics of the transistor, and deals only with singular channel operation.

❑ The MOS11/PSP models:

PSP is a physical surface-potential-based standard compact model for bulk MOSFET transistors that gives an accurate description of currents, charges and their higher order derivatives, an important feature for RF circuit designs [19]. The core model includes the computation of the surface potentials at both the source and drain ends, the terminal charges, and the intrinsic drain current accounting for the symmetric linearization method. The model contains two distinguished set of parameters, a global parameters set that consider the geometry dependencies, and a local parameters set that facilitates the parameter extraction procedure [20].

The PSP model is only certified down to the temperature of 218.15K. For liquid helium temperature operation, temperature-dependent model parameters need to be extracted from the LHT measurements, some of which need to be modified while others are zeroed. For the MOS11 models on the other hand, and besides the temperature-dependent extraction, additional electrical components need to be implemented in the model in order to consider LHT device behavior. For example, the kink effects observed in the 160 nm CMOS devices are taken into account by adding a non-linear resistor in series with the bulk terminal [21], [22]. Two examples of the output characteristics are shown in Figure 6, and a summary of the modified parameters is listed in Table 1. ❑ The BSIM model:

Berkeley Short-channel IGFET Model (BSIM) is an industry standard compact SPICE model.
BSIM is in fact a family of models that evolved along with the evolution of transistor's structure. BSIM-IMG (Independent Multi-Gate) is a surface potential-based for SOI/FDSOI MOSFETs, where the front and back surface potentials are derived simultaneously by solving Poisson's equation for different combinations of biases. BSIM-CMG (Common Multi-Gate) is a surface-potential based model for multi-gate MOSFETs (like FinFETs) that considers arbitrary channel cross-section shape. Once the computation of the surface potentials is performed, it is followed by the computation of the correspondent charges, capacitances, and terminal currents.

The BSIM models were reformulated to allow their operation at cryogenic temperatures. A new temperature dependent charge density model is proposed, taking into account the band tail states, threshold voltage, mobility, and current saturation. The results are shown in Figures 7 & 8 [23]. As for the drawbacks, the BSIM compact models are based on Maxwell-Boltzmann statistics, the impact of Fermi-Dirac statistics is emulated by considering an effective temperature instead, and a threshold voltage correction is added to the gate voltage term.

❑ The L-UTSOI model:

L-UTSOI is a standard surface-potential-based compact model dedicated to FDSOI MOSFET technologies with low-doped channel, developed at CEA-Leti and previously named Leti-UTSOI. The model accounts for the creation of a strong inversion layer at both interfaces of the silicon body, the so-called "dual-channel operation", and gives an accurate description of the currents, charges and their derivatives in all bias-configurations [24]. Moreover, the model reproduces accurately the normalized transconductance efficiency 𝑔 𝑚 𝐼 𝑑 ⁄ , a valuable feature for RF applications [24], [25]. Solving the dual-channel operation equations is challenging mathematically, as the integration of Poisson's equation with boundary conditions (the model consider volume inversion of only mobile charges following the Maxwell-Boltzmann statistics) leads to a set of equations involving either hyperbolic or trigonometric functions. The coupling of the two interfaces is expressed through a charge-dimensioned quantity, that can be either real (hyperbolic mode) if the interfaces are actually coupled, or imaginary (trigonometric mode) if the interfaces are decoupled [25].

At cryogenic temperatures, the L-UTSOI model cannot reproduce the I -V or C -V characteristics, and crashes sometimes if the temperature is below 173K. The temperature dependence of the parameters needs to be cancelled and a temperature-offset parameter is added, the rest of the parameters are fixed to work at 4K only.

Adapting the L-UTSOI model for deep cryogenic temperature on the other hand will result to a model produced C -V characteristics that seem to be good, as shown in Figure 9. However, the model produced transport characteristics especially for positive back-biases which are not accurate, as shown in Figures 10,11 & 12. The kink that the transfer characteristics show in Figure 11 is due to the Intersubband scattering, a phenomenon that is not manifested in higher temperatures but is found to be present at cryogenic temperatures for FDSOI transistors (more on that in Chapter 3). Empirical modifications have been introduced to the mobility law accounting for its degradation. Such an approach allows a better description of the kink effect and the produced transfer curves are more accurate. Moreover, to recapture for the inaccuracy observed in the output characteristics, Figure 10, a non-linear access resistance with symmetrical values at source and drain side is added. Such an approach has been introduced in other models, such as Leti-HSP [26], but its physical origin at cryogenic temperatures is not well understood. To summarize, even though the previous works seem to fit the experimental results of the transistor characterizations at cryogenic temperature, they all share the empirical-adaptation aspect of pre-existing compact models. As a result, Fermi-Dirac statistics, a fundamental aspect of electrons at those temperatures is either absent or post-emulated.

Moreover, neither of the FDSOI cryogenic operation published works demonstrate in an extensive manner the transistor electrostatics characteristics that is to say the electrostatic surface potentials for both interfaces, and the C -V curves. Likewise, the Fermi-Dirac statistics need to be considered in the transport part in the drain current derivation and not only in the electrostatics part. Furthermore, most of the published works maintain the RT mobility laws and adapt them to cryogenic behavior by modifying/adding fitting parameters or considering a constant mobility in the linear regime [17], [27], but none considers introducing a proper mobility law dedicated for MOSFETs operating at cryogenic temperatures [28].

Besides, numerical simulations are an important step that must precede the development of an analytical model and that will serve as a solid ground for the validation of its results/approximations. No numerical simulation has been held in such conditions, whether the simulation of the electrostatics at equilibrium, or electronic transport out of equilibrium. The use of such simulations is very important and necessary for the development of compact models. Firstly, because those simulations take very limited approximation and allow the demonstration of some quantities that are unachievable through experimentation, or they can be very handy when the experimental data required to calibrate the compact models are not available. In the frame of this work, Poisson-Schrodinger simulations of the FDSOI structure at cryogenic temperatures and accounting for Fermi-Dirac statistics were held to exploit and understand the device's electrostatics primarily, and then the simulations were held with the introduction of a quasi-Fermi level term to analyze the transport aspect of the device.

The scope of this thesis therefore, is to exploit the device's physics of FDSOI transistors down to deep cryogenic temperatures and to develop an appropriate core compact model. Such work would be the foundation into a compact model that is suitable for spice simulations.

For such work, we will be facing the following challenges:

• Understanding the underlying physics.

• Considering the dual channel operation of back-biased FDSOI transistor, which comes with its own challenges as we have the co-existence of two coupled interfaces. • The compact model must work for all regions of operation, whether in weak, moderate, or strong inversion, with a smooth transition between them. • The model needs to work for different configurations. That includes geometrical configurations, such as the oxides thicknesses, the silicon channel thickness, the channel length and width. The model needs to work for all biases configuration, in other words for a negative, null or a positive back bias. Finally, even though the model is dedicated for cryogenic operation, it must cover a good range of temperatures.

Chapter two:

The underlying physics at cryogenic temperatures Whereas in the previous chapter we detailed the technological aspects from which emerges the challenges that we are trying to solve with this thesis, the physical features, which should be our starting point, were not discussed. Thusly, the microscopic physical aspects are detailed in this chapter, the exposition of those aspects is put in order; i.e. the electrostatic-based parameters are discussed firstly, followed by the study of the transport-based parameters.

It should be noted that for some physical entities such as the carrier density 𝑛 𝑖𝑛𝑣 , both the physical and the numerical aspects are discussed. Such approach is necessary, as historically low temperature device modeling was considered extremely difficult, implicating, in first instance, our lack of understanding for the physical phenomena that appears at those conditions and the corresponding assumptions that have to be made, and in second instance, the numerical implementation that follows. Both of these two aspects are equally important. Some choices that are made for the rest of this thesis must be discussed in this chapter, such as the statistics used to describe the carriers distribution at cryogenic temperatures, the number of inversion channels involved in the electrostatic control or the transport phenomenon, the mobility law and the scattering mechanisms involved.

Moreover, in this chapter we consider the conduction band edge 𝐸 𝑐 as the energy potential reference, i.e. it will be represented by the point 0 in the equations/graphs hereafter.

Electrostatic-based physical parameters 1.Maxwell-Boltzmann approximation

In this section, the Maxwell-Boltzmann (MB) approximation validity down to cryogenic temperatures is investigated. Whereas some published literature emphasized on the non-validity of Maxwell-Boltzmann approximation down to cryogenic temperatures [1], [2], others argued its usage [3], [4]. To give an overall view of such aspect, the ongoing arguments in the scientific community on the validity of MB approximation at deep cryogenic temperatures will be specified first, followed by the reasons we believe the MB approximation does not hold at cryogenic temperatures and must be replaced by full Fermi-Dirac statistics.

To better address this debate, one must discern between the physical, the numerical, and the practical aspect of this subject. From a numerical viewpoint, the numerical overload in the exponential functions due to the very small 𝑇 is expected and inevitable. In addition, the intrinsic carrier concentration take very small values at deep cryogenic temperatures and becomes null at the 0K limit, following 𝑛 𝑖 = √𝑁 𝑐 . 𝑁 𝑣 . 𝑒𝑥𝑝(-𝐸 𝑔 2. 𝑘. 𝑇 ⁄ ) [5],(note that even 𝑁 𝑐 and 𝑁 𝑣 are temperature dependent, but such detail is irrelevant to our approach as will be shown). The scaling of intrinsic carrier concentrations down to deep cryogenic temperatures is not obvious; a number of solution has been proposed to surmount such problem [6], [7]. The 𝑛 𝑖 arithmetic underflow and the numerical limitations are tackled according to [3], [8] by using a variable arithmetic precision. Such approach is sound only if the non-explosion of the exponential factors is assured. Typically, 𝑛 𝑖 would be zero at the 𝑇 → 0𝐾 limit and takes extremely small values at deep cryogenic temperatures, for example, Figure 1 taken from [9] reveals that 𝑛 𝑖 values can go as low as 10 -278 𝑐𝑚 -3 at deep cryogenic temperatures, which is not meaningful. Jointly, at 0 K temperature, Fermi distribution function 𝑓(𝐸, 𝐸 𝐹 ) = 1 (1 + exp(𝐸 -𝐸 𝐹 𝑘𝑇 ⁄ )) ⁄ would be a step function, but a strict step function would lead to zero mobile carrier, which is inconsistent with the observed correct functioning of MOS transistors at deep cryogenic temperatures. However, the study performed on LDD structures [2] at low temperatures demonstrated the fallacy of such description, as in such structure, the series resistance rise up enormously because of the non-degenerate doping level at the LDD region, hence, the channel is deactivated. Nonetheless, when high enough electric fields are provided, either through the gate or the drain biasing, sufficiently impurity field ionization happens, resulting in a big reduction of the series resistance with biasing [2]. In other words, and in a more general case, even though the low temperature condition does not provide the required thermal energy for the electrons to be activated and open the inversion channel, high enough electric fields take charge of this task, and provide sufficient electrons to open the inversion channel [10], [11].

The physical viewpoint is addressed through two sub-levels; the first one is by inspecting the doping level, while the second one discusses the relative position of the Fermi level. Indeed, the assumption of holding the MB approximation is applicable in some specific cases where the doping level is very high, sometimes just below the degenerate limit and for temperatures down to 100mK when partial ionization is present 1 . Reason being that partial ionization would maintain the non-degeneracy of highly doped semiconductors, thus the MB approximation validity [3], [4].

In our case , the partial ionization argument does not maintain for an apparent reason, as being that though the channel doping is indispensable for bulk MOSFETs, it is not for ultrathin FDSOI architecture, where the doping is only considered as another degree of freedom to control the 𝑉 𝑡ℎ control of the device [12], [13] . Thus, the argument of maintaining the MB approximation due to the presence of partial ionization for highly doped channels does not hold, and the only remaining source of electrons to the channel is through the diffusion from source and drain regions. 1 The doping concentration is presumably fully ionized in room temperature, such assumption does not hold for low temperature simulation. Instead, the partial ionization is described using the next two formulas: Primarily, from a physical viewpoint as well, the best way to settle this debate is by analyzing the relative position of the quasi-Fermi level with respect to the conduction band edge. In such context, we have three main cases; each case is characterized by its own statistics:

𝑁 𝐷 + =
• If the Fermi level is situated deep in the band gap region i.e. 𝐸 𝐹 ≪ 0, then the Maxwell-Boltzmann approximation prevails [14] .

• If the Fermi level crosses the conduction band edge to some extent i.e. 𝐸 𝐹 ≥ 0, then the Maxwell-Boltzmann is not valid anymore and the electrons distribution is described by the Fermi-Dirac statistics [14], [15]. • If the Fermi level is well above the conduction band edge i.e. 𝐸 𝐹 ≫ 0, then the metallic behavior dominates and the electron statistics become fully degenerate [14], [15]. This is better illustrated in Figure 2, where the electron density is shown as a function of the difference between the Fermi level and the conduction band edge. Figure 2 indicates that at very low temperatures the electron density becomes very dependent to small changes of the Fermi level, meaning that MB approximation leads to huge overestimation of carrier densities as soon as the Fermi energy exceeds the band edge energy i.e. 𝐸 𝐹 -𝐸 𝑐 = 0. It is worth noting that the overestimation of carrier density will inevitably propagate to all density dependent quantities, the drain current principally. Another reason to keep the MB statistics is from a practical viewpoint. Considering the Poisson-Boltzmann equation as the standard starting point to calculate the corresponding electrostatics at deep cryogenic temperatures, [4] chooses to stick with the MB approximation to keep the model analytical along with its computational efficiency, with the assumption that Fermi-Dirac statistics would necessitate a numerical integration [4], [8]. In fact, as we will see in Chapter 5, the use of complete Fermi-Dirac statistics would not impede the development of an explicit model. Contrariwise, using Fermi-Dirac statistics inherent to cryogenic consideration have the advantage of the explicit formulation in both strong and weak inversions in a 2D system.

Furthermore, since the simulation at the 0K temperature are impossible with both Fermi-Dirac and Boltzmann statistics due to the infinitely large 1 𝑘𝑇 ⁄ , the Fermi-Dirac integral function should be replaced by a Heaviside function. This approach is coherent physically as the Heaviside function emulates perfectly the fully degenerate metallic statistics.

Quantization of the inversion layer

When a strong electric field is applied perpendicularly to the surface of the silicon channel, the electronic system forms two-dimensional energy bands called "subbands". Each subband corresponds to a quantized motion normal to the surface, with a continuum for motion in the plane parallel to the surface. Meaning that for an Si (100) surface, the perpendicular direction will be assisted by two valleys, which present the highest mass for electrons in motion, whereas the two parallel directions will be assisted by four valleys, which present the lowest masses for electrons in motion [16], [17] . The quantization effect becomes more accentuated for lowtemperature operation compared to room temperature operation [15], [18], [19]. To support such statement we have calculated the subband energetic positions on a 10 nm silicon layer of an FDSOI transistor by solving Schrödinger and Poisson equations self-consistently, without taking into account image charge and exchange-correlation contributions (full analysis will be reported in Chapter 3). The calculations were performed for ten subbands, one ground and nine excited subbands. However, as will be shown in 1.1.6, the electrons reside entirely in the three first subbands at cryogenic temperatures. We can see according to Figure 3, that the lower the temperature is, the higher the separation between the first and second subbands, and respectively between the second and third subbands becomes, signifying a more pronounced quantum effects as stated by [15], [18], [19] . Another subject that should be pointed out in the context of cryogenic temperatures modelling is the thickness of the inversion layer at those temperatures compared to room temperature.

From a classical viewpoint, the turning point, a parameter dependent on the temperature, is much closer to the interface for lower temperatures, following 𝑘 𝐵 𝑇/𝑞𝐹 𝑠 in the simplest approximation [15]. Such formula suggests that the lower the temperature is, the narrower the potential well is, yielding to an accentuated quantization.

A better way of explaining this might be through the introduction of the concept of the electrostatic screening length 𝜆 2𝐷 , a length that characterizes the distance over which any local perturbation of the electric potential can be attenuated by the free-electron charge of the twodimensional inversion layer. According to [20], 𝜆 2𝐷 takes a constant value of about 3𝐴 ° either at low temperature or at high carrier concentration i.e. when the degenerate statistics dominate.

With the quantum confinement, the inversion layer is not located strictly at the interface between the silicon body and the oxide, but a few Angstroms away because of the dark space effect. This has to be accounted for, but it does change the fact that the inversion layer is thin enough for the charge sheet approximation to be valid. Thusly, the charge sheet approximation [21], usually used in room temperature compact models, can also be considered as suitable for cryogenic temperature compact models, promoting a reasonable employment of this approximation in the sections dedicated to the numerical and analytical models (Chapter 4 & 5).

The two-dimensional density of electrons

The next physical parameter to be addressed is the density of states (DOS). As a direct consequence of the existence of a two-dimensional electron system, the attributed density of states is two-dimensional. The two-dimensional DOS is independent of energy 𝐸, and is given by Eq 2.1 and the electrons density is given by the integral demonstrated in Eq 2.2. In the absence of disorder, the density of states function is null below the energy level 𝐸 0 and becomes constant at 𝐸 0 and for higher energies. For the 〈100〉 oriented Silicon film, the valley pair pointing in the 〈100〉 direction, which has a degeneracy of 𝑔 = 2, have a longitudinal mass 𝑚 𝑙 = 0.91. 𝑚 0 and a transverse mass 𝑚 𝑡 = 0.19. 𝑚 0 . Note that, hereafter, when the DOS is treated as constant entity, it will be called 𝐴 2𝐷 .Concretely, if we inject the density of states in the integral given by Eq 2.2 to establish the expression of electrons density, we find Eq 2.3 .

Note that, assuming we are at the absolute zero limit 𝑇 = 0𝐾, the direct analytical development of Eq 2. In spite of this, in reality, the 2D subband is not a mere step function but it exhibits a band tail of states [22], [23]. The origin of such band tails is believed to be due to the potentialfluctuation-induced from crystalline disorder, residual impurities and strain, surface roughness etc. [22]- [24].

The best way to describe the band tail extension ∆𝐸 is through an exponential decrement of the DOS, following Eq 2.4, which is demonstrated through Figure 4. Note that this expression describes continuously the DOS function above and below the conduction band edge. We can see that if the energy level 𝐸 is positive, Eq 2.4 is reduced to the 2D DOS term i.e. 𝐴 2𝐷 , and when the energy level is negative, the exponential decrement of the 2D DOS is present. In contrast to the expression proposed by [23], where two different expressions were given to describe both of the DOS regions.

𝑁(𝐸, ∆𝐸) = 𝐴 2𝐷 1 + 𝑒𝑥𝑝(-𝐸 ∆𝐸 ⁄ )

Eq 2.4

Figure 4. Electronic density of states behavior in the presence of a band tail extension, the 0-point being the band edge.

Actually, the band tail extension ∆𝐸 = 𝑘. 𝑇 𝑠 is a finite characteristic temperature that we cannot go below (a behavior that will be demonstrated in the next section). In other words, this is equivalent to say that our system has an inherent physical limit.

The subthreshold slope saturation

The band tail of states becomes intriguing if one is interested on what happens in the subthreshold region. Such band tails are believed to be the origin of the saturation of the subthreshold slope 𝑆𝑆 = 𝑑𝑉 𝑔 𝑑𝑙𝑛(𝐼 𝑑 ) ⁄ demonstrated at cryogenic temperatures [22], [23]. There are two approaches to calculate the drain current involved in the subthreshold slope derivation. The first one supposes a proportionality between the drain current and the carrier density, such approach is the one used here. The second approach considers directly the electronic conductivity, and the transport in the 2D subband is described using the Kubo-Greenwood formalism [22].

The virtue of the first approach is that, assuming the electronic mobility is constant, the subthreshold slope is directly found as a function of the electron density, i.e. 𝑆𝑆 = 𝑑𝑉 𝑔 𝑑𝑙𝑛(𝐼 𝑑 ) ⁄ = 𝑑𝑉 𝑔 𝑑𝑙𝑛(𝑛 𝑖𝑛𝑣 ) ⁄ . Therefore, the subthreshold slope can be written as in Eq 2.5.

𝑆𝑆 = 𝑑V g 𝑑𝑙𝑛(𝑛 𝑖𝑛𝑣 ) = 𝑑V g 𝑑V s . 𝑑V s 𝑑𝑙𝑛(𝑛 𝑖𝑛𝑣 ) = 𝑑V g 𝑑V s . 𝑛 𝑖𝑛𝑣 𝑞 . 𝑑𝐸 𝐹 𝑑𝑛 𝑖𝑛𝑣 Eq 2.5
As the subthreshold slope is computed in the weak inversion regime where the electron density is neglected, the gate charge conservation of a single gate FDSOI transistor (which is also valid for bulk MOS transistors) is as follows:

C ox . (V g -V fb -V s ) = 𝑄 𝑑 (V s ) + C it . (V s -V 0 ) Eq 2.6
The total derivative of Eq 2.6 gives:

C ox . (dV g -dV s ) = 𝑑𝑄 𝑑 (V s ) + C it . dV s Eq 2.7
Where C it = 𝑞. N it is the interface trap capacitance, and V 0 is the channel Fermi potential.

Considering that by definition the depletion capacitance is C d = 𝑑𝑄 𝑑 (V s ) dV s ⁄ . It should be noted that, in the case of an FDSOI structure, the depletion capacitance is the equivalent capacitance of the silicon film capacitance C si and the back oxide capacitance C ox2 in series, i.e. C d = C si . C ox2 (C si + C ox2 ) ⁄ . So after rearranging the equation, we have:

dV g dV s = C ox + C d + C it C ox Eq 2.8
On the other hand, the expression of 𝑑𝑛 𝑖𝑛𝑣 𝑑𝐸 𝐹 ⁄ can be derived from Eq 2.2, yielding to:

𝑑𝑛 𝑖𝑛𝑣 𝑑𝐸 𝐹 = ∫ 𝑁(𝐸, ∆𝐸).(- 𝜕𝑓 𝜕𝐸 (𝐸, 𝐸 𝐹 )) +∞ -∞ 𝑑𝐸 Eq 2.9
Thus, by substitution in Eq 2.5 we have: for Boltzmann statistics i.e. for the limit 𝑇 ≫ 𝑇 𝑠 .

𝑆𝑆(𝑇) = 1 𝑞 . ∫
Using Eq 2.10, the dependence of the subthreshold slope with temperature is illustrated in Figure 5. We can see the linear dependence at high temperature (𝑇 > 𝑇 𝑠 ) where the band tail influence on the whole DOS is negligible, and the saturation at low temperature (𝑇 < 𝑇 𝑠 ) where the band tail is significant to the DOS with degenerate statistics, where 𝑇 𝑠 represents the transition between the plateau and the linear regions.

Note that, depending on technology, 𝑇 𝑠 lies in the range 30 → 50𝐾. Strictly speaking, for 28 nm bulk CMOS transistor technology it was found to be 46K [3], for 28 nm FDSOI transistor technology it was found to be 35K [23]. Likewise, an effective temperature function can be extracted by multiplying Eq 2.10 by the term q. C ox C ox + C d + C it ⁄ , yielding to Eq 2.11 that gives the real temperature for 𝑇 > 𝑇 𝑠 , and gives the saturation temperature if 𝑇 < 𝑇 𝑠 , as demonstrated in Figure 6. Such expression is advantageous as it reproduces the 𝑆𝑆(𝑇) saturation without the need of introducing a band tailed DOS but rather an effective temperature for the system. Such effective temperature approach is suitable for the development of the analytical model in Chapter 5.

Band gap widening:

The band gap 𝐸 𝑔 is another physical aspect that has to be taken into consideration in the case of deep cryogenic temperatures. Band gap widening at low temperatures is an effect that has been agreed upon in literature [25], [26], though varied expressions has been proposed to obtain its corresponding value, we choose to work with the expression given by Eq 2.12, proposed in [25], though it should be noted that all expressions gives the same value of 1.17 eV at 4K.

𝐸 𝑔 = 1,17𝑒𝑉 + 1,059. 10 -6 𝑒𝑉. (𝑇) -6,05. 10 -7 𝑒𝑉. (𝑇) 2

Eq 2.12

1.1.6 Single subband approximation:

At the zero temperature limit, within the usual gate biasing range, and provided the Fermi level is above the band edge, only the lowest energy level i.e. the ground subband is occupied. This postulate is extended, and sometimes for the sake of simplicity, for non-zero low temperatures. Such approximation is widely accepted and agreed on in literature [15], [20]. To support this interpretation we have calculated the electronic distribution 𝑛 𝑖𝑛𝑣 (𝑥) and the conduction band structure including the subbands positions by solving Schrödinger and Poisson equations selfconsistently. As stated before, the calculations were performed for ten subbands, one ground and nine excited subbands. However, as shown in Figure 7, at the 𝑉 𝑔1 = 1𝑉 limit, for a negative as well as for a null back-bias 100% of the carriers occupy only the ground subband (unprimed), the other subbands are left completely empty. Instead, for a positive back-bias 71.82% of the carriers resides in the ground subbands, 21.92% in the second subband (unprimed), and 6.25% in the third subband (primed). Therefore, for an SOI structure with thick silicon film, when a positive voltage is applied to both the front and back gates, potential wells are generated at both sides and the electrons are confined at both Si-SiO2 interfaces. These electrons contribute to the conductance independently due to their big enough spatial separation. This case is presented in Figure 8. b.

On the other hand, when the silicon film is very thin, both front and back interfaces reinforce the confinement of each other, yielding to a very large energy separation between the subbands. Consequently, the second subband becomes very high energetically and inaccessible, or with a very limited contribution, this case is presented in Figure 8. a. There is a third case presented in Figure 8. c, where the front and back interfaces couple to form the lowest two subbands. In this intermediate regime, the two subbands are energetically close enough, so that strong interactions between them are expected. Moreover, since the energy position of these two subbands is controlled via front and back gate biases, 𝑉 𝑔1 and 𝑉 𝑔2 , depending on their values there will be situations where the two first subbands are very close, generating very strong interactions, as will be discussed in 1.2.4.

Correspondingly, as an additional argument, the spatial distribution of electrons populating each of the subbands for the three structures discussed in Figure 8 are shown in Figure 9. We can see clearly that even though for the 3 nm case the volume inversion is very present [27], the 10 nm and the 16 nm structures demonstrate two separated inversion layers. Provided that the occupation of the two lowest subbands represents 93.74% of the carriers (Figure 9), and that only these two lowest subbands are found below the Fermi level in Figure 8, hereafter, and in the process of the compact model development, only two subbands will be considered. These two lowest subbands will be treated as the front and back layers, such choice will be further approved by the gate-to-channel curves behavior shown in Chapter 3.

Insofar only the electrostatic physical parameters have been discussed. In the second part of the chapter we will address the electronic transport coefficients in a 2D system, i.e. the conductivity, the diffusivity and the mobility. Such parameters are evaluated using the Kubo-Greenwood integral, by summing the parallel contribution of each subband and neglecting the intervalley scattering, at first order [2]. Presumably, the outstanding transport parameter that needs to be discussed is the mobility law, as it is of great importance for the modeling of both the transfer and output characteristics of the MOSFET device.

Transport physical parameters 1.2.1 The conductivity function

The Kubo-Greenwood formalism is a good start for the study of the MOSFET inversion layer transport properties, as it allows the computation of the transport parameters from the metallic regime to the semiconductor regime [20]. Generally, the macroscopic inversion layer sheet conductivity 𝜎 of an electronic system is obtained from the energy conductivity function 𝜎 𝐸 (𝐸) by the Kubo-Greenwood integral, following Eq 2.13 [20]. In this equation, 𝑓 is the Fermi function, and 𝜎 𝐸 (𝐸) is related to the energy mobility function 𝜇(𝐸) according to the Cohen's formula: 𝑑𝜎 𝐸 (𝐸) = 𝑞. 𝜇(𝐸). 𝑁(𝐸). 𝑑𝐸 [28], with 𝐸 being the carrier kinetic energy.

𝜎(𝑇, 𝐸

𝐹 ) = ∫ 𝜎 𝐸 (𝐸). (- 𝜕𝑓 𝜕𝐸 (𝐸, 𝐸 𝐹 )) +∞ -∞ 𝑑𝐸 Eq 2.13
Considering the band tails in the DOS function and a constant mobility 𝜇(𝐸) = 𝜇 0 in the aforementioned Cohen's formula, and integrating the energy conductivity function over all the energy range, the macroscopic conductivity function emerges, following Eq 2.14 .

𝜎(𝐸 𝐹 , ∆𝐸) = 𝑞. 𝐴 2𝐷 . 𝜇 0 . ∆𝐸. 𝑙𝑛(1 + 𝑒𝑥𝑝(𝐸 𝐹 ∆𝐸 ⁄ ))

Eq 2.14

Using Eq 2.14 the conductivity function is plotted in Figure 10, we can see clearly that the conductivity has the same behavior as the density of states in the first approach, as it exhibits an exponential tail below the band edge associated with that of the density of states. Note just (c) as well, the continued increase of the conductivity function above the band edge, an expected behavior due to the increase of the carrier kinetic energy [22]. In an equivalent manner to 1.1.4, the subthreshold slope 𝑆𝑆(𝑇) can be computed by direct consideration of the electronic conductivity. Note that the assumed proportionality between the drain current and the carrier density in the first approach is not valid at low temperatures when the statistics becomes degenerate [22], [START_REF] Ghibaudo | On the diffusion current in a MOSFET operated down to deep cryogenic temperatures[END_REF].

Our starting point for the second approach is the drain current equation of a MOSFET within the gradual channel approximation, given by Eq 2.15, 𝑦 being the space dimension that varies from source to drain.

𝐼 𝑑 = 𝑊. 𝑞. 𝜎(𝑇, 𝐸 𝐹 ).

𝑑𝐸 𝐹 𝑑𝑦

Eq 2.15

Considering 𝐸 𝐹𝑠 the Fermi level at the source point, and 𝐸 𝐹𝑑 the Fermi level at the drain point i.e. 𝐸 𝐹𝑑 = 𝐸 𝐹𝑠 -𝑞𝑉 𝑑𝑠 , 𝑉 𝑑𝑠 being the drain voltage, and after integration over 𝑦 between source and drain, we get:

𝐼 𝑑 = 𝑊 𝑞𝐿 ∫ [∫ 𝜎 𝐸 (𝐸). (- 𝜕𝑓 𝜕𝐸 (𝐸, 𝐸 𝐹 )) +∞ -∞ 𝑑𝐸] 𝐸 𝐹𝑑 𝐸 𝐹𝑠 𝑑𝐸 𝐹 Eq 2.16
Leading to, after some manipulation:

𝐼 𝑑 = 𝑊 𝑞𝐿 [∫ 𝜎 𝐸 (𝐸). 𝑓(𝐸, 𝐸 𝐹𝑠 ) +∞ -∞ 𝑑𝐸 -∫ 𝜎 𝐸 (𝐸). 𝑓(𝐸, 𝐸 𝐹𝑑 ) +∞ -∞ 𝑑𝐸]
Eq 2.17

Following the same steps as before 1.1.4, one obtains the second expression of the 𝑆𝑆(𝑇):

𝑆𝑆(𝑇) = ∫ 𝜎 𝐸 (𝐸). 𝑓(𝐸, 𝐸 𝐹𝑠 ) +∞ -∞ 𝑑𝐸 -∫ 𝜎 𝐸 (𝐸). 𝑓(𝐸, 𝐸 𝐹𝑑 ) +∞ -∞ 𝑑𝐸 ∫ 𝜎 𝐸 (𝐸). (- 𝜕𝑓 𝜕𝐸 (𝐸, 𝐸 𝐹𝑠 )) +∞ -∞ 𝑑𝐸 -∫ 𝜎 𝐸 (𝐸). (- 𝜕𝑓 𝜕𝐸 (𝐸, 𝐸 𝐹𝑑 )) +∞ -∞ 𝑑𝐸 . C ox + C d + C it 𝑞. C ox Eq 2.18
Provided we are in the saturation region of weak inversion i.e. when 𝑉 𝑑𝑠 ≫ 𝑘𝑇 𝑞 ⁄ or 𝑉 𝑑𝑠 ≫ 𝑘𝑇 𝑠 𝑞 ⁄ , the drain-side terms in Eq 2.21 and Eq 2.18 vanish, leading the next reduced form of the subthreshold slope:

𝑆𝑆(𝑇) = ∫ 𝜎 𝐸 (𝐸). 𝑓(𝐸, 𝐸 𝐹𝑠 ) +∞ -∞ 𝑑𝐸 ∫ 𝜎 𝐸 (𝐸). (- 𝜕𝑓 𝜕𝐸 (𝐸, 𝐸 𝐹𝑠 )) +∞ -∞ 𝑑𝐸 . C ox + C d + C it 𝑞. C ox Eq 2.19
The next figure compares the temperature dependence of the subthreshold swing as given by the two approaches in the weak inversion configuration, i.e. by Eq 2.10 and Eq 2.19. The two approaches were calculated with the same band tail extension (∆𝐸 ≈ 5𝑚𝑒𝑉 i.e. 𝑇 𝑠 ≈ 60𝐾) and considering a constant mobility function. Note that in this configuration, both the density of states and the conductivity function are proportional to 𝑒𝑥𝑝(𝐸 ∆𝐸 ⁄ ) yielding to a perfect match between the two approaches, as illustrated in Figure 11. One must not ignore that, all of this is true because we considered a constant mobility, as transport for localized states in the band tail might exist, making Eq 2.18 describing a more general case, as it considers a mobility function of energy. Generally, one shall distinguish between two cases, the case where the disorder is absent, for which the mobility edge 𝐸 𝑐 does coincide with the band edge 𝐸 0 , only delocalized electrons contribute to the conduction. The case where the disorder is present, for which the two edges 𝐸 𝑐 and 𝐸 0 do not coincide, and a mobility edge does exist. Besides the delocalized electrons above 𝐸 0 , we also have localized electrons below 𝐸 0 , which are believed to have a negligible contribution to the total conductivity, but rather, they might contribute to the total current by the so-called hoping conduction [15].

Nevertheless, it is worth noting that, from a practical viewpoint, the band tails are extracted by measuring the drain current variation with 𝑉 𝑔 . On the other hand, it can be extracted for example by integrating the gate-to-channel capacitance over the gate voltage in the weak inversion region, but such approach lacks the sufficient dynamic range for the accurate extraction of the band tail extension from the derivative of 𝑛 𝑖𝑛𝑣 (𝑉 𝑔 ) curves, making it less practical [22].

Although as a consequence of the existence of band tails in the DOS function, the conductivity function takes the same form, no physical condition decree that the two functions has the same length of band tails. In other words, the case where ∆𝐸(𝑛 𝑖𝑛𝑣 ) > ∆𝐸(𝜎 𝐸 (𝐸)) may exist, meaning that some states do not have the same opportunity to contribute to the transport like the rest of the states.

The diffusivity function:

Another mechanism that dominates in the subthreshold region is the electron diffusion. Such parameter is crucial for the understanding of the MOSFET operation in this regime. According to Einstein's relation, and within the assumption of a constant mobility, the electron diffusivity should decrease with temperature, following Eq 2.20. Such relation is valid with the consideration of Boltzmann statistics [START_REF] Ghibaudo | On the diffusion current in a MOSFET operated down to deep cryogenic temperatures[END_REF]. In the case of Fermi-Dirac statistics, the electron diffusivity is given by the so-called generalized Einstein relation, Eq 2.21. Finally, in the case of full degeneracy i.e. metallic statistics, and considering we have a single 2D subband, the electron density can be directly given by 𝑛 = 𝐴 2𝐷 . 𝐸 𝐹 , and the electron diffusivity will be given by Eq 2.22 [START_REF] Ghibaudo | On the diffusion current in a MOSFET operated down to deep cryogenic temperatures[END_REF].

𝑞. 𝐷 = 𝜇. 𝑘𝑇

Eq 2.20

𝑞. 𝐷 = 𝜇. 𝑛.

𝜕𝐸 𝐹 𝜕𝑛

Eq 2.21

𝑞. 𝐷 = 𝜇. 𝐸 𝑓

Eq 2.22

To obtain the variation of the diffusivity with inversion charge density in a single 2D subband, i.e. in the case of Fermi-Dirac statistics, Eq 2.23 is injected in Eq 2.21, yielding to the expression in Eq 2.24. This variation is illustrated in Figure 12, where at high temperatures the electron diffusivity demonstrates a linear dependence with temperature, whereas at low temperatures the electron diffusivity saturates to a constant value given by 𝑞. 𝐷 = 𝜇. 𝐸 𝐹 . Therefore, metallic statistics governs the electron diffusivity at low temperatures, and the classical Einstein relation is inadequate [START_REF] Ghibaudo | On the diffusion current in a MOSFET operated down to deep cryogenic temperatures[END_REF]. Accounting for an energy density of states with an exponential band tail, i.e. Eq 2.3 and the electrons density is given by Eq 2.3, where 𝑇 0 is a characteristic temperature, so as when 𝑇 ≪ 𝑇 𝑠 , the energy density of states and electrons density are always given by Eq 2.3 & Eq 2.3 , whatever the real temperature is. Thus, using the definition given by Eq 2.21, the corresponding diffusivity can be established as in Eq 2.25. Using this equation, Figure 13.a. illustrates the variation of the diffusivity with Fermi level for a constant mobility. We can see clearly that if the Fermi level is well above the conduction band edge i.e. 𝐸 𝐹 ≫ 0, the diffusivity varies linearly with the Fermi level in accordance with the metallic statistics predicted by 𝑞. 𝐷 = 𝜇. 𝐸 𝐹 . Whereas, for 𝐸 𝐹 ≪ 0, the diffusivity saturates at the value 𝐷 = 𝜇. 𝑘𝑇 𝑠 𝑞 ⁄ but with 𝑇 0 playing the role of temperature. The reliability of the assumption of a constant mobility might be debatable, as it seems too simplistic at first sight and need to be examined. To this end, using the relation in Eq 2.24 we

𝑛
(a) (b)
derived the experimental diffusivity from typical experimental effective mobility data obtained by split C-V technique on 28nm FDSOI MOSFETs, to produce the curves in Figure 14. As demonstrated in the figure, even though the effective mobility is not constant versus electron density, the overall variation of the associated Diffusivity exhibits a linear trend versus electron density. Such behavior is well fitted by the metallic statistics limit of 𝑞. 𝐷 = 𝜇. 𝐸 𝑓 computed with a constant mobility. Thus, the assumption of constant mobility for the diffusion computation is justified and can be kept. Hereby, in the section dedicated for the drain current integration in Chapter 5, the mobility is retained constant for the computation of the diffusion current component and considered as a function of inversion charges for the computation of the drift current component.

The mobility function:

The room temperature mobility law has been comprehensively documented for decades and many physical mechanisms were given to explain it [START_REF] Ghibaudo | Low temperature characterization of silicon CMOS devices[END_REF]. With respect to low temperature modeling, almost all of the modeling works that has been published retain the physical laws used for room temperature operation and adapt them for deep cryogenic operation [START_REF] Pahwa | Compact Modeling of Temperature Effects in FDSOI and FinFET Devices down to Cryogenic Temperatures[END_REF], [START_REF] Incandela | Characterization and Compact Modeling of Nanometer CMOS Transistors at Deep-Cryogenic Temperatures[END_REF]. It has been known for decades so far that the mobility laws used for room temperature modeling are not valid for cryogenic temperatures [2], [START_REF] Rais | On the High Electric Field Mobility Behavior in Si MOSFET's from Room to Liquid Helium Temperature[END_REF].

The mobility law 𝜇(𝑄 𝑖𝑛𝑣 ) is given directly by the different scattering mechanisms manifested at any given structure or temperature. This is widely known as the Mathienssen's rule, where the mobility law is the inverse of the sum of different scattering contributions. The main scattering mechanisms that govern the inversion layer mobility in a MOSFET are the Coulomb, surface roughness, and phonon scattering modes [2]. Although the Coulomb and the surface roughness mechanisms prevail at low temperatures, the phonon scattering is not present. The Coulomb scattering rate is inversely proportional to the impurity charge density and has a linear variation with energy [2]. According to Stern [START_REF] Stern | Calculated temperature dependence of mobility in silicon inversion layers[END_REF], combining Coulomb and surface roughness scatterings result in a bellshaped behavior of the mobility law, such thesis was confirmed by the ensuing studies [2], [START_REF] Ghibaudo | Transport in the inversion layer of a mos transistor: Use of kubogreenwood formalism[END_REF], [START_REF] Balestra | Physics and performance of nanoscale semiconductor devices at cryogenic temperatures[END_REF]. The generalized mobility law is given by [2], [START_REF] Rais | On the High Electric Field Mobility Behavior in Si MOSFET's from Room to Liquid Helium Temperature[END_REF], [START_REF] Emrani | Generalized Mobility Law for Drain Current Modeling in Si MOS Transistors from Liquid Helium to Room Temperatures[END_REF]:

𝜇 𝑛 (𝑄) = 𝜇 𝑚𝑎𝑥 ( 𝑄 𝑄 𝑐 ) 𝑛-2 1 + ( 𝑄 𝑄 𝑐 ) 𝑛-1 Eq 2.26
Where 𝜇 𝑚𝑎𝑥 is the maximum effective mobility, 𝑄 𝑐 is the critical inversion charge that characterize the mobility reduction, and the exponent 𝑛 ranges from 3 to 2 as the temperature ranges from very low to room temperature [2].

At cryogenic temperatures, the exponent n is equal to 3, and the mobility law is given by the expression in Eq 2.27. Indeed, the Coulomb scattering that rules the (𝑄 𝑄 𝑐 ⁄ ) behavior, and the surface roughness scattering that rules the (𝑄 𝑐 𝑄 ⁄ ) behavior.

𝜇 𝑛 (𝑄) = 𝜇 𝑚𝑎𝑥 ( 𝑄 𝑄 𝑐 ) + ( 𝑄 𝑐 𝑄 ) Eq 2.27
With the discussion of the mobility and the conductivity of DG-MOSFET or FDSOI-MOSFET, arises the subject of the number of transport layers involved at low temperature. One can say as a direct conclusion from the electrostatic part, that in the case of a back-biased FDSOI there is two layers involved in the conduction, and the overall conduction is the sum, in other words, the total current is the sum of the two front and back currents, such approach will be conducted in Chapter 5.

Yet another scattering mechanism:

Considering Coulomb scattering and surface roughness scattering as the only scattering mechanisms that prevails at deep cryogenic temperatures is largely true for a bulk MOSFET assuming that only one subband involved in transport. However, for a back biased FDSOI, one cannot ignore the interaction between the two subbands, since such interaction gives birth to an additional scattering mechanism, the so-called "inter-subband scattering" [10], [START_REF] Casse | Evidence of 2D intersubband scattering in thin film fully depleted silicon-on-insulator transistors operating at 4 . 2 K Evidence of 2D intersubband scattering in thin film fully depleted silicon-on-insulator transistors operating at 4 . 2 K[END_REF].

The starting point concerning this discussion should be the subband interaction mentioned formerly, as the decrement in the drain current observed in the experimental curves of transfer characteristics for a back-biased 28 nm FDSOI transistor at a temperature of 4.2K, Figure 15, cannot be explained by a mere sum of conduction from two subbands, but suggests that a subbands interplay took place.

According to [START_REF] Colinge | The new generation of soi mosfets[END_REF], for the intersubband scattering to occur two conditions must be satisfied. The first one concerns the temperature which should be low enough compared to ∆𝐸 𝑘 ⁄ . The second condition concerns the drain voltage which should be not much larger than ∆𝐸 𝑞 ⁄ . By definition, the intersubband scattering is absent if only one subband is occupied, which in our case is in the subthreshold region. As the front gate voltage is increased, the onset of the second subband occurs and it starts becoming populated, and the scattering occurs. Therefore, as conclusion, intersubband scattering should occur every time a new subband is populated [START_REF] Colinge | The new generation of soi mosfets[END_REF]. Poisson-Schrodinger simulations allowed us to inspect such behavior from an electrostatic point of view, as illustrated in Figure 16 below. One can clearly see that at the intersubband scattering effect appears when the energetic separation between the first two subbands is lower than 20𝑚𝑒𝑉, supporting the subbands interaction interpretation.

Chapter Three: Poisson-Schrödinger simulations

Every built model is based on some approximations that are meant to make the calculations feasible, easier, faster, or are used in some cases where incomplete information are known about a physical process. Depending on the level of modeling, such approximations make the predictions of the model differ from real measurements. Strictly speaking, each time an approximation is added, the model needs to be examined to justify such choice.

In our study, three levels of modeling are performed. The first one is of a numerical nature, which is based on the self-consistent solution of Poisson and Schrödinger equations. Such model is based on the so-called "effective mass approximation". These numerical simulations give us the chance to understand the electrostatic behavior of the FDSOI transistor exhibited in a C-V measurement. We shall call these numerical simulations simply as "PS simulations", in order to distinguish them from the following level model, of a numerical nature as well. In the second level of modeling, more approximations are made but the model is still on the numerical nature. This model will be referred to hereafter as "the numerical model". The third level of modeling is of an analytical nature, and a multitude of approximations is introduced on this level, the model will be referred to hereafter as "the analytical compact model".

Expressly, the PS simulations are used to validate the charge and current solutions predicted by the numerical model, which in turn is used to validate the same predicted solutions by the analytical compact model. Note that, in the frame of our work, PS simulations are a useful tool for understanding the physics that governs FDSOI transistors down to deep cryogenic temperatures. Nevertheless, for the purpose of not diverging from the aim of the study, none of the numerical aspects of the simulations are detailed, one should refer to [1]- [3] for such end. Furthermore, PS simulation results will also be comparison to collected data from C-V measurement, constructing a first checkpoint in the development process of our analytical compact model.

The simulation procedure:

Our Poisson-Schrödinger (PS) solver is based on a Python program solving self-consistently the Schrödinger and Poisson equations stated subsequently for an undoped FDSOI structure down to deep cryogenic temperatures and accounting for Fermi-Dirac statistics. Indeed, several published works have performed Poisson-Schrödinger simulations (PS) in bulk or FDSOI MOSFETs at room or low temperatures [1], [4], [5], but not down to very low temperatures. Simulations for temperatures less than 2K for example are lacking in literature.

Typical FDSOI structure used for PS simulation is shown in Figure 1, illustrating the stack composed of a front oxide of 1 𝑛𝑚 thickness, the silicon film (body) of 7 𝑛𝑚 thickness, and a buried oxide (BOX) of 25 𝑛𝑚 thickness. The front and back metal gates along with the highly doped regions of source and drain are illustrated as well. Note that, hereafter, the front and back silicon-oxide interfaces will be given the subscripts 1 and 2 respectively in the different mathematical equations. Considering an isotropic and parabolic conduction band, the effective mass approximation states that the behavior of an electron of mass 𝑚 0 in the crystal's potential is akin to an electron of an effective mass 𝑚 𝑒 for which the wave function is a Bloch function. Such electron of mass 𝑚 𝑒 moves in a periodic crystal lattice supposedly empty of ions. In other words, the crystal potential is substituted by the fact that the electron has an effective mass 𝑚 𝑒 different from the free electron mass 𝑚 0 [6].

For a 〈100〉 silicon-oriented structure, the calculations were performed for six valleys, a couple of unprimed valleys parallel to the confinement orientation, and two couples of primed valleys in the transverse plane perpendicular to the confinement orientation. For each series of valleys (primed and unprimed), ten subbands are arbitrarily considered in our case, one ground and nine excited subbands. The barrier height between the channel and both the front and back gate oxides is set to 3.1 𝑉, and the relative confinement effective mass of both the front and back gate oxides is set to 0.5.

The calculations are performed by solving two coupled non-linear partial differential equations (i.e. Poisson and Schrodinger equations) along with their boundary conditions. For most boundary value problems, the exact analytical solution is very difficult if not impossible to find except for some very special cases. Mainly, numerical methods based on iteration methods are used to find the solution.

O y

x

The Poisson equation is a differential equation of the second order that relates the electrostatic potential to the total carrier density. Such equation is solved using two boundary conditions with respect to the displacement field, as stated in Eq 3.1:

{ 𝛻(𝜀 𝑟 𝛻𝑉 𝑠 ) = -𝑞. 𝑛 𝑖𝑛𝑣 (𝑥) 𝜀 0 ⁄ 𝜀 𝑜𝑥1 (𝑑𝑉 𝑠 𝑑𝑥 ⁄ ) 0 -= 𝜀 𝑠𝑖 (𝑑𝑉 𝑠 𝑑𝑥 ⁄ ) 0 + 𝜀 𝑜𝑥2 (𝑑𝑉 𝑠 𝑑𝑥 ⁄ ) (𝑡 𝑜𝑥1 +𝑡 𝑠𝑖 ) -= 𝜀 𝑠𝑖 (𝑑𝑉 𝑠 𝑑𝑥 ⁄ ) (𝑡 𝑜𝑥1 +𝑡 𝑠𝑖 ) + Eq 3.1
The front and back surface potentials are deduced using the boundary definitions 𝑉 𝑠 (𝑡 𝑜𝑥1 ) = 𝑉 𝑠1 , and 𝑉 𝑠 (𝑡 𝑜𝑥1 + 𝑡 𝑠𝑖 ) = 𝑉 𝑠2 , and the boundary conditions of the electrostatic potential in the gates are given by 𝑉 𝑠 (0) = 𝑉 𝑔1 -∆𝜙 𝑔1 , and 𝑉 𝑠 (𝑡 𝑜𝑥1 + 𝑡 𝑠𝑖 + 𝑡 𝑜𝑥2 ) = 𝑉 𝑔2 -∆𝜙 𝑔2 . Note that in this work, we consider for the sake of simplicity midgap gates i.e. we choose a null front and back gate work functions ∆𝜙 𝑔1,2 with respect to midgap channel material.

The Schrödinger equation is solved in both the silicon film and oxides regions, with boundary conditions guaranteeing the electron wave function continuity at the Si/SiO2 interfaces, as stated in Eq 3.2. Such equation is solved for each subband within the effective mass approximation considering the electron envelope wave function [6].

{ - ℏ 2 2 ∂ 𝜕𝑥 ( 1 𝑚 𝑐𝑜𝑛𝑓,𝑣 (𝑥) ∂ 𝜕𝑥 𝜓 𝑠 𝑣,𝑙 (𝑥)) -𝑞𝑉(𝑥)𝜓 𝑣,𝑙 (𝑥) = 𝐸 𝑣,𝑙 𝜓 𝑣,𝑙 (𝑥) lim 𝑥→-∞ 𝜓 𝑣,𝑙 (𝑥) = 0 lim 𝑥→+∞ 𝜓 𝑣,𝑙 (𝑥) = 0 ∫ |𝜓 𝑣,𝑙 (𝑥)| 2 𝑑𝑥 +∞ -∞ = 1 Eq 3.2
It should be noted that the x coordinates start from the oxide edge and not from the oxide/film interface, and that an electron wave function penetrates in the oxides.

Finally, in order to finalize the package, one more equation is needed; namely, the total charge density, which is computed from the sum of contribution of all subbands, following:

𝑛 𝑖𝑛𝑣 (𝑥) = ∑ ∑ 𝑔 𝑣 𝐴 2𝑑𝑣 𝑘𝑇. 𝐹 0 ( 𝐸 𝐹 -𝐸 𝑣,𝑙 𝑘𝑇 ) . 9 𝑙=0 1 𝑣=0 |𝜓 𝑣,𝑙 (𝑥)| 2 Eq 3.3
The self-consistency of the calculations is achieved by using the following algorithm:

1. The Poisson equation is first solved with 𝑛 𝑖𝑛𝑣 (𝑥) = 0 as an initial condition.

2. The resulting potential distribution is fed into the Schrödinger equation to calculate the electronic wave functions and their corresponding energy levels, which represent the Eigen states and the Eigen values of Schrödinger's equations respectively. 3. Using this information, the electron concentration 𝑛 𝑖𝑛𝑣 (𝑥) is calculated using Eq 3.3. 4. The electron concentration is then re-introduced in the Poisson equation and a Newton-Raphson iteration process is used until convergence of the electron concentration is attained. 5. An error calculation step is performed between the electronic density obtained from the previous and the current iterations, if the convergence criterion for a variation of electronic density is satisfied, the convergence criterion is reached. Otherwise, the algorithm starting from step 2 is repeated until convergence is attained.

Appropriate keys for the interpretation of simulation results:

It was noticed that the overall electrostatic behavior exhibited by the different physical entities such as the surface potentials, the inversion charge and its derivative i.e. the gate to channel capacitance, can be attributed to the principal effects discussed in this section. These appealing effects are better stated beforehand the simulation results demonstration, as they provide the necessary tools to explain the different behaviors and tendencies in the subsequent sections. The present section is divided into two parts, a first part where we discuss back and front channel openings along with the evolution of back biasing. And a second part where we discuss the so-called "through the silicon coupling", that characterizes thin silicon films in the case of positive back bias configurations.

Back and front channels openings:

The openings of the back and front channels respectively are two main events that occur recurrently throughout this type of study and are the source of the major reported behaviors of the electrostatic parameters. Certainly, the existence of both channels is only true for positive back bias configurations, as for null and negative back bias configurations only the activation of the front channel is envisioned.

Comprehensively, the back channel opening corresponds to the onset of the first subband, and the front channel opening corresponds to the creation of the front potential well. Such attributions are only exclusively true for cryogenic temperature operation, as for higher temperatures several subbands are populated.

Such interpretation is confirmed by the exposition of electrostatic parameter curves along with the corresponding conduction band diagrams, to formulate a good approach for understanding these two main events. For example, one of the direct consequences of these two main events is the two-plateau behavior exhibited by the 𝐶 𝑔𝑐 curves in the case of positive back bias configurations. With two different elevations in two different polarizations, the variations of the gate-to-channel capacitance 𝐶 𝑔𝑐 (𝑉 𝑔1 ) = 𝑑𝑄 𝑖𝑛𝑣 𝑑 ⁄ 𝑉 𝑔1 with front gate voltage for various back gate bias, reveals the early onset of back inversion channel followed by the front channel opening for 𝑉 𝑔2 = +3𝑉, as shown in Figure 2. The correspondence between the polarizations of the back and front channel openings and the dual elevations in the 𝐶 𝑔𝑐 curves substantiate our interpretation. Further, this front channel opening should never be confused with the second subband onset, for those two events happen in two different polarizations.

Planely, one can see clearly from the comparison of Figure 2 and Figure 3 below, that the onset of the first subband of the unprimed valley 𝐸 0,0 (presented by the blue color) corresponds to the first elevation of the gate-to-channel capacitance curve. Note that, without loss of generality, midgap gates are considered in this study (as stated before), our potential's reference is indeed the position of Fermi level in the source, and the population of different subbands is activated once they reach the Fermi level and go beneath it.

The second elevation on the other hand is attributed to the creation of the front potential well, allowing the displacement of a portion of the electron gas to the indicated well. Such front potential well creation is followed by the onset of the second subband of the unprimed valley 𝐸 0,1 , presented by the orange color in the Figure .   The onset of the third subband on the other hand is presented by the undulation noticed beyond the polarization 𝑉 𝑔1 = 0.8 𝑉, such undulation is due to the difference of the valley masses since the third subband 𝐸 1,0 belongs to the primed valleys series characterized by a different mass.

Gate-to-channel capacitance curves as a function of front gate bias and for different back bias (t ox =1nm, t box =25nm, t si =10nm, T=4K)

.

First subband onset

Front potential well creation Otherwise, in the case of a negative or a null back bias, the electrostatic configuration permits the creation of a single potential well at the front interface, in which the totality of the inversion charge is confined when the onset of the ground subband is achieved, making the gate-tochannel curves exhibit only one plateau. Note that the roundedness feature exhibited by the 𝐶 𝑔𝑐 (𝑉 𝑔1 ) curves between the two plateaus is due to through-the-silicon coupling between the front gate and the back channel inversion layer discussed in the next section. Such interpretation will be further confirmed with the posterior study of the gate-to-channel curves performed for different silicon thicknesses.

In a related manner, we analyze the behavior of the remaining electrostatic parameters when variating the back bias polarization, as such behaviors are consistent with the formerly displayed Figures. Firstly, in Figure 4 we display the inversion charge curves as a function of front gate voltage at temperature 𝑇 = 4 𝐾, for a negative, null and a positive back bias. Note above all, the threshold voltage dependence on back biasing, which is basically due to the onset of the first subband. Expressly, the onset of the first subband happens first in the case of positive back bias because the positive polarization from the two sides (front and back) lowers the first subband, in comparison to the null or negative back biases.

Second subband onset Third subband onset

Moreover, Figure 5 demonstrates the front surface potential as a function of front gate voltage and for different back biases. Naturally, the positive back bias curve present three limbs compared to the two limbs presented by the null and negative back bias curves. Note as well that, once the onset of the first subband occurs for null and negative back bias curves they join up the positive back bias curve in the strong inversion region. In addition to the onset of the different subbands, the behavior of the electron gas inside the silicon film is addressed in this section. In Figure 7 is illustrated the band diagram across the stack and the electron density profile in the channel obtained from PS simulation at T=4K and for a given bias condition. Severely, only the ground subband (green line) is presented here. From the previous Figure, one can discern that when the applied back bias is negative, the structure presents a single quantum well and the electron gas resides mostly in the front side of the film. When the applied back bias is null, the physics is very similar except that a portion of the electron gas is not pulled away from the backside due to the absence of the negative polarization and can now reside in the backside of the film. When on the other hand a positive back bias is applied, the structure presents a two coupled quantum wells and since the applied polarization is positive from the two sides, the electron gas is pulled in the two directions. It should be noted that this displacement of the electron gas from the back side to the front side in the confinement direction is due to the "less strict" confinement behavior exhibited by the device in the case of the co-existence of two quantum wells i.e. in the case of positive back biasing.

Previously we showed in 2.1.6 that for extremely thin silicon films (𝑡 𝑠𝑖 = 3 -5 𝑛𝑚) only one inversion layer is manifested in the so-called volume inversion mode, and that for moderately thin and larger silicon films (𝑡 𝑠𝑖 ≥ 6 𝑛𝑚) the device displays two inversion layers that could be overlapped or not, depending on the film thickness. Yet, no particular indication was given for the overlap of the two inversion layers present in moderately thin silicon film geometries. Such overlap of the two inversion charges present for positive back biases is the source of the exhibited roundness by the gate-to-channel curves between the two plateaus. Since, thinner silicon films are at once characterized by stronger through-the-silicon front and back channel coupling and an accentuated roundness between the two plateaus. Evolution of gate-to-channel curves with silicon thickness listed thereafter will further confirm this interpretation. As shown in Figure 8 captured in the electrostatic symmetry configuration, when the silicon film is wide enough (𝑡 𝑠𝑖 = 16 𝑛𝑚), the two inversion charges are spatially separated and the overlap is nonexistent. In the case of a thin silicon film (𝑡 𝑠𝑖 = 7 𝑛𝑚) the two inversion layers are so close and overlap in the mid-region of the film. In first instance, throughout this chapter, we will study the effect of the existence of this overlap region on the electrostatic behavior of the device. Subsequently, for the numerical and the analytical model, the charge dwelling in the overlap region will be considered as a coupling charge, resulting from the capacitive coupling between the two interfaces.

Simulation results for different temperatures and silicon thicknesses:

To get the most of PS simulations, an extensive exhibition of the influence of different configurations, such as the silicon channel thickness in this context, as well as the temperature, are presented in this section. Generally, the electrostatic parameters demonstrate correspondent responses to each configuration's variation. Such responses can be consistently explained using the interpretations discussed in 1.2.

The evolution with temperature:

In a correlated manner to section 1.1, our PS simulations were even made possible at the 𝑇 → 0 𝐾 limit by replacing the 𝐹 0 Fermi-Dirac integral function by a Heaviside function. Such choice will allow us to emulate the fully degenerate metallic statistics. For temperatures other than the 𝑇 → 0 𝐾 limit, we retain the classical 𝐹 0 Fermi-Dirac integral function. Figure 9 suggests that the lower the temperature is, the steeper the subthreshold slope is, which in turn guarantees a faster transition of the transistor between the off and on states. Another direct by-product of the steeper subthreshold slope is the shallow moderate inversion region that lessens as the temperature lowers, until its total vanishing at the 𝑇 → 0 𝐾 limit.

Finally, it should be noted that this study does not consider the band tail states exhibited by the density of states function as discussed in 2.1.3 (previous chapter), allowing therefore the attainment of an ideal subthreshold slope in the 𝑇 → 0 𝐾 limit. Undeniably, such perfect subthreshold swing will never be observed in reality due to the subthreshold slope saturation discussed in the preceding chapter.

Supplementary information that PS simulations could provide is the charge centroid position evolution as the temperature lowers. Figure 10 shows such evolution from room temperature to deep cryogenic temperature T = 1 K. Since such information goes along with the dark space width evolution when the temperature decreases, such behavior confirms the expected emphasized confinement at lower temperatures. Note that, the dark space width will be thereafter (in Chapter 5) an inherent key to describe the quantum effects throughout an appropriate quantum shift function. The elbow characterizes the opening of the front channel and the consequent transition of the electron gas (populating the first subband) from the back channel to the front one. Likewise, this elbow characterizes this transition only at very low temperatures, due to the fact that at such temperatures only one subband is populated, compared to higher temperatures where several subbands are populated, as shown in Figure 12, making such transition to occur in a smoother way.

Note the shallower elbow for 𝑡 𝑠𝑖 = 7 𝑛𝑚 silicon film compared to the 𝑡 𝑠𝑖 = 10 𝑛𝑚, that is due in first place to the lower barrier that the electron gas has to overreach to attain the front side for the 𝑡 𝑠𝑖 = 7 𝑛𝑚 case. The height of the barrier exhibited by the conduction band in the middle of the silicon film depends fundamentally on two factors, the silicon film thickness, and the applied front and back polarization.

Note as well that, this transition occurs only for positive back bias configurations, because in those configurations the electron gas is less confined in the front potential well (relatively to the null and negative back bias configurations), and its displacement in the confinement direction is allowed once the front potential well is created. Thusly, the elbow is definitely not observed for null or negative back biases for at these configurations only the front potential well is created, and the electron gas is more confined in the corresponding potential well. Finally, one should point to the fact that this elbow will be the source of many numerical problems, at first instance, on the computation of the gate to channel capacitance in Chapter 4. The transcending of such numerical problems requires a special inspection phase, as will be discussed thereafter.

Figure 13 demonstrates the influence of temperature on the 𝐶 𝑔𝑐 (𝑉 𝑔1 ) curves; one can see the influence of temperature on the shape of the gate-to-channel curves. Manifestly, the lower the temperature, the more straight elevations we have, and the corners become more abrupt. Such effect can be due equivalently to the inversion charges overlapping in the middle of the film one more time, since this overlapping in the middle of the film is more present for higher temperatures as shown in Figure 14. In other words, lowering the temperature has a similar effect (with an inferior scale) to making the silicon film larger. One can see from analyzing the previous Figure, that in the case of a positive applied back bias, the front surface potential exhibits three limbs and two elbows, which manifest in the equivalent polarizations to the corresponding onset of first subband and creation of front channel elbows exhibited by the inversion charge in Figure 11. Note that for higher temperatures, these elbows vanish, and the surface potential drops due to the less present inversion charge in the close-tointerface region (as, for higher temperatures the device exhibit more volume inversion behavior).

If on the other hand the applied back bias is null, no creation of a second well is expected, thusly the curves exhibit only one elbow, which corresponds to the onset of the first subband.

Equivalently, Figure 16 demonstrates the back surface potential as a function of front gate voltage for three different temperatures. At odds with the front surface potential, the back surface potential exhibits two limbs. Note that the onset of the first subband designates the threshold i.e. the limit between the weak inversion region, presented by the linear trend, and the strong inversion region, presented by the saturated trend. Note as well, the smoother transition between the two trends compared to the abrupt one at cryogenic temperature. Such behavior is principally due to the number of subband involved in the transition. For instance, only one subband at the cryogenic temperature, and several ones for room temperature.

The evolution with silicon thickness:

Remaining in the positive back bias configuration, the influence of variating the silicon thickness on the behavior of the inversion charge is illustrated in Figure 17. Note the first subband onset dependence on silicon thickness. This is an expected effect, due to the energetic rise of the first subband when the film thickness is reduced. Moreover, as shown in Figure 17, the creation of the front potential well happens at the same polarization independently of the silicon thickness. Such effect is foreseeable as well, due to the fact that the bending of the conduction band depends only on the applied electric field and not on the film thickness. In a similar manner to Figure 19, the impact of the silicon channel thickness on the front surface potential at T = 4 K and for a positive back bias 𝑉 𝑔2 = +3𝑉 is reported in Figure 20.

One can see from analyzing Figure 20, that the onset of the first subband happens earlier for larger silicon films compared to the thinner ones, yielding to a shorter first limb of the 𝑉 𝑠1 (𝑉 𝑔1 ) curves. On the other hand, engendering a shorter second limb for thinner silicon films is due to the fact that the third limb stays rather unaffected to the variation of silicon film thickness. This can be explained by the drop of silicon film capacitance for larger silicon films, making the back surface potential less affected by the front gate bias. Note that, in such geometrical configuration, similar behavior is expected by the front surface potential if we study the influence of variating the back gate bias while keeping the front gate bias positive. Similarly, Figure 21 reports the effect of variating the silicon channel thickness on the back surface potential at T = 4 K and for a positive back bias 𝑉 𝑔2 = +3𝑉. Note the dependence of the threshold point (the limit between the linear and the saturated trends) on the silicon film thickness, i.e. the threshold happens earlier to larger films compared to the thinner ones. Note as well that the curve corresponding to 𝑡 𝑠𝑖 = 7𝑛𝑚 presents a minute increase before joining the other curves at the same polarization that characterizes the opening of the front potential well. Such particular trend is due to the strong through-the-silicon coupling exhibited at this silicon thickness, allowing the perception of the opening of the front potential well by the back surface potential.

Comparison to the C-V measurements:

In order to give credence to the performed PS simulations, their corresponding 𝐶 𝑔𝑐 (𝑉 𝑔1 ) results need to be validated by comparison to experimental data. For this reason that the pre-last section of this chapter should be the exposition of the gate-to-channel behavior predicted by the model to the C-V measurements. Note that in accordance with the measurements, the simulations were performed for a GO1 FDSOI MOSFET, by considering an equivalent front oxide thickness EOT1 of 1.2 nm. Note as well that the gate work functions, which were considered null earlier, are modified in accordance with measurements. The front gate function is found to be -0.1 𝑉, and the back gate function is found to be 0.5 𝑉, which is consistent with the well located below the buried oxide in the tested capacitance.

On top of that, the chosen AC level for the measurements becomes a significant parameters at the cryogenic temperatures range, as its influence becomes ostensible on the slope of the first elevation of the 𝐶 𝑔𝑐 (𝑉 𝑔1 ) curves. For the present measurements an AC level of 40 mV was chosen, with a frequency of 100 KHz.

Expanding the 1-D PS solver to the non-equilibrium conditions:

At this point, one should introduce the notion of the quasi-Femi level 𝜙 𝑖𝑚 that describes the population of the electrons throughout the channel when the system is displaced from its equilibrium condition, i.e. when an external voltage is applied to the drain terminal. Such notion is considerably practical, as it allows using the same equations that describes the electron density in the equilibrium condition cases, to the non-equilibrium ones.

Expanding the 1-D PS solver goes back to introducing the quasi-Fermi level in Eq 3.3, such that it emerges as in Eq 3.4, and the inversion charge is now computed in several points throughout the silicon channel.

𝑛(𝑥) = ∑ ∑ 𝑔 𝑗 𝐴 2𝑑𝑗 𝑘𝑇. 𝐹 0 ( 𝐸 𝐹 -𝐸 𝑙,𝑣 -𝜙 𝑖𝑚 𝑘𝑇 ) . Moreover, the procedure of obtaining the current information consists of integrating the computed inversion charge at each step of the quasi-Fermi level vector, as shown in Eq 3.5. Nonetheless, such information is still valuable for validating the drift current component obtained by the numerical model in Chapter 4.

𝐼 𝑑 𝑡𝑜𝑡𝑎𝑙 = 𝑊 𝐿 𝜇 𝑛 ∫ 𝑄 𝑖𝑛𝑣 (𝑉 𝑔1 , 𝑉 𝑔2 , 𝜙 𝑖𝑚 ).

𝑉 𝑑 0 𝑑𝜙 𝑖𝑚 Eq 3.5

Note that, the mobility law is indeed not included in this procedure. Thusly we consider a constant mobility and, for the sake of simplicity and without loss of generality, we consider (𝑊 𝐿 ⁄ ). 𝜇 𝑛 = 1. We dedicate this chapter for the development of a numerical model that is at odds with the prior PS simulation of simplistic nature. Such model halfway between the PS simulations and the analytical compact model is crucial for the development of the latter, as it establishes the general equations involved in the charge and current computations.

Our numerical model is dedicated to undoped long channel FDSOI devices operating at cryogenic temperatures and includes certain approximations, namely the charge sheet, the triangular well, and the gradual channel ones. Such approximations are meant to make the model simpler and computationally efficient. The criterion for each approximation once included in the model, is the model ability to predict the device characteristics with enough accuracy.

Expressly, this introductory section is dedicated to the discussion of each of the abovementioned approximations, their corresponding justifications and advantages/limitations. Following, a brief introduction of a smoothing method meant for connecting functions with respect to a particular abscissa point is presented.

The charge sheet approximation compresses the inversion layer into a conducting plane of zero thickness, where the conductivity is controlled by the electron density per unit area, denoting the loss of all information that describe the charge spatial distribution in the x direction [1]. Such approximation has the advantage of leading into a very simple algebraic formula concerning the charge or the current of long channel devices that applies across all regimes, i.e. from weak to moderate and to strong inversion regime, without any clamping or parameter changing, and without any manifesting discontinuities in the charges, currents, or their respective derivatives. The second advantage is its straightforward extension to two-dimension device analysis; as such analysis will be needed for the implementation of short channel effects that will be addressed in chapter 5.

In contrast to the numerical solution of the Schrödinger differential equation where its application to any potential profile 𝑉(𝑥) is feasible, for a simpler numerical model, the potential shape 𝑉(𝑥) which originally has a semi-parabolic shape, has to be approximated. Our choices are limited to triangular or square potential well approximations. One can deduce based on the potential profiles demonstrated in the previous chapter that the front and back interfaces are better approximated using the asymmetric triangular well approach. Square well approximations on the other hand are more suitable for ultra-thin film or heterostructures of different semiconductors [2].

The triangular potential well approximation is widely used in the MOSFET modelling community. Such approximation considers an asymmetric triangular potential well, with an infinite vertical barrier at 𝑥 = 𝑡 𝑜𝑥1 , and a barrier that varies linearly for 𝑥 > 𝑡 𝑜𝑥1 in the front interface case, equivalently, an infinite barrier at 𝑥 = 𝑡 𝑜𝑥1 + 𝑡 𝑠𝑖 and a barrier that varies linearly for 𝑥 < 𝑡 𝑜𝑥1 + 𝑡 𝑠𝑖 in the back interface case. According to [3], owing to the absence of a depletion depth, such approximation is even more justified in the case of undoped films.

Such approach approximates the gradually changing slope of the electrostatic potential 𝑉(𝑥) into a constant slope using the notion of effective electric field 𝐸 𝑒𝑓𝑓1,2 . Whereas in the case of Poisson-Schrödinger simulation, the electron energies and states are established by solving Schrödinger equation based on the corresponding boundary conditions, in the case of triangular well approximation approach, these eigenstates are given by the analytical solutions to the Airy functions, as expressed in Eq 4.1 [3], [4]:

𝐸 𝑖1,2 = ( ℏ 2 2𝑚 𝑒 ) 1 3 ⁄ ( 3𝜋𝑞𝐸 𝑒𝑓𝑓1,2 2 (𝑖 + 3 4 )) 2 3 ⁄ Eq 4.1
Where the index 𝑖 corresponds to the different populated subbands. Indeed, in our case only the ground subbands are considered, furthermore, such equation can be expressed by the mean of the gate charge densities, as in Eq 4.2. Indeed, this approximation does not provide the carrier distribution profile nor the inversion layer centroid. However, it is rather advantageous in our case because it is computationally simple to implement, as it involves only the power function of the gate charge function [5].

𝐸
The gradual channel approximation assumes that the potential along the channel varies gradually and in a suitable way for the 1-D electrostatic solution to be valid. Accordingly, such approximation enables the partition of the 2-D problem encountered in device modeling into two separate 1-D problems, the electrostatic 1-D problem in the 𝑥 direction, and the conduction down the channel 1-D problem in the 𝑦 direction [1].

Supplementary, and for proper operation in circuit simulators, the expressions used in our model should be at least 𝐶 1 continuous, or ideally 𝐶 ∞ continuous. However, in practice, the current and charge expressions for different bias, temperature or geometrical configurations can be different due to the differing physics in each configuration and thus the chosen approximations to describe that physics. On top of this, different operational inversion regions, i.e. weak-moderate and strong inversion regions, may have different expressions as well.

On that account, employing piecewise smooth functions is a deliberate way in order to obtain globally continuous functions, such approaches are found consistently in the field of semiconductor device modelling. For instance, in the published work [6], a generalized formalism generating such piecewise smooth functions is presented. Such smoothing method which is meant for connecting globally continuous piecewise functions 𝑓 𝑖 (𝑥), in the vicinity of abscise coordinates 𝑥 𝑖 , in order to get an infinitely differentiable functions, i.e. the resulting piecewise function can be evaluated at ]-∞, +∞[ and are 𝐶 ∞ , and the final function 𝐹(𝑥) can be obtained by simply summing up all piecewise functions [6]. Ultimately, a state-of-the-art implementation of such method will be demonstrated in section 1.1.4.

The Charge numerical model: 1.The classical starting set of transcendent coupled equations:

One of the inherent features of SOI devices is that the electrical properties of the back interface are influenced by the ones at the front one [7]. Such coupling between the two interfaces will result in a feedback mechanism, where the surface potential is controlled not only by the adjacent gate but also by the opposite one [8], [9]. That is to say, modifying the back-gate bias directly changes the back surface potential and indirectly changes the front interface properties; reciprocally, since the front interface potential was altered, this in turn will cause an additional change of the back-channel potential [9]. In our mathematical formalism, this can be interpreted by the fact that the surface potential at one interface is a function of the applied bias on the opposite interface and vice versa i.e. 𝑉 𝑠1 = 𝑓(𝑉 𝑔1 , 𝑉 𝑔2 ) and 𝑉 𝑠2 = 𝑓(𝑉 𝑔1 , 𝑉 𝑔2 ). This is expressly due to the electrostatics that generates two-coupled quantum wells DCQW, which consist of two 2D electrons gases that are spatially close enough and energetically separated by a relatively narrow barrier. Thusly, the studied system will have an additional degree of freedom [10], [11]. Such additional degree of freedom produced by the DCQW will be manifested in our system of fundamental equations, as an extra degree of complexity, i.e. an extra capacitive coupling term must be taken into consideration, in addition to the gates and inversion charges terms in order to fully model the system.

From an electrostatic point of view, the coupling term can be derived in the subthreshold regime considering the capacitive scheme of the stack, as illustrated in Figure 1. This simple scheme considers the three capacitances of the system in series, i.e. the front oxide, silicon film, and back oxide capacitances. Such capacitive scheme is initially derived in the weak inversion mode, and will be rather then extended to the moderate and strong inversion modes [9].

Figure 1. The capacitor scheme of a long channel FDSOI transistor

Therefore, in our approach we model this coupling charge presented in the system using the silicon film capacitance. Such simple method is quite advantageous for a compact modelling approach, as it describes in a reasonable manner the coupling of the two interfaces, and allows to model the coupling with linear terms, as will be demonstrated in the starting set of equations [9]. Indeed, our choice of the coupling term is convenient, as it allows a good description of the coupling charges for thin silicon films where 𝐶 𝑠𝑖 becomes large, and equivalently for larger silicon films where 𝐶 𝑠𝑖 becomes small reflecting the negligible amount of coupling charge present in these configurations. Finally, it should be noted that such coupling of surface potentials is undoubtedly reflected as a coupling of threshold voltages as well [7].

Accordingly, if we consider the typical scheme of an FDSOI structure presented in Chapter 3 and apply the charge conservation principle to the first and second interfaces respectively, we will acquire that the charge in the front gate equals the inversion charge of the first interface, plus the coupling charge. Equivalently, the charge in the second gate will be equal to the inversion charge of the second interface minus the coupling charge, as featured in Eq 4.3:

{ 𝑄 𝑔1 = -𝑄 𝑖𝑛𝑣1 + 𝑄 𝑐𝑝𝑙 𝑄 𝑔2 = -𝑄 𝑖𝑛𝑣2 -𝑄 𝑐𝑝𝑙 Eq 4.3
Regarding the closed-form mathematical expressions of the gate charge terms, and the coupling charge term, they are expressed using:

{ 𝑄 𝑔1,2 = 𝐶 𝑜𝑥1,2 . (𝑉 𝑔1,2 -𝑉 𝑓𝑏1,2 -𝑉 𝑠1,2 ) 𝑄 𝑐𝑝𝑙 = 𝐶 𝑠𝑖 . (𝑉 𝑠1 -𝑉 𝑠2 ) Eq 4.4
Concerning the inversion charge expression we preliminarily choose its formulation within the classical assumption, as described in Eq 4.5. Such expression suggests that we consider having two delta functions representing the front and back inversion charges, which start populating once the edge of the conduction band taps the Fermi level.

-𝑄 𝑖𝑛𝑣1,2 = 𝑞𝑘𝑇𝐴 2𝑑 𝑙𝑛 (1 + 𝑒𝑥𝑝 ( 𝑉 𝑆1,2 -𝑉 0 -𝜙 𝑖𝑚 𝑘𝑇 ))

Eq 4.5

The flat band voltage 𝑉 𝑓𝑏1,2 takes into account the potential drop inside the corresponding oxide due to the difference of work functions between the silicon channel and the gate metal, along with the consideration of the different interface trap charges. One should emphasize that, without loss of generality, we consider the midgap, represented here by the term 𝑉 0 , as our chosen potential reference.

Essentially, by solving the system of coupled equations we mean finding the corresponding values of the Fermi level, or in our context of the front and back surface potentials (𝑉 𝑠1 , 𝑉 𝑠2 ), for each variation of the applied bias to the front and back gates respectively, i.e. the couple (𝑉 𝑔1 , 𝑉 𝑔2 ).

Note that in our system of coupled equations, we have three terms that are linear with respect to 𝑉 𝑆1 and 𝑉 𝑆2 , the corresponding 𝑄 𝑔1 , 𝑄 𝑔2 , and 𝑄 𝑐𝑝𝑙 , and two non-linear terms with respect to 𝑉 𝑆1 and 𝑉 𝑆2 , the corresponding 𝑄 𝑖𝑛𝑣1 and 𝑄 𝑖𝑛𝑣2 . Considering the transcendental nature or our equations, they require the application of some numerical method in order to be solved, such approach is the one we follow in this chapter. Subsequently, in Chapter 5 we seek into succeeding to find a procedure that permits their transformation into algebraic equations.

Our system of coupled equations is solved using a python script that calls the fsolve function from the scipy.optimize package. Such function employs Powell's method for roots findings, a procedure that is performed for each bias configuration until the convergence is achieved [12]. Manifestly, in order to converge to the roots, a starting set of initial values for both the front and back surface potentials is needed, a suitable initial guess in our case is the silicon film midgap implemented in the form of the (𝑉 0 , 𝑉 0 ) list.

Thus far, only the classical form of our system of coupled equations is considered, which can be plainly solved and provide the classical solutions. Nonetheless to obtain the exact solutions, a further step must be completed, the incorporation of the quantum shift function.

The incorporation of the quantum shift function:

Such quantum shift function incorporation is achieved by the simple substitution of Eq 4.5 by Eq 4.6 in our system of coupled equations i.e. Eq 4.3.

- where the prefactor 𝛽 𝑄𝑀 encompasses all the constant coefficients that appear in Eq 4.2 . Such approach suggests that we consider having two delta functions representing the front and back inversion charges. However, in contrast to the classical approach, the condition for different subbands to start being populated is for the Fermi level to reach the energy levels of these subbands.

Primarily, owing to our lack of information of the quantum shift function behavior in the negative gate charge region an even function of the gate charge is commonly considered in Comprehensibly the 2 3 ⁄ power law is not sufficient to describe the manifested physics in the negative gate charge region, nor around the zero-gate-charge configuration within the moderate inversion region. One of the employed solutions proposed in literature [13] is the shown sensibility of the model outputs on the attributed precision to the rational power of the Airy function in its float format. Such values depend as well on the designated interface, as different values can be given to the front and back interfaces. Such modification of Airy's function is explained in literature [13] by the degree of penetration of the electronic wave function in the oxide barrier.

It was found indeed that such approach is functional, though not practical. As even though we get rid of the singularity presented around the zero electric field configuration by modifying the power values, those same values depend on the oxide and silicon thicknesses. Such dependence is undesirable for us, firstly because the model is not robust geometrical-configurations-wise and secondly because in our study the symmetry front-back-interface wise must be preserved.

As an initial remedy to the numerical pathology that one envisages is the inclusion of gate charge offset in the 𝛥𝑉(𝑄 𝑔1,2 ) expression, following Eq 4.9. Such infinitesimal offset charge becomes functional once we pass by the null-gate-charge point.

𝛥𝑉(𝑄 𝑔1,2 ) = 𝛽 𝑄𝑀 . |𝑄 𝑔1,2 + 𝑄 𝑜𝑓𝑓𝑠𝑒𝑡 | 2/3
Eq 4.9

Figure 3 comprehends the different electrostatic quantities as functions of front gate voltage computed using the Equations: Eq 4.9, Eq 4.6, and Eq 4.3. Note the vanished numerical pathology ascribed to the existence of the charge offset term 𝑄 𝑜𝑓𝑓𝑠𝑒𝑡 .

Figure 3. Front and back surface potentials, total inversion charge density, and gate-to-channel capacitance curves as

functions of front gate bias computed using the Equations: Eq 4.9, Eq 4.6, and Eq 4.3.

Moreover, as an alternative remedy we examined the effect of the inclusion of the 𝛥𝑉(𝑄 𝑔1,2 ) expressions in the coupling term 𝑄 𝑐𝑝𝑙 , as in Eq 4.10.

𝑄 𝑐𝑝𝑙 = 𝐶 𝑠𝑖 . [(𝑉 𝑠1 -𝛥𝑉(𝑄 𝑔1 )) -(𝑉 𝑠2 -𝛥𝑉(𝑄 𝑔2 ))]

Eq 4.10

In Figure 4 we present the different electrostatic quantities as functions of the front gate voltage using the Equations: Eq 4.10, Eq 4.7, Eq 4.6, and Eq 4.3. Noticeably, regardless of the aptitude of such approach into eliminating the numerical pathology and allowing more numerical stability around the zero-gate-charge point, it provides a better description of the roundness of the gate-to-channel curves manifested between the back and front channel openings. However, it was found that, either attempts intended to overcome the pathology does not fully solve the problem, as the model remains susceptible for numerical pathologies when we change the geometrical configuration, for instance, the silicon film thickness or the front oxide thickness, implying the non-validity of such approaches. Hence, in the next section we follow a special approach in order to find an inherent remediation for such pathology.

The extended quantum shift function:

In order to transcend such numerical pathology, we rely on a reverse path approach stemmed from the PS simulation data. Methodically, starting form PS simulation results, and from Eq 4.6 we can practically predict the behavior of the quantum shift function 𝛥𝑉(𝑄 𝑔1,2 ). Such function can be expressed in terms of the front/back surface potentials and front/back inversion charge densities respectively as in Eq 4.11.

𝛥𝑉(𝑄 𝑔1,2 ) = 𝑉 𝑆1,2 -𝑉 0 -𝑘𝑇. 𝑙𝑛 (𝑒𝑥𝑝 ( -𝑄 𝑖𝑛𝑣1,2 𝑞𝑘𝑇𝐴 2𝐷 ) -1)

Eq 4.11

Note that, whereas the front/back surface potentials information along with the front/back gate charges information are provided by PS simulation, the latter gives the total inversion charge density without discerning the front and back components. Accordingly, in order to get that information, the inversion charge needs to be apportioned between the front and back interfaces with respect to the gate charges densities. Generally, we should effectively have the following three cases:

• If 𝑄 𝑔1 > 0 and 𝑄 𝑔2 < 0, then: 𝑄 𝑖𝑛𝑣1 = 𝑄 𝑖𝑛𝑣 and 𝑄 𝑖𝑛𝑣2 = 0.

• If 𝑄 𝑔1 < 0 and 𝑄 𝑔2 > 0, then: 𝑄 𝑖𝑛𝑣1 = 0 and 𝑄 𝑖𝑛𝑣2 = 𝑄 𝑖𝑛𝑣 .

• If 𝑄 𝑔1 > 0 and 𝑄 𝑔2 > 0, then: 𝑄 𝑖𝑛𝑣1 = -𝑄 𝑔1 and 𝑄 𝑖𝑛𝑣2 = -𝑄 𝑔2 .

Therefore, the front and back inversion charge densities can be expressed in terms of front and back gate charges following Eq 4.12: Noticeably, the plotted function manifests two distinctive regions, a linear region associated with the negative gate charges, and a 2/3 power region associated with the positive gate charges. Whereas the 2/3 power behavior can be directly attributed to the presence of an asymmetric triangular potential well in the strong inversion region, allowing the use of Airy's function solutions [3], [4], the linear behavior is seemingly inherent to the weak inversion charge where the gate charge densities are minor.

Therefore, in order to establish a global quantum shift function that describes the whole region of gate charge densities, we must primarily define two functions to describe the linear and the 2/3 power regions separately, then succeed into building the whole block through the implementation of the formalism described in [6].

First, the linear region can be expressed plainly by introducing a prefactor to the gate charge densities as in Eq 4. 13.The slope of the linear region is found to be a function of the silicon film thickness, capacitance, and an extra smoothing parameter 𝑑𝑘𝑠, which represents the dark space. In this context, the parameter 𝛾 is equivalent to the notion of effective silicon thickness 𝐶 𝑠𝑖 𝑒𝑓𝑓 = 𝜀 𝑠𝑖 𝑡 𝑠𝑖 𝑒𝑓𝑓 ⁄ , where 𝑡 𝑠𝑖 𝑒𝑓𝑓 is the on-hand silicon thickness after excluding the dark space i.e. 𝑡 𝑠𝑖 𝑒𝑓𝑓 = 𝑡 𝑠𝑖 -𝑑𝑘𝑠 , leading after some rearrangements to the expression in Eq 4. 14.

𝛥𝑉 𝑛𝑒𝑔 (𝑄 𝑔1,2 ) = (1 𝛾 ⁄ ). 𝑄 𝑔1,2
Eq 4.13

𝛾 = 𝐶 𝑠𝑖 (1 -(𝑑𝑘𝑠 𝑡 𝑠𝑖 ⁄ )) ⁄ Eq 4.14
Ensuing, the 2/3 power dependence of the quantum shift function is given directly using the Airy's function solution, following Eq 4.15. In addition, and in order to better fit the PS simulations, 𝛽 𝑄𝑀 can be moderately reduced through the multiplication by a factor ≈ 0.8. 𝛥𝑉 𝑝𝑜𝑠 (𝑄 𝑔1,2 ) = 𝛽 𝑄𝑀 . 𝑄 𝑔1,2 2/3 Eq 4. 15 Heretofore, we have two distinctive mathematical expressions to describe the quantum shift function in both regions. Such piecewise bi-regional function has to be connected smoothly in order to obtain a fully differentiable global function.

In order to build our global function, the first step is the smooth clamping of each primary function. Hence, we use the assistance of two auxiliary functions that incorporate the smoothing parameter 𝛿𝑄 𝑔 , following Eq 4. 16. Factually, 𝑄 𝑛𝑒𝑔 (𝑄 𝑔1,2 , 𝛿𝑄 𝑔 ) smoothly clamps 𝑄 𝑔1,2 to the upper limit 𝑄 𝑔1,2 = 0 when 𝑄 𝑔1,2 increases from -∞ to +∞. On the other hand, 𝑄 𝑝𝑜𝑠 (𝑄 𝑔1,2 , 𝛿𝑄 𝑔 ) smoothly clamps 𝑄 𝑔1,2 to the lower limit 𝑄 𝑔1,2 = 0 when 𝑄 𝑔1,2 increases from -∞ to +∞.

{ 𝑄 𝑛𝑒𝑔 (𝑄 𝑔1,2 ) = 1 2 . (𝑄 𝑔1,2 -√𝑄 𝑔1,2 2 + 𝛿𝑄 𝑔 2 ) 𝑄 𝑝𝑜𝑠 (𝑄 𝑔1,2 ) = 1 2 . (𝑄 𝑔1,2 + √𝑄 𝑔1,2 2 + 𝛿𝑄 𝑔 2 )
Eq 4.16

Decidedly, the two primary functions 𝛥𝑉 𝑛𝑒𝑔 (𝑄 𝑔1,2 ) and 𝛥𝑉 𝑝𝑜𝑠 (𝑄 𝑔1,2 ) have to be connected with respect to the abscissa axis in the vicinity of null gate charge point i.e. 𝑄 𝑔1,2 = 0. the approach proposed by [6] allows such finality.

Expressly, in our case, the global function is 𝛥𝑉(𝑄 𝑔1,2 ), the two primary functions are 𝛥𝑉 𝑛𝑒𝑔 (𝑄 𝑔1,2 ) and 𝛥𝑉 𝑝𝑜𝑠 (𝑄 𝑔1,2 ) respectively, the connection of the two functions is made in the vicinity of the null gate charge point i.e. 𝑄 𝑔1,2 = 0, the same smoothing parameter 𝛿𝑄 𝑔 is used for both lower and upper limit functions. The final function is defined as the direct sum of the two piecewise functions, following Eq 4.17:

𝛥𝑉(𝑄 𝑔1,2 ) = 𝛥𝑉 𝑛𝑒𝑔 (𝑄 𝑔1,2 ) + 𝛥𝑉 𝑝𝑜𝑠 (𝑄 𝑔1,2 )

Eq 4.17

Furthermore, the constant terms, (𝛿𝑄 𝑔 2 ⁄ ) for the linear region, and -(𝛿𝑄 𝑔 2 ⁄ ) Finally, using Eq 4.18, we can reproduce the quantum shift function behavior in the case of a positive back bias configuration predicted by PS results using the surface potentials and the gate charge densities provided by the numerical model. A comparison between the produced plots is illustrated in Figure 6 for two different silicon thicknesses.

Plainly, the function described in Eq 4.18 does well in reproducing the trend exhibited by the function 𝛥𝑉(𝑄 𝑔1 ) derived from PS simulations, especially regarding the smooth transition between the linear and the 2/3 power behaviors. Before initiating the section dedicated to the comparison between the numerical model results and the ones provided by PS simulations using two approaches that are distinctive to the aforementioned one, it should be noted that since the PS simulations solve Schrodinger equation in the entire meshed system, i.e. the silicon film alongside the wave function penetration zones, such computation return the total inversion charge profile in the silicon film, without any prior distinction between the front and back inversion charge densities.

In order to compute these two densities starting from PS data we have two approaches:

• The first approach splits the silicon film into two equilateral halves, where the inversion charge that dwells between 𝑥 = 𝑡 𝑜𝑥1 and 𝑥 = 𝑡 𝑜𝑥1 + 𝑡 𝑠𝑖 2 ⁄ is purported to be the front inversion charge, and the inversion charge that dwells between 𝑥 = 𝑡 𝑜𝑥1 + 𝑡 𝑠𝑖 2 ⁄ and 𝑥 = 𝑡 𝑜𝑥1 + 𝑡 𝑠𝑖 is purported to be the back inversion charge.

• The second approach exploit the front and back surface potentials provided by the PS simulations and evaluate the front and back inversion charge densities using the gate charge and coupling charge terms as indicated in the equations Eq 4.3 and Eq 4.4.

Figure 7 incorporate the front and back inversion charges curves as a function of front gate bias evaluated using both of these approaches. Indeed, owing to the coupling term the front and back inversion charge densities are misestimated using the gate charge densities. In this comparative study we adhere to the second approach as it includes the coupling term conjecture. This choice is due principally to the fact that the coupling term will be included in the analytical compact model starting set of equations as well and secondarily to the fact that either the numerical or the analytical compact models do not provide the profile of electrons throughout the silicon film.

Numerical charge model validation by comparison to PS results:

Undoubtedly, a comparison stage where the results given by the numerical model are exposed to the ones provided by PS simulations in order to validate the former is vital in such study. Thusly, in this section we explore such comparison regarding the various electrostatic quantities. The comparison is performed for two different silicon thicknesses 𝑡 𝑠𝑖 = 7 𝑛𝑚 and 𝑡 𝑠𝑖 = 10 𝑛𝑚 in order to establish the numerical robustness of the numerical model regarding various geometrical configurations.

Decidedly, in the list of Figures, Figure 8 to Figure 12, we present the front surface potential, back surface potential, front inversion charge, back inversion charge, and gate-to-channel capacitance curves as functions of front gate voltage obtained by numerical computation is exposed to the curves provided by PS simulations.

Such Figures show an overall agreement between the numerical results and PS simulation results

, further underlying the pertinence of our numerical model. Some inaccuracies with respect to PS simulation results are exhibited by the numerical model mostly in the moderate and strong inversion regions, which can be attributed to the series of approximations that had to be done in order to build such model, in addition to the smoothing functions and parameters that can play a role in such inaccuracies as well.

The Current numerical model:

We describe in this part the long channel drain current calculation, assuming an infinite saturation velocity and in the frame of the gradual channel approximation. Such calculations are computed for the front and back channels separately based on the prior separation of inversion charges, then the total current will be straightforwardly the sum of the front and back currents.

The drift-diffusion numerical model:

Similarly to [14], we start from the general drain current equation expressed in Eq 4. 19. The mobility located originally inside the integral can have two forms in our study. Firstly, we consider an effective constant mobility 𝜇 𝑒𝑓𝑓 averaged throughout the inversion charge and the channel length, in which case the mobility term 𝜇 𝑒𝑓𝑓 can be displaced outside the integral. Ensuing, we consider the bell-shaped mobility law exhibited in Chapter 2. Using Eq 4.6 the closed-from expression of the term 𝑑𝜙 𝑖𝑚 𝑑𝑦 ⁄ can be derived, as in Eq 4. 20.

𝐼

𝑑𝜙 𝑖𝑚 𝑑𝑦 = 𝑑𝑉 𝑠1,2 𝑑𝑦 - 𝑑𝛥𝑉(𝑄 𝑔1,2 ) 𝑑𝑦 - 𝑘𝑇 𝑞 𝑁 2𝐷 1 1 -𝑒 -𝑄 𝑖𝑛𝑣1,2 𝑞 𝑁 2𝐷 𝑑𝑄 𝑖𝑛𝑣1,2 𝑑𝑦 Eq 4.20
Hereupon, we identify the term 𝜒 𝑠 1,2 = 𝑉 𝑠 1,2 -∆𝑉(𝑄 𝑔1,2 ) representing the potential of the corresponding subband, doing so, the drift and diffusion terms are uncovered, as in Eq 4. 21.

𝑑𝜙 𝑖𝑚 𝑑𝑦 = 𝑑𝜒 𝑠 1,2 𝑑𝑦 - 𝑘𝑇 𝑞 𝑁 2𝐷 1 1 -𝑒 -𝑄 𝑖𝑛𝑣1,2 𝑞 𝑁 2𝐷 𝑑𝑄 𝑖𝑛𝑣1,2 𝑑𝑦 Eq 4.21
The next step will be to inject Eq 4.21 in the general drain current equation Eq 4.19, followed by switching the integral boundaries initially defined along the electrical length of the channel into the corresponding evaluations of 𝜒 𝑠 For both the effective mobility and the mobility law conditions a numerical algorithm is made to perform the numerical calculations. Such algorithm considers the one-dimensional numerical integration along a quasi-Fermi level vector that increases monotonously along the channel with uniform segments, starting from the source point where 𝜙 𝑖𝑚 = 0, to the drain point where 𝜙 𝑖𝑚 = 𝑉 𝑑 . In fact, the subband potential and the inversion charge that correspondingly vary gradually along the channel, are computed in each step 𝑗.

Regarding the drift component of the current, the computation is performed using the trapezoidal rule, where for each quasi-Fermi level segment Nevertheless, the integral concerning the diffusion component is performed with respect to the dimensionless quantity 𝑢 1,2 = 𝑄 𝑖𝑛𝑣1,2 𝑞 𝑁 2𝐷 ⁄ , which is evaluated at the source and drain ends just as well, following the expression in Eq 4. 23. The integration is performed with the help of the quad function from the scipy.integrate package; such function uses a technique from the FORTRAN's library Quadpack. Since in our case we compute a finite integral, the integration is performed using a Clenshaw-Curtis method, which according to the official Scipy open source website it uses Chebyshev moments [12].

{ 𝐼 𝑑𝑟𝑖𝑓𝑡1,2 ≅ - 𝑊 𝐿 . 𝜇 𝑒𝑓𝑓 1,2 . ∑ (𝜒 𝑠1,2 𝑗+1 -𝜒 𝑠1,2 𝑗 ) 𝑗 𝑚𝑎𝑥 𝑗=0 . 𝑄 𝑖𝑛𝑣1,2 𝑗,𝑚𝑒𝑑 𝐼 𝑑𝑖𝑓𝑓1,2 = 𝑊 𝐿 . 𝑘𝑇𝑞 𝑁 2𝐷 . 𝜇 𝑒𝑓𝑓1,2 . 𝑘𝑇𝑞 𝑁 2𝐷 . ∫ 𝑢 1,2 1 -𝑒 -𝑢 1,2 𝑑𝑢 𝑢 1,2𝑑𝑟𝑛 𝑢 1,2𝑠𝑟𝑐 Eq 4.23
In the case of a bell-shaped mobility law, the expressions laid out in Eq 4. 23 As the introduction of the mobility law is peculiar to the present numerical model so far (such law is not incorporated in the PS simulations), we first present the computed transfer and output characteristics using the bell-shaped mobility law in this section, permitting the performance of comparison to PS results in the next one.

Figure 13 illustrates the computed transfer characteristics 𝐼 𝑑 = 𝑓(𝑉 𝑔1 ) in the linear regime, i.e. for 𝑉 𝑑 = 0.05 𝑉, and the saturated regime, i.e. for 𝑉 𝑑 = 1𝑉, as well as the output characteristics 𝐼 𝑑 = 𝑓(𝑉 𝑑 ), obtained by the numerical model employing the bell-shaped mobility law. Equivalently, the respective derivatives of the transfer and output characteristics obtained by the numerical model i.e. the conductance and transconductance curves are illustrated in Figure 14. Note the appealing shape of the transconductance in the linear regime, which resembles the one established by the gate-to-channel curves. In summary, we presently have a starting set of equations suitable for an accurate description of the electrostatic quantities in our system, such as the surface potentials, inversion charges, and gate-to-channel capacitance. The numerical robustness of this charge model is guaranteed by the employment of an engineered form of the quantum shift function which allows a better description of its behavior in the negative and positive gate charge regions, in addition to a smooth transitions between the two regions. The model is subsequently compared to the obtained results from PS simulations for two different silicon thicknesses in order to be validated. Justifiably, the built numerical charge model is suitable for all geometrical and backbias configurations.

Moreover, and on this basis, we developed a drift-diffusion numerical current model using the gradient of the quasi-Fermi level along the channel, and assuming initially a constant effective mobility, a configuration that allowed us to compare our results to the pre-obtained PS ones, ensuingly, the mobility bell-shaped law was employed in the model permitting its potential us to validate the analytical compact model results.

Indeed, the present chapter allowed us to establish a strong and consistent ground that we can use to build our analytical compact charge and current models in the next one.

This chapter will present the final stage of our study, in which we seek to develop an analytical charge and current models for FDSOI MOS transistors operating at cryogenic temperatures basing ourselves on the established numerical model in Chapter 4.

There are essentially three approaches to the compact modeling of MOS transistors [1]:

• The threshold-voltage-based approach which had its wide usage previously in models such as MOS Model 9, BSIM 3, and BSIM 4. • The charge-based approach, employed for instance in the EKV model • The surface-potential-based approach, employed for instance in the PSP and the L-UTSOI models. Such approach has become the conventional approach employed by the compact modelling community for both bulk and SOI devices.

Regarding the surface-potential-based approach, which is the approach we follow in the present work, the traditional scheme comprises the integration of the Poisson equation in the channel region, which with the help of the boundary conditions at the Si-SiO2 interface will generate the surface potential equation. Subsequently, the surface potentials are used to evaluate the terminal charges and the drain current, respectively [1].

Moreover, a good compact model must incorporate many attributes that have been elaborated meticulously in the renowned book within the compact modeling community [2]. Here, we cite for instance the ones we will aim to satisfy in the present work:

• The model should provide continuous results for currents and charges and their respective derivatives regarding each terminal voltage, guaranteeing thusly the numerical convergence of the model inside the circuits simulators where the non-linear Kirchhoff current law equations are solved commonly using the Newton-Raphson algorithm which requires smooth models to function effectively. • The equations involved in the formalism of the model as well as the model parameters need to be physically based as much as possible. • It should guarantee a smooth transition between the weak-, moderate-, and stronginversion operation regions. • It should meet the demands over large bias ranges.

• It should meet the demands over the temperature range of interest.

• The model should have as few parameters as possible, which should be linked to the device geometry and fabrication process parameters as strongly as possible. • It should be able to predict the behavior accurately using any combination of channel width and length values. • The model should be symmetric for symmetric devices; for instance, with this notion for an FDSOI device, we refer to the source-drain symmetry and the front-back interface symmetry. Whereas the structure of the device is not symmetric in the x direction (as front and back oxides does not have the same thickness), the model has to be flexible for an eventual switch of front-back gates usage.

• Finally, the model should be computationally efficient and numerically robust. Furthermore, the development of any analytical compact model passes by two big stages:

• The derivation of the core model: we mean by that the charge model as well as the current model for long channel devices • The implementation of various effects like add-ons or corrections. Those postimplemented effects can be very wide, such as small transistor effects, self-heating effect, the effect of doping, access resistances, parasitic capacitances, saturation velocity, channel length modulation, mechanical stress effects, impact ionization effect, gate leakage, inner fringe capacitance, GIDL… Note that, generally the core model remains unaltered, and the implementation of the effects happens gradually. In this respect, the compact model has to be improved in a continuous manner.

In this framework, we partition this chapter into two parts, a first part where we explain the buildout steps for a long channel model, starting by the elaboration of an initial guess for the surface potentials, from which the exact solution is engendered after a number of error correction steps, paving the way for a robust charge model. Subsequently, an analytical long channel drift-diffusion drain current is derived for both configurations, the constant effective mobility one and the bell-shaped mobility law one, taking into consideration the manifestation of the pinch-off points. Built on that, the incorporation of the short channel effects is reached in the second part. Such effects include the velocity saturation, the parasitic resistances, the DIBL, the threshold voltage roll-off, and the subthreshold slope degradation effects. Naturally, the results of the presented model are compared to experimental data as a final validation step in each of the two parts.

Note that, from now onwards, since the surface potentials 𝑉 𝑠1,2 are given the subscript 𝑠, we refer to the potential at the source end using the subscript 𝑠𝑟𝑐 as in 𝑉 𝑠𝑟𝑐 , and to the potential at the drain end using the subscript 𝑑𝑟𝑛 as in 𝑉 𝑑𝑟𝑛 . Equivalently, the front and back inversion charges at the source and drain terminals will be referred to as 𝑄 𝑖𝑛𝑣1,2 𝑠𝑟𝑐 and 𝑄 𝑖𝑛𝑣1,2 𝑑𝑟𝑟𝑛 respectively.

The Charge analytical model: 1.The surface potential initial guess derivation:

In this section we detail the necessary steps to develop a suitable analytical formulation of the front and back surface potentials. To this end, we start from our system of coupled equations recalled below in Eq 5.1 and written in its spread format as in Eq 5.2. It should be noted that at this stage, we seek the initial guesses of surface potentials in its classical form implying that the front and back quantum shift functions are not involved in Eq 5.2.

{ 𝑄 𝑔1 = -𝑄 𝑖𝑛𝑣1 + 𝑄 𝑐𝑝𝑙 𝑄 𝑔2 = -𝑄 𝑖𝑛𝑣2 -𝑄 𝑐𝑝𝑙 Eq 5.1 { 𝐶 𝑜𝑥1 . (𝑉 𝑔1 -𝑉 𝑓𝑏1 -𝑉 𝑠1 ) = 𝑞𝑘𝑇𝐴 2𝑑 𝑙𝑛 (1 + 𝑒𝑥𝑝 ( 𝑉 𝑆1 -𝑉 0 -𝜙 𝑖𝑚 𝑘𝑇 )) + 𝐶 𝑠𝑖 . (𝑉 𝑠1 -𝑉 𝑠2 ) 𝐶 𝑜𝑥2 . (𝑉 𝑔2 -𝑉 𝑓𝑏2 -𝑉 𝑠2 ) = 𝑞𝑘𝑇𝐴 2𝑑 𝑙𝑛 (1 + 𝑒𝑥𝑝 ( 𝑉 𝑆2 -𝑉 0 -𝜙 𝑖𝑚 𝑘𝑇 )) -𝐶 𝑠𝑖 . (𝑉 𝑠1 -𝑉 𝑠2 ) Eq 5.2
As indicated formerly, using the Fermi-Dirac statistics inherent to cryogenic consideration have the advantage of the explicit formulation in both strong and weak inversions, since in the weak inversion mode the terms labeling the inversion charge densities are abolished, leading to the expressions illustrated in Eq 5.3:

𝐶 𝑜𝑥1,2 . (𝑉 𝑔1,2 -𝑉 𝑓𝑏1,2 -𝑉 𝑠1,2 ) = ±𝐶 𝑠𝑖 . (𝑉 𝑠1 -𝑉 𝑠2 ) Eq 5.3
Correspondingly, in the strong inversion mode the asymptotic behavior of these terms can be depicted using the following pattern:

𝑞𝑘𝑇𝐴 2𝑑 𝑙𝑛 (1 + 𝑒𝑥𝑝 ( 𝑉 𝑆1,2 -𝑉 0 -𝜙 𝑖𝑚 𝑘𝑇 )) → 𝑞𝐴 2𝑑 (𝑉 𝑆1,2 -𝑉 0 -𝜙 𝑖𝑚 )
Eq 5.4

Yielding to the expressions illustrated in Eq 5.5, note that presently both mathematical expressions are composed of linear terms exclusively. 𝐶 𝑜𝑥1,2 . (𝑉 𝑔1,2 -𝑉 𝑓𝑏1,2 -𝑉 𝑠1,2 ) = 𝑞𝐴 2𝑑 (𝑉 𝑆1,2 -𝑉 0 -𝜙 𝑖𝑚 ) ± 𝐶 𝑠𝑖 . (𝑉 𝑠1 -𝑉 𝑠2 )

Eq 5.5

Moreover, and accounting for the dual channel operation, each of the front and back interfaces can be in the weak or strong inversion mode, thusly, four combinations are conceivable depending on the applied bias on each gate, as detailed below:

• 1 st case: when the front interface is in the weak inversion mode, and the back interface is in the weak inversion mode as well, this configuration will be entitled simply as the "weak-weak" configuration and given the subscript "ww". • 2 nd case: when the front interface is in the strong inversion mode, and the back interface is in the weak inversion mode, this configuration will be entitled simply as the "strongweak" configuration and given the subscript "sw". • 3 rd case: when the front interface is in the weak inversion mode, and the back interface is in the strong inversion mode, this configuration will be entitled simply as the "weakstrong" configuration and given the subscript "ws". • 4 th case: when the front interface is in the strong inversion mode, and the back interface is in the strong inversion mode as well, this configuration will be entitled simply as the "strong-strong" configuration and given the subscript "ss".

Consonantly, the mathematical statements describing each of the aforementioned asymptotic cases are:

• In the weak inversion / weak inversion case:

{ 𝐶 𝑜𝑥1 . (𝑉 𝑔1 -𝑉 𝑓𝑏1 -𝑉 𝑠1 ) = 𝐶 𝑠𝑖 . (𝑉 𝑠1 -𝑉 𝑠2 ) 𝐶 𝑜𝑥2 . (𝑉 𝑔2 -𝑉 𝑓𝑏2 -𝑉 𝑠2 ) = -𝐶 𝑠𝑖 . (𝑉 𝑠1 -𝑉 𝑠2 )
Eq 5.6

• In the strong inversion / weak inversion case:

{ 𝐶 𝑜𝑥1 . (𝑉 𝑔1 -𝑉 𝑓𝑏1 -𝑉 𝑠1 ) = 𝑞𝐴 2𝑑 (𝑉 𝑆1 -𝑉 0 -𝜙 𝑖𝑚 ) + 𝐶 𝑠𝑖 . (𝑉 𝑠1 -𝑉 𝑠2 ) 𝐶 𝑜𝑥2 . (𝑉 𝑔2 -𝑉 𝑓𝑏2 -𝑉 𝑠2 ) = -𝐶 𝑠𝑖 . (𝑉 𝑠1 -𝑉 𝑠2 )
Eq 5.7

• In the weak inversion / strong inversion case:

{ 𝐶 𝑜𝑥1 . (𝑉 𝑔1 -𝑉 𝑓𝑏1 -𝑉 𝑠1 ) = 𝐶 𝑠𝑖 . (𝑉 𝑠1 -𝑉 𝑠2 ) 𝐶 𝑜𝑥2 . (𝑉 𝑔2 -𝑉 𝑓𝑏2 -𝑉 𝑠2 ) = 𝑞𝐴 2𝑑 (𝑉 𝑆2 -𝑉 0 -𝜙 𝑖𝑚 ) -𝐶 𝑠𝑖 . (𝑉 𝑠1 -𝑉 𝑠2 )
Eq 5.8

• In the strong inversion / strong inversion case:

{ 𝐶 𝑜𝑥1 . (𝑉 𝑔1 -𝑉 𝑓𝑏1 -𝑉 𝑠1 ) = 𝑞𝐴 2𝑑 (𝑉 𝑆1 -𝑉 0 -𝜙 𝑖𝑚 ) + 𝐶 𝑠𝑖 . (𝑉 𝑠1 -𝑉 𝑠2 ) 𝐶 𝑜𝑥2 . (𝑉 𝑔2 -𝑉 𝑓𝑏2 -𝑉 𝑠2 ) = 𝑞𝐴 2𝑑 (𝑉 𝑆2 -𝑉 0 -𝜙 𝑖𝑚 ) -𝐶 𝑠𝑖 . (𝑉 𝑠1 -𝑉 𝑠2 )
Eq 5.9

For the purpose of obtaining the analytical expressions of the front and back surface potential initial guesses, labeled 𝑉 𝑠1 0 and 𝑉 𝑠2 0 respectively, we first need to write all the asymptotic cases in a more compact form. In order to do so, we will have to define two auxiliary terms 𝐶 𝑜𝑥1,2 * and 𝑉 𝑔1,2 * , which will have different expressions depending on whether the corresponding interface is in weak or strong inversion, following Eq 5.10 in the case where the front or back interface is in weak inversion, and Eq 5.11 in the case where the front or back interface is in strong inversion.

{ 𝐶 𝑜𝑥1,2 * = 𝐶 𝑜𝑥1,2 𝑉 𝑔1,2 * = 𝑉 𝑔1,2 -𝑉 𝑓𝑏1 Eq 5.10 { 𝐶 𝑜𝑥1,2 * = 𝐶 𝑜𝑥1,2 + 𝑞𝐴 2𝑑 𝑉 𝑔1,2 * = 𝐶 𝑜𝑥1,2 (𝑉 𝑔1,2 -𝑉 𝑓𝑏1 ) + 𝑞𝐴 2𝑑 (𝑉 0 + 𝜙 𝑖𝑚 ) 𝐶 𝑜𝑥1,2
Such arrangement will result in the direct analytical compact expressions of 𝑉 𝑠1 0 and 𝑉 𝑠2 0 , as depicted in Eq 5.13:

{ 𝑉 𝑠1 0 = ( 1 𝐶 𝑠𝑖 + 1 𝐶 𝑜𝑥2 * ) 𝑉 𝑔1 * + 1 𝐶 𝑜𝑥1 * 𝑉 𝑔2 * 1 𝐶 𝑠𝑖 + 1 𝐶 𝑜𝑥2 * + 1 𝐶 𝑜𝑥1 * 𝑉 𝑠2 0 = ( 1 𝐶 𝑠𝑖 + 1 𝐶 𝑜𝑥1 * ) 𝑉 𝑔2 * + 1 𝐶 𝑜𝑥2 * 𝑉 𝑔1 * 1 𝐶 𝑠𝑖 + 1 𝐶 𝑜𝑥1 * + 1 𝐶 𝑜𝑥2 * Eq 5.13
At this stage, we have obtained four surface potential expressions corresponding to the four asymptotic cases described previously. Correspondingly, the initial guess is obtained as a smooth minimum function of the four asymptotic cases described in Eq 5. 13.

Thereafter, a sequence of error corrections is applied to the surface potential terms in order to converge to the exact solutions. The iterative scheme presented in Eq 5.14 permits to get a new value of 𝑉 𝑠1 and 𝑉 𝑠2 each time we apply the corrections using 𝜀 1 (𝑉 𝑠1 ) and 𝜀 2 (𝑉 𝑠2 ). In our case, such operation needs to be done four times before the convergence to the exact values is achieved. Note that for the very first calculation the value of the terms 𝑉 𝑠1 𝑜𝑙𝑑 and 𝑉 𝑠2 𝑜𝑙𝑑 is plainly the initial guesses 𝑉 𝑠1 0 and 𝑉 𝑠2 0 respectively. Moreover, it should be noted that such iterative procedure does not use any loops, but rather its analytical expressions is repeated four times until convergence to the exact solution is achieved.

{ 𝑉 𝑠1 𝑛𝑒𝑤 = 𝑉 𝑠1 𝑜𝑙𝑑 + 𝜀 1 (𝑉 𝑠1 𝑜𝑙𝑑 ) 𝑉 𝑠2 𝑛𝑒𝑤 = 𝑉 𝑠2 𝑜𝑙𝑑 + 𝜀 2 (𝑉 𝑠2 𝑜𝑙𝑑 )
Eq 5.14

Decidedly, our next quest is to derive the analytical expressions of the errors 𝜀 1 and 𝜀 2 . In order to do so, we reinstate the quantum shift functions within the inversion charge density terms, into which a First-order Taylor expansion is applied. Correspondingly, we replace all the 𝑉 𝑠1 and 𝑉 𝑠2 terms in our system of coupled equations by the relationship 𝑉 𝑠1,2 = 𝑉 𝑠1,2 0 + 𝜀 1,2 , yielding to the terms 𝑄 𝑖𝑛𝑣1,2 (𝑉 𝑠1,2 0 + 𝜀 1,2 ) into which the first-order Taylor expansion is applied in the vicinity of 𝜀 1,2 = 0 and with the assumption that the error terms 𝜀 1,2 are infinitesimal quantities relatively to the 𝑉 𝑠1,2 terms, as depicted in :

𝑄 𝑖𝑛𝑣1,2 (𝑉 𝑠1,2 0 + 𝜀 1,2 ) = 𝑄 𝑖𝑛𝑣1,2 (𝑉 𝑠1,2 0 ) + 𝜕𝑄 𝑖𝑛𝑣1 𝜕𝑉 𝑆1 (𝑉 𝑠1,2 0 ). 𝜀 1,2
Eq 5.15

d∆V(Q g1,2 ) dQ g1,2 = (1 γ ⁄ ) 0.5 (1 - Q g1,2 Q pos (Q g1,2 ) -Q neg (Q g1,2 ) ) + 2 3 . β QM . 0.5 (1 + Q g1,2 Q pos (Q g1,2 ) -Q neg (Q g1,2 ) ) . 1 Q pos 1/3 (Q g1,2 )
Eq 5.27

Finally, using the definitions depicted in Eq 5.19 and Eq 5.18 and the expressions presented in Eq 5.17, we can get to the compact closed-form of the error corrections 𝜀 1,2 as illustrated in Eq 5. 28.

{ 𝜀 1 = 𝐶 𝑠𝑖 2 𝑅 2 + 𝐶 𝑅 2 𝑅 1 𝐶 𝑅1 𝐶 𝑅2 -𝐶 𝑠𝑖 2 𝐶 𝑠𝑖 1 𝜀 2 = 𝐶 𝑠𝑖 1 𝑅 1 + 𝐶 𝑅 1 𝑅 2 𝐶 𝑅2 𝐶 𝑅1 -𝐶 𝑠𝑖 1 𝐶 𝑠𝑖 2
Eq 5. 28 Accordingly, the kinetics of the convergence of the front/back surface potentials starting from the initial guesses and reaching the exact numerical solution is depicted in Figure 1.Note that, the small segment of the curves that does not change along with iterations is attributed to the δQ g parameter that becomes significant around the zero-gate-charge point.

Analytical charge model validation by comparison to numerical results:

As per usual, before proceeding further, our analytical charge model needs to be validated first by confronting its results, whether it is for the surface potentials or the inversion charge densities, to the exact solution obtained through numerical calculations.

In Figure 2 we expose the front and back surface potential curves as a function of front gate voltage and for different back biases obtained by analytical calculations to the numerical solution. Equivalently, in Figure 3 we expose the front and back surface potential densities as functions of front gate voltage and for different back biases obtained by analytical calculations to the numerical solution. Analogously, Figure 4 demonstrates the good agreement between the numerical 𝐶 𝑔𝑐 curves (solid red lines) and the ones obtained through analytical calculations (blue dashed curves) for three different back biases. Clearly, we have a notable good agreement between the curves analytically computed and the ones obtained through numerical calculations for all back bias configurations, which validates our analytical charge model and paves the way for us to go to the next section.

The analytical drain current model:

In this section, we describe the long channel drain current calculation within an infinite saturation velocity assumption. Equivalently to the procedure detailed in Chapter 4, the general form of a drain current equation is first recalled by Eq 5. Prior to the derivation of the analytical expression for the diffusion and drift integrals, one additional information is essential and needs to be discussed here, the so called "pinchoff/saturation" point. Expressly, when the drain bias is relatively small, i.e. we reside in the linear region, the inversion charge density at the drain end of the channel is moderately lower than the inversion charge density at the source end. As we increase the drain bias (for a fixed gate bias), the current increases until it reaches its maximum value, the drain current saturation limit 𝐼 𝑑,𝑠𝑎𝑡 , but the inversion charge density at the drain side decreases until finally it vanishes when the applied drain bias reaches the value 𝑉 𝑑,𝑠𝑎𝑡 . In other words, the surface potential saturates at the drain end of the channel when the drain current saturation occurs. This phenomenon is called pinch-off and is illustrated in Figure 5 [3]. When 𝑉 𝑑 increases beyond saturation, the pinch-off point moves toward the source, but the drain current remains essentially the same. This is because for 𝑉 𝑑 > 𝑉 𝑑,𝑠𝑎𝑡 the voltage at the pinch-off point remains at 𝑉 𝑑,𝑠𝑎𝑡 and the current will always follow the equivalence illustrated by Eq 5.30 [3]: Note that, according to the generic current equation and in order to maintain the current continuity throughout the channel, the decrease of the inversion charge 𝑄 𝑖𝑛𝑣 at the drain side must be accompanied by an increase of the term 𝑑𝜙 𝑖𝑚 𝑑𝑦 ⁄ . When 𝑉 𝑑 reaches 𝑉 𝑑,𝑠𝑎𝑡 , we have 𝑄 𝑖𝑛𝑣 = 0 and correspondingly 𝑑𝜙 𝑖𝑚 𝑑𝑦 ⁄ = ∞. This implies that the electric field in the y direction changes more rapidly than the field in the x direction, and the gradual channel approximation breaks down in this region [3]. Extensively, beyond the pinch-off point electrons are no longer confined to the surface channel, and a two-dimensional analysis of the device is necessary for the region between the pinch-off point and the drain point [3].

∫
Expressly, one should point out that since for the dual channel operation we have two pinchoff points, the designated point considered hereafter in the calculation of 𝑄 inv1,2 𝑠𝑎𝑡 and 𝜒 s1,2 𝑠𝑎𝑡 corresponds to the pinch-off voltage of the strongest channel.

The diffusion current component computation:

Regarding the diffusion term integral the mobility is maintained constant. Such choice is justified in our study as we already established in Section 1.2.2. Therefore, the hereabouts derivation of the diffusion current component will be implemented in either derivation of drain current, the one with effective mobility and the one with the mobility function. Firstly, we recall our established formula of the diffusion integral: Considering the indefinite nature of the above-stated integral, the best approach is to approximate it by a definite one, to which an analytical solution exists, following: At this point we need to implement an additional approximation that allows us to simplify the inversion charge expression involved in the drift current integral. Expressly, to realize an analytical calculation of the drift integral described in Eq 5.34, we need to analyze the dependence of -𝑄 𝑖𝑛𝑣 with 𝜒 𝑠 , as detailed in the next section.

𝐼 𝑑𝑖𝑓𝑓1,
∫ 𝑢 1,2 1 -𝑒 -𝑢 1,2 𝑑𝑢 𝑢 1,2𝑑𝑟𝑛 𝑢 1,2𝑠𝑟𝑐 ≅ ∫ (𝑢 + 𝑒 -𝑢 𝑟 ) 𝑑𝑢 𝑢 1,2𝑑𝑟𝑛 𝑢 1,2𝑠𝑟𝑐 = 𝑢 1,2𝑑𝑟𝑛 2 2 - 𝑢 1,2𝑠𝑟𝑐 2

The inversion charge linearization technique:

Figure 7 illustrates the dependence of inversion charge densities -𝑄 𝑖𝑛𝑣1,2 with the corresponding subband potential 𝜒 𝑠1,2 for three different back biases. We can see from Figure 7 that the dependence is typically linear in the case of null and negative back biases. Whereas, due to the dual channel operation in the case of positive back bias configuration it presents two distinguishable slopes behavior. Generally, the conventional methodical procedure to address the drift current computation challenge in literature is the inversion charge linearization technique. Such technique is referred to in [4] as the "Symmetric linearization method" and its comprehensive form is used in the formulation of the SP and PSP models.

According to [5], the inversion charge linearization technique approximates the inversion charge by its First-order Taylor expansion in the vicinity of the potential middle point, as exhibited in Eq 5.35: 𝑉 𝑠𝑚 is the surface potential midpoint that allows the approximation of the inversion charge by a one-slope line throughout the quasi-Fermi level. Note that the linear inversion charge approximation is appropriate to the analytical model in addition to the already taken approximation for the numerical model development, implying a relative error between the analytical drain current and the numerical computed one, as we will see later.

𝑄 𝑖𝑛𝑣1,2 = 𝑄 𝑖𝑛𝑣1,2 𝑚 + 𝛼
In contrast to the aforementioned works [4], [5] where the linearization of the inversion charge is performed with respect to the surface potentials 𝑉 𝑠1,2 , in our case such linearization is performed with respect to the subband potentials 𝜒 𝑠1,2 .

Moreover, in the case of double channel operation, due to the explicit two-slope behavior in the positive back configuration, we cannot use the mere all-over linearization of the inversion charge to obtain a drain current expression that is valid in all regimes. Instead, and correspondingly to the L-UTSOI model [6], we compute the slope of the -𝑄 𝑖𝑛𝑣 (𝜒 𝑠 ) curves in the source and saturation points respectively, then we approximate -𝑄 𝑖𝑛𝑣 for the integral calculation as the maximum value of two linear functions of 𝜒 𝑠 , named -𝑄 𝑖𝑛𝑣 𝑙𝑖𝑛 and -𝑄 𝑖𝑛𝑣 𝑙𝑖𝑛 ̃ and defined respectively in the source and saturation ends by the expression in Eq 5.37. Therefore, in order to perform the drift current piece-wise calculation the value of the subband potential at the intersection point needs to be established. The closed-form of 𝜒 s1,2 𝑖𝑛𝑡 can be derived directly by satisfying the equality of the two expressions in Eq 5.37 since such point corresponds to the intersection of the two slopes. Such equality which will convey us the expression in Eq 5.40: Therefore, the calculation of the drift component of the drain current consists of the following sequence:

{ -𝑄

𝜒 s1,2 𝑖𝑛𝑡 = (𝑄
• Computation of the subband potentials and inversion charge densities at the source end for both the front and back interfaces. • Computation of the subband potentials and inversion charge densities at pinch-off end for both the front and back interfaces. • Computation of the inversion charges derivative 𝜆 1,2 (the charge linearization slopes) at the source and saturation points respectively. • Computation of the drift current using Eq 5.39.

Indeed, as can be seen in the aforementioned procedure sequence, this approach requires the computation of the derivative of -𝑄 𝑖𝑛𝑣 with respect to 𝜒 𝑠 at both the source and drain sides, such derivation will be detailed in the next section.

Derivation of the inversion charge linearization slopes:

In this part, our focal point would be to obtain the closed-form expressions of the inversion charge linearization slopes, following our approach of computing the slope at the source and saturation points respectively. It should be noted that such derivation is established for a fixed front and back biasing values, leaving 𝜙 𝑖𝑚 as the only variable (variating thusly 𝜒 𝑠1,2 ) throughout the channel, i.e. what we compute factually is the term

𝑑𝑄 𝑖𝑛𝑣1,2 𝑑𝜙 𝑖𝑚 ⁄ 𝑑𝜒 𝑠1,2 𝑑𝜙 𝑖𝑚 ⁄
.

Foremostly, and for the sake of preserving the compact form of our equations at will, we designate the arguments of the exponential function, originally allocated to the front and back interfaces separately, as simply 𝐴𝑟𝑔 1 and 𝐴𝑟𝑔 2 following: Which can be rearranged to be written as: We call the terms residing between the brackets 𝐶 1 and 𝐶 2 respectively:

𝐴𝑟𝑔 1,2 = 𝑉 𝑠 1,2
𝑑𝑄 𝑖𝑛𝑣1,2 𝑑𝑉 𝑠1,2 = 𝑞 𝐴 2𝐷 (1 - 𝑑∆𝑉(𝑄 𝑔1,2 ) 𝑑𝑉 𝑠1,2 - 1 𝑑𝑉 𝑠1,2 𝑑𝜙 𝑖𝑚 ⁄ ) 𝑒 𝐴𝑟𝑔 1,2 1 + 𝑒 𝐴𝑟𝑔
{ 𝐶 1 = C ox1 + C si + 𝑞 𝐴 2𝐷 (1 - 𝑑∆𝑉(𝑄 𝑔1 ) dV s1 ) 𝑒 𝐴𝑟𝑔 1 1 + 𝑒 𝐴𝑟𝑔 1 𝐶 2 = C ox2 + C si + 𝑞 𝐴 2𝐷 (1 - 𝑑∆𝑉(𝑄 𝑔2 ) dV s2 ) 𝑒 𝐴𝑟𝑔 2 1 + 𝑒 𝐴𝑟𝑔 2
Eq 5.48

And the terms representing activation functions as 𝑠 1 and 𝑠 2 respectively:

{ 𝑠 1 = 𝑞 𝐴 2𝐷 𝑒 𝐴𝑟𝑔 1 1 + 𝑒 𝐴𝑟𝑔 1 𝑠 2 = 𝑞 𝐴 2𝐷 𝑒 𝐴𝑟𝑔 2 1 + 𝑒 𝐴𝑟𝑔 2
Eq 5.49

We thusly can write our system of equations in the next form: We thereby have attained the closed-form expression of both the terms 𝑑𝑉 𝑠1 𝑑𝜙 𝑖𝑚 ⁄ and 𝑑𝑉 𝑠2 𝑑𝜙 𝑖𝑚 ⁄ . Finally, the closed-form expression of the linearization slopes can be directly expressed following Eq 5.53:

{ 𝐶 1 𝑑𝑉 𝑠1 𝑑𝜙 𝑖𝑚 = 𝑠 1 + C si . 𝑑𝑉 𝑠2 𝑑𝜙 𝑖𝑚 (1)
{ 𝑑𝑄 𝑖𝑛𝑣1 𝑑𝑉 𝑠1 = 𝑠 1 (1 - 𝑑∆𝑉(𝑄 𝑔1 ) 𝑑𝑉 𝑠1 - 𝐶 1 𝐶 2 -C si 2 𝑠 1 𝐶 2 + C si 𝑠 2 ) 𝑑𝑄 𝑖𝑛𝑣2 𝑑𝑉 𝑠2 = 𝑠 2 (1 - 𝑑∆𝑉(𝑄 𝑔2 ) 𝑑𝑉 𝑠2 - 𝐶 1 𝐶 2 -C si 2 𝑠 2 𝐶 1 + C si 𝑠 1 )
Eq 5.53

Employing the closed-form expressions of the linearization slopes and the linear form of the inversion charge densities presented in Eq 5.37, the front and back inversion charges linearization is verified in Figure 8. We can see that the applied expressions describe adequately the inversion charge throughout the subband potential range and for all back biases.

saturation point for the second part of integration, are supplanted using the sequence from the source point to the first effective point for the first part of integration, and from the first effective point to the second effective point for the second part of integration.

Correspondingly, and considering a constant effective mobility along the channel, the final analytical expression of the drift current can be directly derived, yielding Eq 5.56: Considering on the other hand the bell shape mobility law described in Eq 5.57, the analytical expression of the drift current with a mobility function will be as in Eq 5.58: Finally, it should be noted that based on the discussion performed in Chapter 2, we include also the effective temperature using the expression illustrated by Eq 5.59:

𝜇 𝑛 (𝑄 𝑖𝑛𝑣
𝑇 𝑒𝑓𝑓 (𝑇, 𝑇 0 ) = 𝑇 0 . (1 + 𝛼. 𝑙𝑛 (1 + exp ( 𝑇 -𝑇 0 𝛼. 𝑇 0 )))
Eq 5.59

Thusly, we have discussed the different aspects concerning the long channel transistor model development, allowing us to address the short channel effects that would be addressed as addons to the core model as we will see in the next section.

Short channel MOSFET current model:

Certain physical phenomena are negligible in large dimensions devices but become more significant in determining the behavior of the MOSFET in the case of devices with reduced dimensions. As the transistor channel length is reduced, the electrostatic control of the source and drain zones increases and preponderate that of the gates. Such 2-D electrostatic effects generate a degradation of the transistor subthreshold slope, a linear threshold voltage decrease, also known as the roll-off effect, and an increased sensitivity of the threshold voltage to the drain bias, also known as the DIBL effect.

All these aforementioned effects can be modeled by considering the transistor as a simple capacitive network, such capacitive modeling is true in the weak inversion regime, then, this approach was extended as an approximative one into the strong inversion regime as well, such choice is justified by our need for fast computation for the consideration of the 2-D analysis, which necessitates numerical computations otherwise.

Note that, in order to preserve the consistency of our model, the aim is to always keep the long channel model (the core model) and implement the short channel effects as suitable modifications. Thusly, the so-far designated short channel effects are introduced by modifying the device geometry and applied biases according prior to surface potential calculations as it will be demonstrated afterwards.

Moreover, when we apply an electric field upon the silicon channel, the mobile carriers are accelerated and gain a drift velocity that overlays their random thermal motion [3]. Note that, such velocity of the electrons does not increase indefinitely under field acceleration, since they are scattered repeatedly and lose their gained energy after each inelastic collision. At low electric fields, the drift velocity 𝑣 𝑑 is proportional to the electric field strength 𝐸, with the mobility 𝜇 as the proportionality constant 𝑣 𝑑 = 𝜇. 𝐸, where the mobility is proportional to the time interval between collisions and inversely proportional to the effective mass of electrons [3].

Nonetheless, note that the above-stated linear velocity-field relationship is valid only when the electric field is not too high. For at high fields, the average carrier energy increases, and carriers lose their energy by optical phonon emission nearly as fast as they gain it from the field, engendering a decrease of the mobility as the field increases until eventually the drift velocity reaches a limiting value, such mechanism is called velocity saturation [3].

Implementation of velocity saturation effect:

Comprehensively, when the drain voltage increases in a long channel device, the drain current first increase, then becomes saturated at a voltage equals to 𝑉 𝑑,𝑠𝑎𝑡 with the onset of the pinchoff at the drain side.

Comparatively, in a short channel device, the saturation of drain current may occur at a much lower voltage due to this velocity saturation effect, implying a saturation current 𝐼 𝑑,𝑠𝑎𝑡 that is detached from the 1/𝐿 dependence depicted in Eq 5.60 and Eq 5.61 for long channel devices. In other words, the drain current saturates due to either pinch-off or velocity saturation at the drain, the reported value of the velocity saturation in literature is around 𝑣 𝑠𝑎𝑡 ≈ 7 -8 . 10 4 𝑚/𝑠.

According to [7], the drain current saturation can be derived through the velocity saturation by considering the maximum attainable value by the conductivity in the channel, as depicted by Eq 5.60: 𝜎 𝑚𝑎𝑥1,2 = 𝑊. 𝐶 𝑜𝑥1 . 𝑣 𝑠𝑎𝑡 Eq 5.60

In the present work, we choose to carry out an equivalent description to define an upper limit of the drain current as featured in Eq 5.61, where the designated charge density is the uttermost value in the silicon channel i.e. the one at the source side.

𝐼 𝑑1,2 𝑠𝑢𝑝 = 𝑊. 𝑄 𝑠𝑟𝑐1,2 . 𝑣 𝑠𝑎𝑡 Eq 5.61

Such upper limit of the drain current 𝐼 1,2 𝑠𝑢𝑝 is then implemented into the model through the application of the Matthiessen rule between it and the initial drain current of the corresponding interface obtained from the long channel supposition, as in Eq 5.62:

𝐼 𝑑1,2 𝑓𝑖𝑛𝑎𝑙 = 𝐼 1,2 𝑠𝑢𝑝 . [1 + (𝐼 𝑑1,2 𝑙𝑜𝑛𝑔 𝐼 𝑑1,2 𝑠𝑢𝑝 ⁄ ) -𝑚 ] -1 𝑚 ⁄
Eq 5.62

Implementation of the DIBL and charge sharing effects:

In this part, for the sake of simplification, we consider our system of equations in the weak inversion configuration i.e., without the mobile charges. Considering firstly, the familiar 1-D system of equations:

{ 𝐶 𝑜𝑥1 . (𝑉 𝑔1 -𝑉 𝑓𝑏1 -𝑉 𝑠1 ) = 𝐶 𝑠𝑖 . (𝑉 𝑠1 -𝑉 𝑠2 ) 𝐶 𝑜𝑥2 . (𝑉 𝑔2 -𝑉 𝑓𝑏2 -𝑉 𝑠2 ) = 𝐶 𝑠𝑖 . (𝑉 𝑠2 -𝑉 𝑠1 )
Eq 5.63

After some rearrangement, we can write the 1-D system of equations in the format illustrated in Eq 5.64. This format of the system of equations will serve as a reference to which will be comparing the equivalent system of equations in when 2-D electrostatics are applied. If we consider the 2D system of equations on the other hand, the charge conservation principle in this case must be rather applied to the capacitive scheme presented in Figure 9. Where the correspondent charge conservation equations are given by:

{ C ox1 (V g1 -V fb1 -V s1 ) + [C s1 (V src + V 0 ) -C s1 V s1 ] + [C d1 (V d + V 0 ) -C d1 V s1 ] = C si (V s1 -V s2 ) C ox2 (V g2 -V fb2 -V s2 ) + [C s2 (V src + V 0 ) -C s2 V s2 ] + [C d2 (V d + V 0 ) -C d2 V s2 ] = C si (V s2 -V s1 )
Eq 5.65

Where within the supposition of a perfectly symmetric scheme, the four capacitances C s1 , C s2 , C d1 , and C d2 are identical, and adhere to the following definition:

C s1 = 𝜀 𝑠𝑖 𝑊 𝑡 𝑠𝑖 2 ⁄ 𝐿 2 ⁄ = 𝜀 𝑠𝑖 𝑊 𝑡 𝑠𝑖 𝐿 Eq 5.66
Thusly, if we compare the 2-D electrostatic short channel system of equations given by Eq 5.65 to its equivalent 1-D long channel one given by Eq 5.64, after nullifying the four capacitances C s1 , C s2 , C d1 , and C d2 , we recognize that the surface potentials can still be computed using the long channel approach, rather, the front and back oxide capacitances, and the front and back applied gate biases, are altered into new effective values given by Eq 5.67, where source bias is set as a reference i.e. V src = 0:

{ 𝐶 𝑜𝑥1,2 𝑒𝑓𝑓 = 𝐶 𝑜𝑥1,2 . (1 + C s1,2 + C d1,2 𝐶 𝑜𝑥1,2 ) 𝑉 𝑔1,2 𝑒𝑓𝑓 = 𝑉 𝑔1,2 + C d1,2 𝐶 𝑜𝑥1,2 . (V d + V 0 ) + C s1,2 𝐶 𝑜𝑥1,2 . V 0 (1 + C s1,2 + C d1,2 𝐶 𝑜𝑥1,2 )
Eq 5.67

Correspondingly, the implementation of these effects into our model is made through the use of these effective oxide capacitances and applied gate biases. Note that this modification must be included in the pinch-off point computation as well as this one depends on the corresponding values of 𝑉 𝑔1,2 and 𝐶 𝑜𝑥1,2 as well.

Note furthermore that, the method presented in this section will encompass the impact of the DIBL effect along with the subthreshold slope degradation, and the threshold voltage roll-down. Such collation of the impact of supposedly several effects is due to the fact that these effects have the same physical origin.

Implementation of the parasitic resistance effect:

Each discussion we have raised in the present study considered only the intrinsic part of the transistor which leads the behavior of the transistor; parasitic resistances on the other hand belong to the extrinsic part that connects the intrinsic part to other devices of an integrated circuit [2]. The impact of this extrinsic part can be seen in a general decreasing of the device performance for instance; such impact becomes palpable as the dimension of the device are reduced [2]. Several components form the extrinsic part, such as the parasitic resistances and the parasitic capacitances; the discussion of the parasitic resistances is alone considered in the present work.

In our approach, we will consider a simple scheme with two access parasitic resistances, one at the source side and the other at the drain side as in Figure 10, where the intrinsic part of our device is depicted by the rectangle in the middle. In this simplistic scheme we consider the effect of the parasitic resistances is delimited to the front gate oxide alone.

There are several components to the source and drain series resistances such as the contact resistance 𝑅 𝑐𝑜 , the current spreading resistance 𝑅 𝑠𝑝 that depends on the source/drain junction depth and the inversion layer thickness, the channel resistance 𝑅 𝑐ℎ , and the sheet resistance 𝑅 𝑠ℎ that is usually small in comparison to 𝑅 𝑐ℎ . In the next paragraph we will cover a direct method in order to implement the parasitic resistances effect in our model. Such method ranges from the linear regime for the implementation of the corresponding resistances, which will be then generalized to other regimes.

Therefore, starting from the linear regime one can assume that such introduction of parasitic resistances transforms the actual applied bias into altered ones as in (𝑉 𝑔𝑠 , 𝑉 𝑑𝑠 ) Where the drain current can be written as a function of the rest quantities as in:

𝐼 𝑑 = 𝐼 𝑑 0 1 + 𝑊 𝐿 𝜇 𝑛 (𝑄)𝑄 𝑖𝑛𝑣 2𝑅 𝑠 Eq 5.72
Where 𝐼 𝑑 0 = 𝑊 𝐿 𝜇 𝑛 (𝑄)𝑄 𝑖𝑛𝑣 𝑉 𝑑 is the obtained drain current prior to the implementation of the parasitic resistances.

Accordingly, Eq 5.72 can be rearranged to be written as in Eq 5.73, as this format can support the allegation of an altered mobility law. Correspondingly, the altered mobility law 𝜇 𝑛 (𝑄) ′ will be given by Eq 5.74:

𝜇 𝑛 (𝑄) ′ = 𝜇 𝑛 (𝑄) 1 + 𝑊 𝐿 𝜇 𝑛 (𝑄)𝑄 𝑖𝑛𝑣 2𝑅 𝑠 Eq 5.74
Such altered mobility law can be also written in the format given by Eq 5.75: Note that this approach is Comparable to the one used in both the PSP and the UTSOI1 models, where the parasitic resistance effect is implemented in via the inclusion of an additional term to the mobility degradation factors i.e., the Coulomb scattering and surface roughness scattering degradation factors in our case. Finally, note that 𝑅 𝑠 is an adjustable parameter, which equivariantly to L-UTSOI model, takes the typical values of the order 𝑅 𝑠 = 60Ω for a channel width of 1𝜇𝑚 in the present study. Note as well that recognizing that this is just a first order approximation which is not necessarily valid in high nonlinear regimes, we choose to adopt this 𝑅 𝑠 inclusion approach and extend it to the nonlinear regime of operation, as in PSP and L-UTSOI models for simplification.

𝜇 𝑛 (𝑄) ′ =

Comparison of the analytical model results to the experimental data:

In this section, we will be exposing the analytical model results to the collected experimental data for three different FDSOI structures with the same width 𝑊 = 1𝜇𝑚 and channel lengths of 𝐿 = 1𝜇𝑚, 𝐿 = 120 𝑛𝑚 and 𝐿 = 30 𝑛𝑚 respectively.

Expressly, the experimental data fitting in our case is essentially based on the choice of a set of the parameters 𝑉 𝑓𝑏1,2 , 𝜇 𝑚𝑎𝑥1,2 , 𝑄 𝑐1,2 for long channel transistors and the parameters 𝑉 𝑓𝑏1,2 , 𝜇 𝑚𝑎𝑥1,2 , 𝑄 𝑐1,2 , 𝑅 𝑠 , 𝑣 𝑠𝑎𝑡 , 𝑚 for short channel transistors. Outstandingly, in some cases as we will see thereafter, we will be bounded to an attribution of a set of parameters to each back bias configuration. Such customized treatment could have been redeemed by a standard parameter extraction phase; however such task would necessitate an extensive study for this aim, which is beyond the scope of the present work.

Regarding the exposition of the analytical model results to experimental data for long channel transistor, the first challenge that we will be facing is finding a way to emulate the intersubband scattering that appears in the linear 𝐼 𝑑 (𝑉 𝑔1 ) characteristics. Such challenge is not evident as such effect is ascribed to subbands interaction (as already illustrated) appearing during the transport phase, making it an appropriate quantum mechanical phenomenon. Engineering a compact model formula that perfectly describes such effect in a physics-based way is imperceptible and has never been addressed in literature. For such reasons, in this part we will be seeking the approach described in the next paragraph to emulate such effect.

As an overall analysis regarding the intersubband scattering, one can say that a big portion of the electrons that were originally in the back side of the film, swing to the front side of the film, and as the mobility in the back is improved compared to the one in the front (due to high-k oxide in the front), the charge that at first see the back improved mobility, see after the transition the front deteriorated mobility, leading to a decrease of the total mobility of the system, this is what explains the drain current decrease and the negative transconductance. Plainly, to explicit this effect in our model we considered that the maximum value of the mobility in the back channel is in itself dependent on the mean value of the front inversion charge, in a way that the more front inversion charge we have the less back mobility we have.

Thusly, to depict the signature of such mechanism through the alteration of the back mobility in the present work we use the expression presented in Eq 5.77:

𝜇 𝑚𝑎𝑥2 *
= 𝜇 𝑚𝑎𝑥1,2 . exp ( -𝑄 1 𝑚𝑒𝑑 𝑄 𝑟𝑒𝑓 )

Eq 5.77

Where 𝜇 𝑚𝑎𝑥2 * is the reduced back-channel mobility, 𝑄 1 𝑚𝑒𝑑 = 𝑄 1 𝑠𝑟𝑐 + 𝑄 1 𝑑𝑟𝑛 2 is the median value of front inversion charge, and 𝑄 𝑟𝑒𝑓𝑏 is a fitting parameter that controls the hump.

Note that, whereas the front and back flat band voltages, front and back upper values of the mobilities along with the inversion charge parameter 𝑄 𝑟𝑒𝑓 take values that are independent of the back bias, in this case 𝑉 𝑓𝑏1 = -0.09 𝑉, 𝑉 𝑓𝑏2 = 0.65 𝑉, 𝜇 𝑚𝑎𝑥1 = 0.95 𝑚 2 𝑉. 𝑠 ⁄ , 𝜇 𝑚𝑎𝑥2 = 4.84 𝑚 2 𝑉. 𝑠 ⁄ ,𝑄 𝑟𝑒𝑓 = 10.6 * 10 -4 𝐶/𝑚 2 ; 𝑄 𝑐1 , 𝑄 𝑐2 are treated as functions of the back-bias 𝑉 𝑔2 as depicted in the following plots: Note as well that, such typical approach of parameter attribution is unique for the linear transfer characteristics long channel transistor where the intersubband scattering effect is present and is conserved for the saturated transfer characteristics and the output characteristics with weaker ensuing accuracy. Indeed, the parameter values given in this last section are not claimed to be the ultimate ones, but rather such choice of parameter values is only meant to give us the opportunity to inspect the current outputs of our analytical model, offering a semi-quantitative analysis of our model results with respect to the experimental data. Further parameter optimization would require more work, going beyond the term of this thesis.

To summarize, we succeeded in the course of the present chapter to give a detailed step-by-step demonstration of our analytical model, starting from the derivation of the surface potentials initial guesses to the iterative process of error corrections yielding to the corresponding exact solutions. Accordingly, we exhibited an adapted two-slope method of the inversion charge linearization technique in order to be utilized in the drift current computation. Consistently, we derived the closed-form analytical solutions of both the diffusion and drift current integrals. Subsequently, the short channel effects were implemented as add-ons to the long channel model, allowing the exposition of its results to the collected experimental data for both the short and long channel transistors.

Conventionally, quantum computers consist of two parts: a quantum processor that comprises a set of qubits, and a classical electronic interface part required to perform the control and readout of quantum states. In order to acquire a more compact and more reliable system the classical logic interfaces needs to be established in the cryogenic chamber, resulting in the desirable aspects of an enhanced clock speed, an improved noise performance, a reduced signal latency/timing errors, and larger bandwidth. Based on the qubit sensitivity, the classical electronic interface needs to meet some precise conditions such as low temperature functionality, signal accuracy, high speed, low transmitted noise, and the generation of specific microwave bandwidths. In comparison to bulk CMOS, FDSOI transistors could offer the ideal cryogenic device performance already, making it a very good candidate for such task. However, process-design kits lack models describing the operation of MOS devices at cryogenic temperatures. Building compact models for cryogenic operation is crucial and needs to be tackled urgently. This is where the present thesis emerges as a venture to satisfy such need.

To tackle such task two approaches have been adopted by the research community. The first approach consists of adapting existing standard models for deep cryogenic operation. Such adaptation is made using empirical solutions in order to improve their accuracy and predictability in cryogenic operation. Such approach comes along with a lot of advantages, since these standard models contain already all the additional effects and features and are numerically robust. However, these advantages come along with the limitations that these standard models are not adapted for cryogenic operation, namely because they consider a 3D gas of electrons and use MB statistics to describe their distribution. The second approach consists of building a physics-based models aimed for cryogenic operation from scratch. Such choice allows to overpass the limitations of the first approach. For instance, 2-D electron gas along with the use of FD statistics can be considered initially. Nevertheless, these models are not as of now mature enough to be implemented in PDK's.

Addressing the task of building compact models suitable for cryogenic operation has to deal with two challenges, our lack of understanding for the physical phenomena that appears at those conditions and the numerical management of different mathematical expressions that describes such physical phenomena.

Accordingly, in the course of the present thesis, we discussed the Maxwell-Boltzmann approximation validity down to cryogenic temperatures via the exposition of the reasons we believe the MB approximation does not hold at cryogenic temperatures and must be traded by full Fermi-Dirac statistics. Indeed the choice of maintaining the MB approximation is applicable in some specific cases where the doping level beneath the degenerate limit and was preserved for numerical reasons considering that the implementation of Fermi-Dirac statistics would necessitate a numerical integration. Nonetheless, considering the relative position of the quasi-Fermi level which could traverse slightly the conduction band edge i.e. 𝐸 𝑓 ≥ 0, the MB approximation becomes inaccurate, and the electrons distribution is described by the Fermi-Dirac statistics.

Concerning the numerical integration argument, we demonstrated in the course of this thesis that the use of complete Fermi-Dirac statistics does not imply a numerical integration process due to the 2D subband systems, restraining thusly the development of an explicit model. If anything, using Fermi-Dirac statistics inherent to cryogenic consideration have the advantage of the explicit mathematical formulation in both strong and weak inversions.

In addition to the existence of a two-dimensional electron system allowing the use of the charge sheet approximation, we demonstrated the 2D subband is not a mere step function, but it exhibits a band tail of states, proposing thusly a suitable continuous expression that describes the exponential decrement of the 2D DOS. Accordingly, two approaches to compute the subthreshold slope saturation were presented, along with a description of the conductivity function employing the Kubo-Greenwood integral using with the diffusivity function in the degenerate statistics regime. Primarily, in the course of this thesis, we exposed the bell-shape mobility law that involves only the Coulomb and surface roughness mechanisms, not including the phonon scattering that does not prevail at cryogenic temperatures. Correspondingly, based on a solid electrostatic ground provided by PS simulations the origin of the intersubband scattering effect observed experimentally on the linear transfer characteristics of back-biased FDSOI structures was demonstrated. It was found that such effect is attributed to the narrow energetic separation that exists between the first two subbands, promoting thusly the subbands interaction hypothesis.

Moreover, based on Poisson-Schrödinger simulations, we confirmed our assumption of considering only the population of two subbands for FDSOI structures operating at deep cryogenic temperatures. Along with that, electrostatic parameter curves were exposed with the corresponding conduction band diagrams and an in-situ analysis of the population of different subbands for the FDSOI structure. Such approach allowed us to portrait the behavior manifested by the different electrostatic parameters, into two main event, namely the openings of the back and front channels. One appealing consequence of these two main events is the two-plateau behavior exhibited by the 𝐶 𝑔𝑐 curves in the case of positive back bias configurations. Furthermore, the performed gate-to-channel PS simulations were validated by comparison to collected C-V measurements.

Additionally, we presented an expanded version of the 1-D PS solver via the introduction of an extra dimension representing the quasi-Fermi level. Distinctly, we conducted PS simulations as well at the 𝑇 → 0𝐾 limit. Such task was made possible through the replacement of the Fermi-Dirac integral function by a Heaviside function, since such function emulates perfectly the fully degenerate metallic statistics.

Consistently, in the framework of this thesis, we presented a system of two coupled charge equations that involve a charge coupling term and a quantum shift function using the Airy solution function. The use of Airy's solution function commonly comes along with an inherent numerical pathology manifested around the null gate charge configuration. Hence, an extended form of the quantum shift function is established through the application of a technique generating a globally continuous function originating from piecewise smooth functions. Accordingly, the numerical current outputs were presented based on numerical integration formalisms that consider Fermi-Dirac statistics.

Starting from the presented system of coupled equations, the surface potentials solution is derived through a step-by-step technique ranging from initial guesses via the application of a number of error correction steps. Such approach is a good proof of the noncompulsory of the numerical integration when Fermi-Dirac statistics are applied. Moreover, closed from analytical expressions were demonstrated for diffusion and drift current computations along with a twoslope inversion charge linearization technique applicable for back-bias structures that considers the computation of the respective slopes at the source and saturation ends. Finally, a few short channel effects were implemented to the core model such as the velocity saturation, the DIBL and charge sharing, and the parasitic resistance effects. In each step, the model results were confronted to either PS simulations or experimental data in order to be validated.

Accordingly, one should emphasize that during the required years for the development of the presented model, the coherence of the model was privileged over the accuracy. Expressly, no non-homogenous expressions have been coerced in the model as customized truncation to get to describe certain operation regimes. Furthermore, the consolidation of each single step in the numerical and analytical model development process was a priority, in other word, strengthening the numerical stability of the model along with attributing a physical significance to the different involved elements was a requisite. Note as well that the number of the involved parameters (in exception with the adopted expression to describe the intersubband scattering effect) was kept limited, a selection that is commonly considered as a bonus feature for compact models. Namely, such choice of consolidation process implicates more devotion of time and efforts for such consolidation process, which could be inconspicuous in exposition works of this nature, in comparison to other more appealing processes that could have expeditious results.

Above all, the model presented here could be claimed as suitable for different geometrical, back-biasing and temperature configurations, particularly the dual channel operation mode manifested in the case of positive back biases that could be very challenging. Indeed, it should be pointed out that since the model development presented in the course of this thesis is a first attempt, it is not considered as a mature model yet. Nevertheless, for a first endeavor, it is the choice of privileging the coherence and the numerical stability of the model that should have been taken, since it allows the establishment of solid foundations, based on which upcoming efforts and improvements could be built.

In this context, the aforementioned two approaches must certainly not be treated as independent endeavors, as the quintessential approach to achieve mature compact models for cryogenic operation is to ensure a mutual feedback between the two approaches mentioned formerly. That is to say, the second approach can provide new understanding of physical phenomena at its very fundamental level, enabling the derivation of mathematical formalisms that are consistent with the physics of cryogenic operation on firsthand, and are numerically stable for these conditions on second hand. In reverse, the first approach can supply some convenient and ready-to-use formalisms to be implemented for cryogenic operation, certain short channel effects for instance.

Furthermore, and as manifested in the course of the present thesis, some conveyed cryogenic physical phenomena of quantum mechanical nature, such as the intersubband scattering, still need supplementary fundamental investigation. As the fulfillment of such fundamental studies paves the way for the attainment of closed-form analytical expressions that are independent of any device geometry or back biasing configuration.

Based upon the present thesis, additional efforts are foreseeable, namely the study of the PMOS device behavior and the implementation of additional short channel effects. Additionally, future efforts to enhance the accuracy of the model with respect to the experimental data are needed as well. Furthermore, the model presented here was built on the Python environment which implies the need to transform it into the Verilog-A language in order for it to be subject of some genuine tests in circuit simulators.

Introduction……………………………………………………………………………………. 1 Chapter 1 :

 11 Introduction ........................................................................................................................... 1.1 State-of-the-art of standard compact models at cryogenic temperatures ....................................... • The EKV MOSFET model…………………………………………………………………………………………………………7 • The MOS11/PSP models…………………………………………………………………………………………………………8 • The BSIM model…………………………………………………………………………………………………………………….9 • The L-UTSOI model………………………………………………………………………………………………………………10

Figure 2 .

 2 Figure 2. (a) Scheme of the CEA-Leti qubit [5]. (b) Colorized transmission electron microscopy image of the device along a longitudinal cross-sectional plane. Scale bar, 50 nm [7].

Figure 3 .

 3 Figure 3. Schematic representation of the two approaches.

Figure 4 .

 4 Figure 4. 100 K measured and simulated drain-to-source current versus gate-to-source and drain-to-source voltage curves of the N-MOSFET with W/L = 15.6 μm/0.16 μm [14].

Figure 5 .

 5 Figure 5. Transfer characteristics of a long (to the left) and a 28 nm short nMOS down to 4.2K, the simulations are done using the EKV-MOSFET model.

Figure 6 .

 6 Figure 6. Output characteristics of a 160 nm nMOS (left), and 40 nm nMOS (right). Measurements at 300K (dotted line), measurement at 4K (solid lines), and the MOS11/PSP model (dashed lines).

Figure 7 .

 7 Figure 7. Transfer characteristics of a 28 nm FDSOI MOSFET (EOT = 1.7 nm, and back gate EOT = 25 nm), the simulations were done using the BSIM model.

Figure 8 .

 8 Figure 8. Output characteristics of a FINFET for different gate biases, T = 8 k for the left and T = 77 K for the right.

Figure 9 .

 9 Figure 9. Gate to channel capacitance as a function of front gate bias for a 10 µm X 10 µm FDSOI nMOS (EOT1 = 1.2 nm, EOT2 = 25 nm), at 4K. L-UTSOI simulations are in solid lines and experimental data in dotted lines.

Figure 10 .

 10 Figure 10. Output characteristics for a 10 µm X 10 µm FDSOI nMOS (EOT1 = 1.2 nm, EOT2 = 25 nm), at 4K. L-UTSOI simulations are in solid lines and experimental data in dotted lines.

Figure 11 .

 11 Figure 11. Transfer characteristics (𝑉 𝑑𝑠 = 0.05 𝑉) for a 10 µm X 10 µm FDSOI nMOS (EOT1 = 1.2 nm, EOT2 = 25 nm), at 4K. L-UTSOI simulations are in solid lines and experimental data in dotted lines.

Figure 12 .

 12 Figure 12. Transfer characteristics (𝑉 𝑑𝑠 = 0.9 𝑉) for a 10 µm X 10 µm FDSOI nMOS (EOT1 = 1.2 nm, EOT2 = 25 nm), at 4K. L-UTSOI simulations are in solid lines and experimental data in dotted lines.

Figure 1 .

 1 Figure 1. The evolution of intrinsic carrier concentration 𝑛 𝑖 in silicon with temperature [9].

Figure 2 .

 2 Figure 2. Electron density as a function of the relative position of Fermi level with respect to the conduction band edge for different temperatures. Solid lines: Fermi-Dirac statistics, dashed lines: Maxwell-Boltzmann approximations. The fully degenerate limit (T = 0K) is plotted for comparison [1].

Figure 3 .

 3 Figure 3. The energetic separation of the first three subbands as a function of temperature for a back-biased FDSOI structure.

Figure 5 .

 5 Figure 5. Temperature dependence of the subthreshold slope for two different band tails.

Figure 6 .

 6 Figure 6. Effective temperature function plot for two different band tails.

Figure 7 .

 7 Figure 7.The fraction of the inversion charge populating each of the first four subbands of a 10 nm silicon film of an FDSOI structure for: (𝑎):𝑉 𝑔2 = -3𝑉, (𝑏) ∶ 𝑉 𝑔2 = 0𝑉, (𝑐) ∶ 𝑉 𝑔2 = +3𝑉

Figure 8 .

 8 Figure 8.The conduction band structure along with the energetic positions of the first four subbands for very large, very thin, and medium silicon thicknesses𝑉 𝑔1 = 1𝑉, 𝑉 𝑔2 = 3𝑉 , (𝑎):𝑡 𝑠𝑖 = 3𝑛𝑚, (𝑏) ∶ 𝑡 𝑠𝑖 = 10𝑛𝑚, (𝑐) ∶ 𝑡 𝑠𝑖 = 16𝑛𝑚, Fermi level is presented by the dashed line.

Figure 9 .

 9 Figure 9. The spatial electron distribution attributed to the four first subbands along with the total electron distribution 𝑉 𝑔1 = 1𝑉, 𝑉 𝑔2 = 3𝑉, (𝑎):𝑡 𝑠𝑖 = 3𝑛𝑚, (𝑏) ∶ 𝑡 𝑠𝑖 = 10𝑛𝑚, (𝑐) ∶ 𝑡 𝑠𝑖 = 16𝑛𝑚

Figure 10 .

 10 Figure 10. The conductivity function with exponential band tail (𝜇 0 = 0.1 𝑚 2 𝑉𝑠 ⁄ , ∆𝐸 = 3𝑚𝑒𝑉)

Figure 11 .

 11 Figure 11. The subthreshold slope dependence with temperature as given by the carrier density approach and the conductivity approach.

Figure 12 .

 12 Figure 12. the diffusivity function dependence on temperature in both the linear scale (a) and log scale (b) of the temperature axes, the simulation is made for three different 2D electrons densities (as shown in the legend).

Figure 13 .

 13 Figure 13. The diffusivity function dependence on 2D electrons density (a) and on the Fermi level (b).

Figure 14 .

 14 Figure 14. a) Experimental (symbols) variation of effective mobility 𝜇 𝑒𝑓𝑓 with 2D carrier density n at T = 4.2 K and T = 200 K. Typical 𝜇 𝑒𝑓𝑓 experimental data from 10𝜇𝑚 × 10𝜇𝑚 28 nm FDSOI MOSFET with 1.8 nm EOT gate oxide, 7 nm silicon and 30 nm bottom oxide. b) Corresponding experimental (symbols) variation of diffusivity D with 2D carrier density n at T = 4.2 K and T = 200 K. The diffusivity D is calculated using Eq 2.24. The dashed lines show the linear trends obtained for a constant mobility (Resp. μ= 0.108 and 0.05 m²/Vs) using the metallic statistics limit of Eq 2.22.

Figure 15 .

 15 Figure 15. The transfer characteristics of a back-biased FDSOI transistor measured at temperature T = 4.2 K, the demonstrated oscillations corresponds to the Intersubband scattering phenomena.

Figure 16 .

 16 Figure 16. The corresponding energetic separation between the first subband plotted along with the output characteristic in the case of a back bias of 𝑉 𝑔2 = 3𝑉.

Figure 1 .

 1 Figure 1. Typical scheme of an FDSOI structure.

Figure 2 .

 2 Figure 2. Gate-to-channel capacitance curves as a function of front gate bias and for different back bias (t ox =1nm, t box =25nm, t si =10nm, T=4K).

Figure 3 .

 3 Figure 3. Captured moment of the onset of the first, second and third subband, as well as the creation of the front potential well for: Vg2=3V t ox =1nm, t box =25nm, t si = 7 nm, T = 4K.

Figure 4 .

 4 Figure 4. Inversion charge as a function of front gate voltage and for three different back biases, in both the linear and the log scale.

Figure 5 .

 5 Figure 5. Front surface potential as a function of front gate voltage and for three different back biases.Likewise, Figure6demonstrates the back surface potential as a function of front gate voltage and for different back biases. Despite the expected two limbs behavior, where the first limb characterizes the weak inversion region before the first subband onset, and the second limb characterize the strong inversion region after the first subband onset, note the earlier saturation of the curve corresponding to positive back bias, ensued by the one corresponding to the null than then negative back bias.

Figure 6 . 1 . 2 . 2

 6122 Figure 6. Back surface potential as a function of front gate voltage and for three different back biases

Figure 7 .

 7 Figure 7. Typical band diagram and total electron distribution from PS simulation for a FDSOI structure (Vg1=1V, tox1=1nm, tox2= 25nm, tsi= 7nm, Vg2 = -3, 0, +3V, T=4K).

Figure 8 .

 8 Figure 8. Electron distrubution in the channel for a thin silicon film: t si =7nm presented by the red line, and a thick silicon film: t si = 16 nm presented by the blue dashed line (t ox =1nm, t box =25nm, T=4K, V g1 =0.73V, V g2 =+3V).

Figure 9

 9 Figure9shows the variations of the inversion charge 𝑄 𝑖𝑛𝑣 in the silicon channel as a function of front gate voltage 𝑉 𝑔1 with a back bias 𝑉 𝑔2 = +3𝑉, obtained from PS simulations for various temperatures T = 0 -60 K.

Figure 9 .

 9 Figure 9. Inversion charge as a function of front gate voltage for different temperatures (t ox1 =1nm, t ox2 = 25nm, t si = 7nm, V g2 = +3V).

Figure 10 .

 10 Figure 10. Charge centroid position in the silicon film as a function of temperature (t ox1 =1nm, t ox2 =25nm, t si =10nm, V g2 =+3V, V g1 =0.3V) Moreover, Figure 11 demonstrates 𝑄 𝑖𝑛𝑣 (𝑉 𝑔1 ) curves for different temperatures and for two different silicon thicknesses at positive back bias. Indeed, the curves corresponding to the silicon thickness 𝑡 𝑠𝑖 = 10 𝑛𝑚 show an elbow at 𝑉 𝑔1 = 0.6 𝑉 which starts becoming visible for temperatures lower than 20K.

Figure 11 .

 11 Figure 11. Inversion charge curves as a function of front gate voltage, computed for different temperatures, and for two different silicon thicknesses.

Figure 12 .

 12 Figure 12. Demonstration of multi-subband population at T = 300 K.

Figure 13 .

 13 Figure 13. Gate to channel curves as a function of front gate voltage, and for different temperatures, (t ox1 =1nm, t ox2 =25nm, t si =10nm, V g2 =+3V).

Figure 14 .

 14 Figure 14. Electron gas profile for different temperatures, and for a positive back biasIn a similar manner, the front and back surface potentials as functions of the front gate voltage are extracted from the simulations results and presented in Figure15.

Figure 15 .

 15 Figure 15. Front surface potential as a function of front gate voltage, for three different temperatures and for a null and a positive back bias.

Figure 16 .

 16 Figure 16. Back surface potential as a function of front gate voltage, for three different temperatures, and for a null and a positive back bias.

Figure 17 .

 17 Figure 17. Inversion charge as a function of front gate voltage and for four different silicon thicknesses.

Figure 20 .

 20 Figure 20. Front surface potential as a function of front gate voltage and for four different silicon thicknesses.

Figure 21 .

 21 Figure 21. Back surface potential as a function of front gate voltage and for four different silicon thicknesses.

Figure 22 .

 22 Figure 22. Experimental validation of the simulated gate to channel curves as a function of front gate voltage and for different back biases and for a temperature of 4.2 K.

Figure 22

 22 Figure 22 demonstrates the good agreement between PS simulations (solid lines) and the measured data (symbols) for several back biases 𝑉 𝑔2 = -3, -2 … 4 𝑉. Indeed, the back bias influence can be well described by the simulations down to the temperature T = 4.2 K and the two-plateau behavior predicted earlier by the PS simulations is well embodied in the measurements as well. Substantially, the threshold voltage dependence on back biasing along with the change in the 𝐶 𝑔𝑐 (𝑉 𝑔1 ) curves slope are well described by PS simulations.

Figure 23 .

 23 Figure 23. Experimental validation of the simulated gate to channel curves as a function of front gate voltage and for different back biases and for a temperature of 20K.

Figure 24 .

 24 Figure 24. Experimental validation of the simulated gate to channel curves as a function of front gate voltage and fordifferent back biases and for a temperature of 50K.

Figure 25 .

 25 Figure 25. Experimental validation of the simulated gate to channel curves as a function of front gate voltage and for different back biases and for a temperature of 100K.

Figure 23 illustrates

 23 Figure 23 illustrates the computed transfer characteristics 𝐼 𝑑 = 𝑓(𝑉 𝑔1 ) in the linear regime, i.e. for 𝑉 𝑑 = 0.05 𝑉, and the saturated regime, i.e. for 𝑉 𝑑 = 1𝑉, as well as the output characteristics 𝐼 𝑑 = 𝑓(𝑉 𝑑 ), for the indicated FDSOI transistor at the temperature 𝑇 = 4 𝐾.

Figure 23 .

 23 Figure 23. Computed Transfer and output characteristics of FDSOI transistor at T = 4K.

Figure 2 .

 2 Figure 2. Front and back surface potentials, total inversion charge density, and gate-to-channel curves as functions of frontgate bias; the singularity is indicated by the red dotted ellipses.

Figure 4 .

 4 Figure 4. Front and back surface potentials, total inversion charge density, and gate-to-channel capacitance curves asfunctions of front gate bias computed using the Equations: Eq 4.10, Eq 4.7, Eq 4.6, and Eq 4.3.

{

  𝑄 𝑖𝑛𝑣1 = (𝑄 𝑔1 + 𝑄 𝑔2 ).𝑚𝑎𝑥(𝑄 𝑔1 , 0) 𝑚𝑎𝑥(𝑄 𝑔1 , 0) + 𝑚𝑎𝑥(𝑄 𝑔2 , 0)𝑄 𝑖𝑛𝑣2 = (𝑄 𝑔1 + 𝑄 𝑔2 ). 𝑚𝑎𝑥(𝑄 𝑔2 , 0)𝑚𝑎𝑥(𝑄 𝑔1 , 0) + 𝑚𝑎𝑥(𝑄 𝑔2 , 0)Eq 4.12Finally, Figure5demonstrates the comprehensive trend of the front quantum shift function as a function of the front gate voltage ∆𝑉(𝑄 𝑔1 ) in the case of a positive back bias configuration.

Figure 5 .

 5 Figure 5. Front Quantum shift expression as a function of front gate charge in the case of a positive back bias, derived numerically and from PS data.

Figure 6 .

 6 Figure 6. Front Quantum shift expression as a function of front gate charge in the case of a positive back bias, derived numerically and from PS data, for two different silicon thicknesses.

Figure 7 .

 7 Figure 7. Front and back inversion charges curves as a function of front gate bias evaluated using both of the aforementioned approaches.

Figure 8 .

 8 Figure 8. Front surface potential curves as functions of front gate voltage for three different back biase, two different silicon thicknesses, and for a temperature of 4 K, the solid lines represent the numerical model results and symbols represent the PS simulation results.

Figure 9 .

 9 Figure 9. Back surface potential curves as functions of front gate voltage for three different back biases, two different silicon thicknesses, and for a temperature of 4 K, the solid lines represent the numerical model results and symbols represent the PS simulation results.

Figure 10 .

 10 Figure 10. Front inversion charge density curves as functions of front gate voltage for three different back biases, two different silicon thicknesses, and for a temperature of 4 K, the solid lines represent the numerical model results and symbols represent the PS simulation results.

Figure 11 .

 11 Figure 11. Back inversion charge density curves as functions of front gate voltage for three different back biases, two different silicon thicknesses, and for a temperature of 4 K, the solid lines represent the numerical model results and symbols represent the PS simulation results.

Figure 12 .

 12 Figure 12. Gate-to-channel capacitance curves as functions of front gate voltage for three different back biases and two different silicon thicknesses, the solid lines represent the numerical model results and symbols represent the PS simulation results.

Figure 13 .

 13 Figure 13. Transfer and output characteristics computed for three different back biases using the numerical model with the incorporation of the bell-shaped mobility law.

Figure 15 .

 15 Figure 15. Transfer and output characteristics curves computed using a constant mobility and for three different back biases and two different silicon thicknesses, the solid lines represent the numerical model results and symbols represent the PS simulation results.

Figure 1 .

 1 Figure 1. the kinetics of the convergence of the front/back surface potentials throughout the four steps of error correction.

Figure 2 .

 2 Figure 2. Front and back surface potential curves as functions of front gate voltage for three different back biases, the solid red lines represent the numerical model results, and the blue dashed ones represent the analytical model results.

Figure 3 .

 3 Figure 3. Front and back inversion density curves as functions of front gate voltage for three different back biases, the solid red lines represent the numerical model results, and the blue dashed ones represent the analytical model results.

Figure 4 .

 4 Figure 4. Gate-to-channel capacitance curves as functions of front gate voltage for three different back biases, the solid red lines represent the numerical model results, and the blue dashed ones represent the analytical model results.

Figure 5 .

 5 Figure 5. scheme Illustrating the pinch-off point.

Figure 6

 6 Figure 6 compares the numerical and the analytical approximated solutions of the diffusion integral with 𝑟 = 1.8, the constant 𝑟 is chosen meticulously with the aim of minimizing the relative error between the analytical approximated solution of the integral and the numerical exact one.

Figure 6 .

 6 Figure 6. Comparison between numerical (red lines) and analytical (blue triangles) calculations of the diffusion term integral.

Figure 7 .

 7 Figure 7. The dependence of inversion charge densities -𝑄 𝑖𝑛𝑣1,2 with the corresponding subband potential 𝜒 𝑠1,2 for three different back biases.

Figure 9 .

 9 Figure 9. The capacitive scheme of a 2-D electrostatic analysis for an FDSOI structure.

Figure 10 .

 10 Figure 10. simple scheme representing two access parasitic resistances and the corresponding altered biases.

Figure 11 .

 11 Figure 11. Parameters 𝑄 𝑐1 and 𝑄 𝑐2 as functions of the back bias.

  Accordingly, the following Figures illustrate the analytical

  model results in comparison to the experimental data for long channel transistor of channel length 𝐿 = 1𝜇𝑚, for different back biases, with respect to the linear and saturated transfer characteristics, and the output characteristics.

Figure 12 .

 12 Figure 12. The linear transfer characteristics of long channel transistor (𝐿 = 1𝜇𝑚 and 𝑊 = 1𝜇𝑚) for different back biases: experimental data (to the left) compared to the analytical model results (to the right).

Figure 13 .

 13 Figure 13. The saturated transfer characteristics of long channel transistor (𝐿 = 1𝜇𝑚 and 𝑊 = 1𝜇𝑚) for different back biases: experimental data (to the left) compared to the analytical model results (to the right).

Figure 14 .

 14 Figure 14. The output characteristics of long channel transistor (𝐿 = 1𝜇𝑚 and 𝑊 = 1𝜇𝑚) for different back biases: experimental data (to the left) compared to the analytical model results (to the right).On the other hand, such effect is not present in the short channel transistors i.e. 𝐿 = 120 𝑛𝑚 and 𝐿 = 30 𝑛𝑚 in the present study. For instance, regarding the 𝐿 = 120 𝑛𝑚 short channel transistor, the analytical model results produced in Figure15, Figure16, and Figure17using the following set of parameters: 𝑉 𝑓𝑏1 𝑉 = -0.09 𝑉, 𝑉 𝑓𝑏2 = 0.65 𝑉, 𝜇 𝑚𝑎𝑥1 = 0.04 𝑚 2 𝑉. 𝑠 ⁄ , 𝜇 𝑚𝑎𝑥2 = 0.25 𝑚 2 𝑉. 𝑠 ⁄ ,𝑄 𝑟𝑒𝑓 = 10 3 𝐶/𝑚 2 ; 𝑄 𝑐1 = 20 * 10 -4 𝐶/𝑚 2 , 𝑄 𝑐2 = 20 * 10 -4 𝐶/𝑚 2 .

Figure 15 .

 15 Figure 15. The linear transfer characteristics of short channel transistor (𝐿 = 120 𝑛𝑚 and 𝑊 = 1𝜇𝑚) for different back biases: experimental data (to the left) compared to the analytical model results (to the right).

Figure 16 .

 16 Figure 16. The saturated transfer characteristics of short channel transistor (𝐿 = 120 𝑛𝑚 and 𝑊 = 1𝜇𝑚) for different back biases: experimental data (to the left) compared to the analytical model results (to the right).

Figure 17 .

 17 Figure 17. The output characteristics of short channel transistor (𝐿 = 120 𝑛𝑚 and 𝑊 = 1𝜇𝑚) for different back biases: experimental data (to the left) compared to the analytical model results (to the right). Moreover, using the following values of parameters: 𝑉 𝑓𝑏1 = -0.09 𝑉, 𝑉 𝑓𝑏2 = 0.65 𝑉, 𝜇 𝑚𝑎𝑥1 = 0.08 𝑚 2 𝑉. 𝑠 ⁄ , 𝜇 𝑚𝑎𝑥2 = 0.46 𝑚 2 𝑉. 𝑠 ⁄ ,𝑄 𝑟𝑒𝑓 = 10 3 𝐶/𝑚 2 ; 𝑄 𝑐1 = 5 * 10 -4 𝐶/𝑚 2 and 𝑄 𝑐2 = 10 -3 𝐶/𝑚 2 , the next plots are produced for short channel transistor of 𝐿 = 30 𝑛𝑚 channel length, as depicted in Figure 18, Figure 19, and Figure 20.

Figure 18 .

 18 Figure 18. The linear transfer characteristics of short channel transistor (𝐿 = 30 𝑛𝑚 and 𝑊 = 1𝜇𝑚) for different back biases: experimental data (to the left) compared to the analytical model results (to the right).

Figure 19 .

 19 Figure 19. The saturated transfer characteristics of short channel transistor (𝐿 = 120 𝑛𝑚 and 𝑊 = 1𝜇𝑚) for different back biases: experimental data (to the left) compared to the analytical model results (to the right).

Figure 20 .

 20 Figure 20. The output characteristics of short channel transistor (𝐿 = 30 𝑛𝑚 and 𝑊 = 1𝜇𝑚) for different back biases: experimental data (to the left) compared to the analytical model results (to the right).

Table 1 .

 1 MOS11 and PSP Modified parameters for 160 nm and 40 nm NMOS[21].

  2/3 for the 2/3 power region, are added to ensure that the result is null when the gate charge density is null, following Eq 4.18. Note that, unreservedly, such expression of the potential quantum shift function is suitable for all geometrical configurations.

	𝛥𝑉(𝑄 𝑔1,2 , 𝛿𝑄 𝑔 )
	= + 𝛽 𝑄𝑀 . (𝑄 𝑝𝑜𝑠 (𝑄 𝑔1,2 ) 1 𝛾 . (𝑄 𝑛𝑒𝑔 (𝑄 𝑔1,2 ) + 2/3 -( 𝛿𝑄 𝑔 ) 2	2 𝛿𝑄 𝑔	)	2/3	)	Eq 4.18

  1,2 and 𝑄 𝑖𝑛𝑣1,2 in the source and drain ends respectively, resulting in Eq 4.22.

	𝐼 𝑑𝑟𝑖𝑓𝑡1,2 = ∫	𝜒 s1,2 𝑑𝑟𝑛	𝜇 𝑛 . 𝑄 𝑖𝑛𝑣 (𝑉 𝑠 )	. 𝑑𝜒 𝑠
		𝜒 s1,2 𝑠𝑟𝑐		Eq 4.22
	{ 𝐼 𝑑𝑖𝑓𝑓1,2 = ∫ 𝑄 inv1,2 𝑑𝑟𝑛 𝑄 inv1,2 𝑠𝑟𝑐	𝑄 𝑖𝑛𝑣 . ( 𝑞 𝑁 2𝐷 𝑘𝑇	1 1 -𝑒 -𝑄 𝑖𝑛𝑣1,2 𝑞 𝑁 2𝐷	𝑑 𝑄 𝑖𝑛𝑣 )

  [𝜙 𝑖𝑚 𝑗 , 𝜙 𝑖𝑚 𝑗+1 ] the integrated function 𝑄 𝑖𝑛𝑣1,2 (𝜒 𝑠1,2 ) is interpolated by a polynomial of degree-one passing through the points (𝜒 𝑠1,2 𝑗 , 𝑄 𝑖𝑛𝑣1,2 (𝜒 𝑠1,2 𝑗 )) and (𝜒 𝑠1,2 𝑗+1 , 𝑄 𝑖𝑛𝑣1,2 (𝜒 𝑠1,2 𝑗+1 )), allowing to approximate the drift integral using the formula in Eq 4.23. The term 𝑄 𝑖𝑛𝑣1,2 𝑗,𝑚𝑒𝑑 designate the median value of the corresponding inversion charge in the segment [𝑗, 𝑗 + 1].

  are plainly reformed through the implantation of the corresponding function, as in Eq 4.24.

		𝐼 𝑑𝑟𝑖𝑓𝑡1,2 ≅ -	𝑊 𝐿	. 𝜇 𝑚𝑎𝑥 1,2 . ∑ 𝑖 𝑚𝑎𝑥 𝑖=0	𝑄 𝑖𝑛𝑣1,2 𝑗,𝑚𝑒𝑑 ⁄ 𝑄 𝑐 1 + ( 𝑄 𝑐 ⁄ ) 𝑄 𝑖𝑛𝑣1,2 𝑗,𝑚𝑒𝑑 2 . (𝜒 𝑠1,2 𝑗+1 -𝜒 𝑠1,2 𝑗 )	. 𝑄 𝑖𝑛𝑣1,2 𝑗,𝑚𝑒𝑑
	{	𝐼 𝑑𝑖𝑓𝑓1,2 =	𝑊 𝐿	. 𝑘𝑇𝑞 𝑁 2𝐷 . 𝜇 𝑚𝑎𝑥1,2 . ∫	𝑄 𝑖𝑛𝑣1,2 𝑗,𝑚𝑒𝑑 ⁄ 𝑄 𝑐 1 + ( 𝑄 𝑖𝑛𝑣1,2 𝑗,𝑚𝑒𝑑 𝑄 𝑐 ⁄ )

  [START_REF] Ghibaudo | On the diffusion current in a MOSFET operated down to deep cryogenic temperatures[END_REF], for which the derivation of the closed form expression of the slope of the quasi-Fermi level gives birth to the diffusion and drift integral terms. The derivation of the analytical expression for integrals is detailed in the next two sections separately.

	𝐼 𝑑1,2 = -	𝑊 𝐿	. ∫ 𝜇 𝑒𝑓𝑓 . 𝑄 𝑖𝑛𝑣1,2 (𝑦). 0 𝐿	𝑑𝑦 𝑑𝜙 𝑖𝑚	(𝑦)	. 𝑑𝑦	Eq 5.29

1.2.2 The drift current component computation:

  At odds with the diffusion term integral, the drift term integral has two different integrals for the effective mobility and the mobility function cases. The general expression of drift current integral is recalled in Eq 5.34.

	𝐼 𝑑𝑟𝑖𝑓𝑡1,2 = ∫ 𝜒 s1,2 𝑑𝑟𝑛	𝑄 𝑖𝑛𝑣 (𝑉 𝑠 )	. 𝑑𝜒 𝑠	Eq 5.34
	𝜒 s1,2 𝑠𝑟𝑐

  1,2 𝑚 . (𝑉 𝑠1,2 -𝑉 1,2 𝑚 )

						Eq 5.35
	Where:
	𝛼 1,2 𝑚 =	𝜕𝑄 𝑖𝑛𝑣1,2 𝜕𝑉 𝑠1,2	|	𝑉 𝑠1,2 =𝑉 s1,2 𝑚	Eq 5.36
	{	𝑉 s1,2 𝑚 =	𝑉 1,2 𝑠𝑟𝑐 + 𝑉 1,2 𝑑𝑟𝑛 2

  𝑖𝑛𝑣 𝑙𝑖𝑛 1,2 = -𝑄 1,2 𝑠𝑟𝑐 + 𝜆 1,2 𝑠𝑟𝑐 (𝜒 𝑠1,2 -𝜒 s1,2 𝑠𝑟𝑐 ) -𝑄 𝑖𝑛𝑣 𝑙𝑖𝑛 1,2 ̃= -𝑄 1,2 𝑠𝑎𝑡 + 𝜆 1,2 𝑠𝑎𝑡 (𝜒 𝑠1,2 -𝜒 s1,2 𝑠𝑎𝑡 )Where the slopes 𝜆 1,2 𝑠𝑟𝑐 and 𝜆 1,2 𝑠𝑎𝑡 are defined as the derivatives of the inversion charge densities with respect to the subband potential at the source and saturation side respectively: Doing so, we now have a two-part piecewise linear inversion charge functions. As a consequence, our drift integral throughout the channel will be performed in two parts, where the first part is computed form the source point the intersection point, and the second part is computed from the intersection point to the saturation point, as described in Eq 5.39:

						Eq 5.37
		𝜆 1,2 𝑠𝑟𝑐 =	𝜕𝑄 𝑖𝑛𝑣1,2 𝜕𝜒 𝑠1,2	|	𝜒 𝑠1,2 =𝜒 s1,2 𝑠𝑟𝑐	Eq 5.38
		{ 𝜆 1,2 𝑠𝑎𝑡 =	𝜕𝑄 𝑖𝑛𝑣1,2 𝜕𝜒 𝑠1,2	|	𝜒 𝑠1,2 =𝜒 s1,2 𝑠𝑎𝑡
	𝜒 s1,2 𝑠𝑎𝑡
	∫ 𝑄 𝑖𝑛𝑣 (𝜒 𝑠 ). 𝑑𝜒 𝑠	= 𝐼 𝑑𝑟𝑖𝑓𝑡 1,2 𝑝𝑎𝑟𝑡1	+ 𝐼 𝑑𝑟𝑖𝑓𝑡 1,2 𝑝𝑎𝑟𝑡2
	𝜒 s1,2 𝑠𝑟𝑐	𝜒 s1,2 𝑖𝑛𝑡				Eq 5.39
	𝐼 𝑑𝑟𝑖𝑓𝑡 1,2 𝑝𝑎𝑟𝑡1	= ∫ (𝑄 1,2 𝑠𝑟𝑐, + 𝜆 1,2 𝑠𝑟𝑐, (𝜒 𝑠1,2 -𝜒 s1,2 𝑠𝑟𝑐 )) . 𝑑𝜒 𝑠
		𝜒 s1,2 𝑠𝑟𝑐 𝜒 s1,2 𝑠𝑎𝑡
	𝐼 𝑑𝑟𝑖𝑓𝑡 1,2 𝑝𝑎𝑟𝑡2	= ∫ (𝑄 1,2 𝑠𝑎𝑡 + 𝜆 1,2 𝑠𝑎𝑡 (𝜒 𝑠1,2 -𝜒 s1,2 𝑠𝑎𝑡 )) . 𝑑𝜒 𝑠
	{	𝜒 s1,2 𝑖𝑛𝑡

  1,2 𝑠𝑎𝑡 -𝑄 1,2 𝑠𝑟𝑐 ) + (𝜆 1,2 𝑠𝑟𝑐 . 𝜒 s1,2 𝑠𝑟𝑐 -𝜆 1,2 𝑠𝑎𝑡 . 𝜒 s1,2 𝑠𝑎𝑡 ) 𝜆 1,2 𝑠𝑟𝑐 -𝜆 1,2 𝑠𝑎𝑡

	Eq 5.40

  -𝑉 0 -𝜙 𝑖𝑚 -∆𝑉(𝑄 𝑔1,2 ) 𝑘𝑇 Since by definition the gate charges are functions of 𝑉 𝑠 1,2 and not of 𝜒 𝑠 1,2 , we need to compute the derivatives with respect to 𝑉 𝑠 1,2 in first place, then transform the output following the equivalence depicted in Eq 5.42:

												Eq 5.41
	𝑑𝑄 𝑖𝑛𝑣1,2 𝑑𝜒 𝑠1,2	=	𝑑𝑄 𝑖𝑛𝑣1,2 𝑑𝑉 𝑠1,2	.	1 𝑑𝜒 𝑠1,2 𝑑𝑉 𝑠1,2 ⁄	=	𝑑𝑄 𝑖𝑛𝑣1,2 𝑑𝑉 𝑠1,2	.	1 -	1 𝑑∆𝑉(𝑄 𝑔1,2 ) ⁄ 𝑑𝑉 𝑠1,2	Eq 5.42
	The 𝑑𝑄 𝑖𝑛𝑣1,2 𝑑𝑉 𝑠1,2 ⁄	term can be formulated in the next manner:
	𝑑𝑄 𝑖𝑛𝑣1,2 𝑑𝑉 𝑠1,2	=	𝑑𝑄 𝑖𝑛𝑣1,2 𝑑𝜙 𝑖𝑚	.	𝑑𝜙 𝑖𝑚 𝑑𝑉 𝑠1,2		Eq 5.43
					= 𝑞 𝐴 2𝐷 (	𝑑𝑉 𝑠1,2 𝑑𝜙 𝑖𝑚	-1 -	𝑑∆𝑉(𝑄 𝑔1,2 ) 𝑑𝜙 𝑖𝑚	)	𝑒 𝐴𝑟𝑔 1,2 1 + 𝑒 𝐴𝑟𝑔 1,2

  1,2 We derive both equations with respect to the quasi-Fermi level 𝜙 𝑖𝑚 :

		𝐶 𝑜𝑥1 . (𝑉 𝑔1 -𝑉 𝑓𝑏1 -𝑉 𝑠1 ) = 𝑞𝑘𝑇𝐴 2𝐷 𝑙𝑛 (1 + 𝑒𝑥𝑝 (	𝑉 𝑆1 -𝑉 0 -𝜙 𝑖𝑚 𝑘𝑇	)) + 𝐶 𝑠𝑖 . (𝑉 𝑠1 -𝑉 𝑠2 )
	{	𝐶 𝑜𝑥2 . (𝑉 𝑔2 -𝑉 𝑓𝑏2 -𝑉 𝑠2 ) = 𝑞𝑘𝑇𝐴 2𝐷 𝑙𝑛 (1 + 𝑒𝑥𝑝 (	𝑘𝑇 𝑉 𝑆2 -𝑉 0 -𝜙 𝑖𝑚	)) + 𝐶 𝑠𝑖 . (𝑉 𝑠2 -𝑉 𝑠1 )	Eq 5.45
			-𝐶 𝑜𝑥1	𝑑𝑉 𝑠1 𝑑𝜙 𝑖𝑚	= 𝑞 𝐴 2𝐷 ( 𝑑𝜙 𝑖𝑚 𝑑𝑉 𝑠1	-1 -	𝑑∆𝑉(𝑄 𝑔1 ) 𝑑𝜙 𝑖𝑚 )	𝑒 𝐴𝑟𝑔 1 1 + 𝑒 𝐴𝑟𝑔 1	+ 𝐶 𝑠𝑖 . (	𝑑𝑉 𝑠1 𝑑𝜙 𝑖𝑚	-	𝑑𝑉 𝑠2 𝑑𝜙 𝑖𝑚	)	Eq 5.46
	{ -𝐶 𝑜𝑥2	𝑑𝑉 𝑠2 𝑑𝜙 𝑖𝑚	= 𝑞 𝐴 2𝐷 ( 𝑑𝜙 𝑖𝑚 𝑑𝑉 𝑠2	-1 -	𝑑∆𝑉(𝑄 𝑔2 ) 𝑑𝜙 𝑖𝑚 )	𝑒 𝐴𝑟𝑔 2 1 + 𝑒 𝐴𝑟𝑔 2	+ 𝐶 𝑠𝑖 . (	𝑑𝑉 𝑠2 𝑑𝜙 𝑖𝑚	-	𝑑𝑉 𝑠1 𝑑𝜙 𝑖𝑚	)
	Then, we replace	𝑑∆𝑉(𝑄 𝑔1,2 ) 𝑑𝜙 𝑖𝑚	by	𝑑∆𝑉(𝑄 𝑔1,2 ) 𝑑𝑉 𝑠1,2	𝑑𝑉 𝑠1,2 𝑑𝜙 𝑖𝑚	and rearrange the equations:
			[𝐶 𝑜𝑥1 + 𝐶 𝑠𝑖 + 𝑞 𝐴 2𝐷 (1 -	𝑑∆𝑉(𝑄 𝑔1 ) 𝑑𝑉 𝑠1	)	𝑒 𝐴𝑟𝑔 1 1 + 𝑒 𝐴𝑟𝑔 1	]	𝑑𝑉 𝑠1 𝑑𝜙 𝑖𝑚	= 𝑞 𝐴 2𝐷	𝑒 𝐴𝑟𝑔 1 1 + 𝑒 𝐴𝑟𝑔 1	+ 𝐶 𝑠𝑖 .	𝑑𝑉 𝑠2 𝑑𝜙 𝑖𝑚	Eq 5.47
	{	[𝐶 𝑜𝑥2 + 𝐶 𝑠𝑖 + 𝑞 𝐴 2𝐷 (1 -	𝑑∆𝑉(𝑄 𝑔2 ) 𝑑𝑉 𝑠2	)	𝑒 𝐴𝑟𝑔 2 1 + 𝑒 𝐴𝑟𝑔 2	]	𝑑𝑉 𝑠2 𝑑𝜙 𝑖𝑚	= 𝑞 𝐴 2𝐷	𝑒 𝐴𝑟𝑔 2 1 + 𝑒 𝐴𝑟𝑔 2	+ 𝐶 𝑠𝑖 .	𝑑𝑉 𝑠1 𝑑𝜙 𝑖𝑚
													Eq 5.44
	By analyzing the expression of Eq 5.44 we can see that the computation of the two intermediary
	d∆V(Q g1,2 ) dV s1,2 ⁄	and 𝑑𝑉 𝑠1,2 𝑑𝜙 𝑖𝑚 ⁄	terms is needed in order to get to the final expression of
	𝑑𝑄 𝑖𝑛𝑣1,2 𝑑𝜒 𝑠1,2 ⁄	, since we already have the closed-form expression of d∆V(Q g1,2 ) dV s1,2 ⁄	from
	the procedure of surface potential error expression procedure, we proceed into the derivation of
	an analytical expression of 𝑑𝑉 𝑠1,2 𝑑𝜙 𝑖𝑚 ⁄	in the next paragraph.
	For the computation of 𝑑𝑉 𝑠1,2 𝑑𝜙 𝑖𝑚 ⁄	we start from our starting set of equations:

  Correspondingly, the 𝑑𝑉 𝑠2 𝑑𝜙 𝑖𝑚 ⁄ term can be obtained by a simple summation of equation (1) from the previous system multiplied by 𝐶 𝑠𝑖 , and equation (2) from the previous system multiplied by 𝐶 1 , giving:

				⇒	𝑑𝑉 𝑠1 𝑑𝜙 𝑖𝑚	=	𝑠 1 𝐶 2 + C si 𝑠 2 𝐶 1 𝐶 2 -C si 2
	𝐶 1 𝐶 2	𝑑𝑉 𝑠2 𝑑𝜙 𝑖𝑚	= C si 𝑠 1 + C 1 𝑠 2 + C si	2 .	𝑑𝑉 𝑠2 𝑑𝜙 𝑖𝑚	Eq 5.52
				⇒	𝑑𝑉 𝑠2 𝑑𝜙 𝑖𝑚	=	𝑠 2 𝐶 1 + C si 𝑠 1 𝐶 1 𝐶 2 -C si 2
							Eq 5.50
			𝐶 2	𝑑𝑉 𝑠2 𝑑𝜙 𝑖𝑚	= 𝑠 2 + C si .	𝑑𝑉 𝑠1 𝑑𝜙 𝑖𝑚	(2)
	At this point the 𝑑𝑉 𝑠1 𝑑𝜙 𝑖𝑚 ⁄	term can be obtained by a simple summation of equation (1) from
	the previous system multiplied by 𝐶 2 , and equation (2) from the previous system multiplied by
	𝐶 𝑠𝑖 , yielding:
	𝐶 1 𝐶 2	𝑑𝑉 𝑠1 𝑑𝜙 𝑖𝑚	= 𝑠 1 𝐶 2 + C si 𝑠 2 + C si	2 .	𝑑𝜙 𝑖𝑚 𝑑𝑉 𝑠1	Eq 5.51

The name "NP-complete" is short for "nondeterministic polynomial-time complete". In this name, "nondeterministic" refers to nondeterministic Turing machines, a way of mathematically formalizing the idea of a brute-force search algorithm. Polynomial time refers to an amount of time that is considered "quick" for a deterministic algorithm to check a single solution, or for a nondeterministic Turing machine to perform the whole search. "Complete" refers to the property of being able to simulate everything in the same complexity class.

. 𝑢 1,2 1 -𝑒 -𝑢 1,2 𝑑𝑢 𝑢 1,2𝑑𝑟𝑛 𝑢 1,2𝑠𝑟𝑐Eq 4.24

Acknowledgments

* (𝐶 𝑜𝑥1 * + 𝐶 𝑠𝑖 )𝑉 𝑔2 * + 𝐶 𝑠𝑖 𝐶 𝑜𝑥1 * 𝑉 𝑔1 * = ((𝐶 𝑜𝑥2 * + 𝐶 𝑠𝑖 )(𝐶 𝑜𝑥1 * + 𝐶 𝑠𝑖 ) -𝐶 𝑠𝑖 2 ) 𝑉 𝑠2 Eq 5.12

Note that although the creation of the second potential well manifests on the same polarization for all silicon thicknesses, the onset of the second subband is accordingly a function of silicon thickness, signifying that the onset of the second subband will happen in a rather higher polarization for thinner silicon films, as shown in Figure 18. Likewise, in Figure 19 is reported the impact of silicon channel thickness on the 𝐶 𝑔𝑐 (𝑉 𝑔1 ) characteristics at T=4K and for a positive back bias 𝑉 𝑔2 = +3𝑉, demonstrating a delayed opening of the gate-to-channel curves for reduced silicon thickness, due to the delayed onset of the first subband.

Note also the significant rounding exhibited before front channel opening for 𝑡 𝑠𝑖 = 7𝑛𝑚, which can be explained by the higher carrier profile overlap between the back and front interface electron distributions as 𝑡 𝑠𝑖 is reduced. In contrast, for larger 𝑡 𝑠𝑖 values, the overlap is decreased, so that front and back channels are better separated.

Chapter 4:

The numerical model literature, for instance, the absolute value function. Undeniably, the resulting function is continuous and differentiable everywhere, however, for the 𝑄 𝑔1,2 = 0 configuration, the corresponding electrostatic quantities and their respective derivative functions present a numerical pathology.

The manifested numerical pathology:

Curiously, it was found that the modelling of quantum mechanical effects with the mere assumption of an infinite triangular potential well at the film-oxide interface, i.e. the general 2/3 power of the quantum shift function, can result in an unappealing behavior of the surface potentials, which will in turn propagates to the inversion charges, and becomes significant in the behavior of the corresponding first-order derivative i.e. the gate-to-channel capacitance. Such singularity is manifested around the point where the slope of the potential profile function 𝑉(𝑥) is reversed, i.e. the zero-gate-charge point. Figure 2 depicts the displayed numerical pathology in various electrostatic quantities. It should be noted that such numerical pathology is manifest primarily in the positive back bias configuration.

In Figure 2 the different electrostatic quantities are presented as functions of front gate voltage. Seemingly, the singularity (designated by the red dotted ellipses) is manifested in the front surface potential, total inversion charge density, and the gate-to-channel curves for a back bias of 𝑉 𝑔2 = +3𝑉.

Numerical current model validation by comparison to PS results:

Equivalently to the section 1.1.5, in this final section and through Figure 15 we expose the computed transfer characteristics 𝐼 𝑑 = 𝑓(𝑉 𝑔1 ) in the linear regime, i.e. for 𝑉 𝑑 = 0.05 𝑉, and the saturated regime, i.e. for 𝑉 𝑑 = 1𝑉, as well as the output characteristics 𝐼 𝑑 = 𝑓(𝑉 𝑑 ), obtained by the numerical model against the ones obtained by PS simulations in the case of a constant mobility. Likewise, we notice an overall agreement between the numerical results and PS simulation results, with some inaccuracies that could be attributed to the relative error manifested in the approximation concerning the drift current integral in Eq 4.23, in addition to the pre-mentioned chosen approximations and smoothing functions.

Chapter 5: The analytical model

Note that, by rearranging our system of equations, we can define two residual error terms, which represent the residual error from each equation in the system as depicted in Eq 5.16: Doing so, after some additional rearrangements, the resulting system is composed of two equations that are both linear to the unknowns 𝜀 1 and 𝜀 2 , as depicted in Eq 5.17:

Eq 5.17

We also choose to label some terms from Eq 5.17 by using the designation described in Eq 5.18 and Eq 5.21, as such choice will allow us to write the analytical expressions of 𝜀 1 and 𝜀 2 in its most compact form: Indeed, the application of First-order Taylor expansion requires the respective derivative of the inversion charge density terms 𝑄 𝑖𝑛𝑣1,2 (𝑉 𝑠1,2 , 𝜙 𝑖𝑚𝑟𝑒𝑓 ) with respect to the surface potentials 𝑉 𝑆1,2 , which will be detailed ensuingly. Moreover, note that such procedure remains of electrostatic nature, therefore the variation of different quantities with respect to the quasi-Fermi level along the channel is disregarded.

Firstly, the expressions of 𝜕𝑄 𝑖𝑛𝑣1,2 (𝑉 𝑠1,2 ) 𝜕V s1,2 ⁄ can be directly given by Eq 5.21, which reveals that the computation of 𝜕𝑄 𝑖𝑛𝑣1,2 (𝑉 𝑠1,2 ) 𝜕V s1,2 ⁄ comes back to the computation of 𝜕∆V(Q g1,2 ) 𝜕V s1,2 ⁄ . Note that subsequently, for the sake of preserving the compact form of our equations at will, we designate the arguments of the exponential function, originally allocated to the front/back interfaces separately, as simply 𝐴𝑟𝑔 1 and 𝐴𝑟𝑔 2 .

𝜕V s1,2 )

Eq 5.20

Correspondingly, for the computation of d∆V(Q g1,2 ) dV s1,2 ⁄ , we can use the equivalence presented in Eq 5.21, considering that initially the quantum shifts are defined as functions of gate charge densities. As doing so we only need to get the derivatives of d∆V(Q g1,2 ) with respect to their plain variable Q g1,2 , then we multiply the out-come by the derivatives of Q g1,2 with respect to V s1,2 i.e. by dQ g1,2 dV s1,2 ⁄ = -C ox1,2 .

Eq 5.21

Consonantly, we call back the definition of ∆V(Q g1,2 ) :

Where the closed-form of the functions Q neg (Q g1,2 ) and Q pos (Q g1,2 ) are given by:

Eq 5.23

Accordingly, the derivative of the expression of ΔV(Q g1,2 ) presented in Eq 5.22 is:

Eq 5.24

With the respective derivatives of Q neg (Q g1,2 ) and Q pos (Q g1,2 ) given by Eq 5.25: Such expressions can be rearranged to be written as in Eq 5.26:

Thusly, by injecting the expressions given by Eq 5.26 into Eq 5.24 we will have the final closedform expression of d∆V(Q g1,2 ) dQ g1,2 ⁄ presented in Eq 5.27: q Moreover, in order to sustain the numerical stability of our model amid the drift drain current calculations, we introduce the concept of effective points, defined as, For the front and back subband potentials: 𝜒 s1,2 𝑒𝑓𝑓1 = min(𝜒 𝑠1,2 , 𝜒 s1,2 𝑖𝑛𝑡 ) , 𝜒 s1,2 𝑒𝑓𝑓2 = min(𝜒 𝑠1,2 , 𝜒 s1,2 𝑠𝑎𝑡 ) Eq 5.54

In a similar manner, we define the effective gate charge densities that will be useful for the computation of the drift integral in the case of a mobility function, following expressions in Eq 5.55.

𝑄 inv1,2 𝑒𝑓𝑓1 = min(𝑄 𝑖𝑛𝑣1,2 , 𝑄 inv1,2 𝑖𝑛𝑡 ) , 𝑄 inv1,2 𝑒𝑓𝑓2 = min(𝑄 𝑖𝑛𝑣1,2 , 𝑄 inv1,2 𝑠𝑎𝑡 )

Eq 5.55

Such effective points are chosen to ensure the logical coherence of the integral boundaries, in case 𝜒 s1,2 𝑖𝑛𝑡 or 𝜒 s1,2 𝑠𝑎𝑡 take irrational values, for instance, due to the denominator term in the expression of 𝜒 s1,2 𝑖𝑛𝑡 term, sometimes when 𝜆 1,2 𝑠𝑟𝑐 → 𝜆 1,2 𝑠𝑎𝑡 , the 𝜒 s1,2 𝑖𝑛𝑡 value might be found beyond the 𝜒 s1,2 𝑠𝑎𝑡 value, such configuration is unappealing. Accordingly, the use of such definition of effective points constitutes a good way to circumvent such numerical instability.

Thusly, the integration boundaries adopted formerly, i.e. between the source point and the intersection point for the first part of integration, and between the intersection point and the

Conclusion and perspectives