F Serra Di

Santa Maria

L Contamin

M Cassé

C Theodorou

Francis Balestra

G Ghibaudo
email: gerard.ghibaudo@grenoble-inp.fr.

F Serra Di

Comprehensive Kubo-Greenwood modelling of FDSOI MOS devices

Keywords: Kubo-Greenwood, mobility, modeling, MOSFET, FDSOI, cryogenic temperature

published or not. The documents may come   L'archive ouverte pluridisciplinaire

Introduction

The Cryogenic electronics is still a key research topic as permitting circuit performance improvements in terms of operation speed, turn-on behavior, thermal noise reduction, punch-through current decrease etc. [START_REF] Gutiérrez | Low temperature electronics: physics, devices, circuits, and applications[END_REF][START_REF] Balestra | Device and Circuit Cryogenic Operation for Low Temperature Electronics[END_REF][START_REF] Wada | Development of Low Power Cryogenic Readout Integrated Circuits Using Fully-Depleted-Silicon-on-Insulator CMOS Technology for Far-Infrared Image Sensors[END_REF][START_REF] Incandela | Characterization and Compact Modeling of Nanometer CMOS Transistors at Deep-Cryogenic Temperatures[END_REF][START_REF] Beckers | Characterization and Modeling of 28-nm Bulk CMOS Technology Down to 4.2 K[END_REF]. It finds application in high speed computing, sensing and detection, space electronics and more recently in readout CMOS electronics for quantum computing [START_REF] Hornibrook | Cryogenic control architecture for large-scale quantum computing[END_REF][START_REF] Maurand | A CMOS silicon spin qubit[END_REF]. In this framework, the characterization and modelling of MOSFETs down to cryogenic temperatures is still a key challenge. Moreover, Kubo-Greenwood formalism has been employed as a powerful tool for the modelling of transport in MOS inversion layers, enabling detailed mobility calculations [START_REF] Ghibaudo | Transport in the inversion layer of a MOS transistor. Use of Kubo-Greenwood formalism[END_REF][START_REF] Dura | Kubo-Greenwood approach for the calculation of mobility in gate-all-around nanowire metal-oxide-semiconductor field-effect transistors including screened remote Coulomb scattering-Comparison with experiment[END_REF][START_REF] Bonno | Effect of strain on the electron effective mobility in biaxially strained silicon inversion layers: An experimental and theoretical analysis via atomic force microscopy measurements and Kubo-Greenwood mobility calculations[END_REF] and recently subthreshold slope calculations [START_REF] Ghibaudo | On the modelling of temperature dependence of subthreshold swing in MOSFETs down to cryogenic temperature[END_REF].

In this work, we propose, for the first time, to use the Kubo-Greenwood approach for the drain current modelling as a function of gate voltage in 28nm FDSOI MOSFETs down to deep cryogenic temperatures. We first validate the single 2D subband approximation for the description of the MOS inversion charge with gate voltage down to liquid helium temperatures. Then, we show that the drain current transfer characteristics can be well modeled within the Kubo-Greenwood formalism down to very low temperatures, provided phonon, Coulomb, neutral and surface roughness scattering limited mobility laws are properly calibrated.

Experiments details

The measurements were completed on 28nm FDSOI MOSFETs with silicon film thickness tsi=7nm and buried oxide (BOX) thickness tbox=25nm from STMicroelectronics. NMOS transistors were processed from (100) handle substrate, with <100>-oriented channel, and a high-k/metal gate Gate-First architecture [START_REF] Planes | 28nm FDSOI technology platform for high-speed low-voltage digital applications[END_REF]. Low-Vth transistors were available with un-doped channel through a doped back plane (NWELL doping NA=10 18 cm -3 ) below the BOX. Thin gate oxide (with equivalent oxide thickness EOT=1.1nm) devices with gate length L varying from 30nm up to 10µm and with gate width W=1µm or 10µm were tested using a cryogenic probe station down to 4.2K.

The gate-to-channel capacitance characteristics Cgc(Vg) were measured with an HP 4284 LCR meter at 1MHz frequency and 10mV AC level using the standard split C-V configuration. The drain current Id(Vg) MOSFET transfer characteristics were recorded in linear region (Vd=30-50mV) with an HP4156 parameter analyzer. All the measurements were made at zero back bias. Note that the influence of back bias on the transfer characteristics has been studied in details elsewhere [START_REF] Paz | Front and back channels coupling and transport on 28 nm FD-SOI MOSFETs down to liquid-He temperature[END_REF] and goes beyond the scope of this work.

Kubo-Greenwood transport and FDSOI MOSFET modelling

The main equations used for the Kubo-Greenwood transport and FDSOI MOSFET modelling are reminded in the following sub-sections.

Kubo-Greenwood transport modelling

The inversion layer density n within a single 2D subband approximation is related to the Fermi level Ef by [START_REF] Ghibaudo | Transport in the inversion layer of a MOS transistor. Use of Kubo-Greenwood formalism[END_REF][START_REF] Ghibaudo | On the modelling of temperature dependence of subthreshold swing in MOSFETs down to cryogenic temperature[END_REF][START_REF] Mott | Electronic Processes In Non-crystalline Materials[END_REF],

𝑛(𝐸 𝑓 , 𝑇) = 𝑘𝑇. 𝐴 2𝑑 . 𝑙𝑛 (1 + 𝑒 𝐸 𝑓 -𝐸 𝑐 𝑘𝑇 ) (1)
where A2D=g.md*/(.ħ) is the 2D density of states with g the subband degeneracy factor, md* the DOS effective mass and ħ the reduced Planck constant and kT is the thermal energy, T being the temperature.

Note that Ef is referred to the subband edge Ec, stated here to zero.

The inversion layer sheet conductivity, , which is a function of Fermi level and temperature, can be computed by integration over energy E of the so-called energy conductivity function E(E) as [START_REF] Ghibaudo | Transport in the inversion layer of a MOS transistor. Use of Kubo-Greenwood formalism[END_REF][START_REF] Ghibaudo | On the modelling of temperature dependence of subthreshold swing in MOSFETs down to cryogenic temperature[END_REF][START_REF] Mott | Electronic Processes In Non-crystalline Materials[END_REF],

(𝐸 𝑓 , 𝑇) = ∫  𝐸 (𝐸) (- 𝜕𝑓 𝜕𝐸 ) 𝑑𝐸 +∞ 0 (2)
where

𝑓 = 1 (1 + 𝑒 𝐸-𝐸 𝑓 𝑘𝑇 ) ⁄
is the Fermi function and,

 𝐸 (𝐸) = 𝑞. 𝐸. 𝐴 2𝑑 . µ(𝐸) (3)
with E being the carrier kinetic energy and µ(E) the energy mobility function.

When several scattering take place in the electronic transport, the energy mobility function µ(E) can be evaluated by the Matthiessen rule after adding the various scattering rates for each energy such as,

µ(𝐸) = ( 1 µ 𝑝ℎ + 1 µ 𝑁 + 1 µ 𝐶 + 1 µ 𝑆𝑅 ) -1 (4)
which accounts for the phonon (µph), neutral (µN), Coulomb (µC) and surface roughness (µSR) limited mobility components, respectively.

Based on Monte Carlo quasi 2D simulations, the phonon limited mobility in an inversion layer has been established being essentially a function of temperature and vertical electric field as [START_REF] Gámiz | A comparison of models for phonon scattering in silicon inversion layers[END_REF],

µ 𝑝ℎ (𝑇, 𝐹) = 1180. [( 𝑇 300 ) 2.11 + ( 𝑇 300 ) 1.7 . ( 𝐹 𝐹 0
)

(𝑇) ] -1 (cm 2 /Vs) (5)
with F being the effective electric field (defined below), (𝑇) = 0. with the Kubo-Greenwood integral of Eq. ( 2) and plotted versus inversion layer density for various temperatures (see Fig. 1). As can be seen from the right figure, the Coulomb effective mobility merges the Boltzmann statistics limit at low carrier density, i.e. µCeffµC(2kT), whereas µCeff reaches the degenerate statistics asymptote at high carrier density, i.e. µCeffµC(Ef)=µC(n/A2d). As shown in Fig. 1 (dashed lines), these two asymptotic regions can be well approximated by the analytical closed-form expression,

µ 𝐶𝑒𝑓𝑓 (𝑛, 𝑇) = [µ 𝐶 (2𝑘𝑇) 𝑎 + µ 𝐶 ( 𝑛 𝐴 2𝑑
) 𝑎 ]

1 𝑎 [START_REF] Maurand | A CMOS silicon spin qubit[END_REF] where the exponent a1.5 has been adjusted to optimize the transition between the two regions. It could be noted that Eq. ( 7) can be regarded as the counterpart of Eqs ( 5) and ( 8) for phonon and surface roughness scatterings.

The surface roughness limited mobility (cm 2 /Vs) is mainly proportional to the square of the reciprocal effective electric field and is well approximated by the expression [START_REF] Villa | A physically-based model of the effective mobility in heavily-doped n-MOSFETs[END_REF]:

µ 𝑆𝑅 (𝑇, 𝐹) = 8.8×10 14 𝐹 2 . 𝑒𝑥𝑝 [-( 𝑇 850 ) 2 ] (cm 2 /Vs). (8)
Finally, the neutral scattering limited mobility is well known to be energy and temperature independent and simply reads [START_REF] Erginsoy | Neutral impurity scattering in semiconductors[END_REF],

µ 𝑁 = 𝐶𝑜𝑛𝑠𝑡𝑎𝑛𝑡 . ( 9
)
where the constant is proportional to the reciprocal neutral defect number. Degen.

T (K)=4.2K

Two options can now be used to calculate the inversion layer sheet conductivity  as a function of inversion layer density and temperature. First, one can rigorously compute the Kubo-Greenwood integral of Eq. ( 2) with the energy mobility law of Eq. ( 4) using the Coulomb scattering energy mobility law of Eq. ( 6) and Eqs ( 5), ( 8)-( 9) for the other scattering mechanisms. In this case, the Matthiessen rule is applied for each energy range in the integration process. Alternatively, the Kubo-Greenwood integral can be overlooked and the conductivity calculated using the explicit mobility laws versus carrier density, electric field and temperature using the Matthiessen rule applied to the mobility laws of Eqs ( 5), ( 7)-( 9) and the carrier density depending on the device biasing condition (see section 3.2).

FDSOI MOSFET modelling

For the modelling of an FDSOI MOSFET device, we consider that the 2D inversion layer is located at the front oxide/silicon channel interface, such that the gate charge conservation equation reads,

𝑉 𝑔 = 𝑉 𝑓𝑏 + 𝑉 𝑠 + 𝑄 𝑖 𝐶 𝑜𝑥 + 𝐶 𝑖𝑡 .(𝑉 𝑠 -𝑉 0 ) 𝐶 𝑜𝑥 + 𝐶 𝑏 .(𝑉 𝑠 -𝑉 𝑏 )
𝐶 𝑜𝑥 [START_REF] Bonno | Effect of strain on the electron effective mobility in biaxially strained silicon inversion layers: An experimental and theoretical analysis via atomic force microscopy measurements and Kubo-Greenwood mobility calculations[END_REF] where The surface potential Vs can be obtained by solving Eq. ( 10) for given front gate voltage and back bias. Thus, the Fermi level Ef can be calculated for any bias as:

𝐸 𝑓 (𝑉 𝑔 , 𝑉 𝑏 ) = 𝑞. [𝑉 𝑠 (𝑉 𝑔 , 𝑉 𝑏 ) -𝑉 0 ] ( 11
)
where V0 is a reference potential depending on channel doping level.

The effective electric field F entering the phonon and surface roughness limited mobility components of Eqs ( 5) and ( 8) is evaluated by the usual expression accounting for the inversion charge and back electric field as, 5), ( 6) and [START_REF] Ghibaudo | Transport in the inversion layer of a MOS transistor. Use of Kubo-Greenwood formalism[END_REF]. It should be mentioned that these values are close to the original ones in [START_REF] Ghibaudo | Transport in the inversion layer of a MOS transistor. Use of Kubo-Greenwood formalism[END_REF][START_REF] Gámiz | A comparison of models for phonon scattering in silicon inversion layers[END_REF][START_REF] Villa | A physically-based model of the effective mobility in heavily-doped n-MOSFETs[END_REF]. The value of µN, mentioned in Fig. 3

. ( 14
)

Results and discussion

Long channel devices

 Vg (V) Vg (V) Vg (V) Vg (V) Y (A 1/2 V 1/2 ) Id (A) gm (A/V) Id (A) T  T  T  T 
The variations of the effective mobility, µeff=/(qn), with inversion layer density deduced from the Kubo-Greenwood modeling are shown in Fig. 4a for various temperatures. They reveal the onset of a bell-shaped mobility behavior at very low temperatures as is usual [START_REF] Balestra | Device and Circuit Cryogenic Operation for Low Temperature Electronics[END_REF][START_REF] Ghibaudo | Transport in the inversion layer of a MOS transistor. Use of Kubo-Greenwood formalism[END_REF][START_REF] Stern | Calculated Temperature Dependence of Mobility in Silicon Inversion Layers[END_REF], due to the dominance of Coulomb and surface roughness scattering processes. Figure 4b displays the change in the temperature dependence of µeff for various carrier densities taken at weak, intermediate and strong inversion.

At low carrier density, µeff increases with temperature due to Coulomb scattering predominance (see Fig. 1) before to decrease at higher temperature due to enhanced phonon scattering. At high carrier density, µeff always decreases with temperature due to phonon diffusion supremacy. These features are also illustrated in Fig. 5, where the variations of µeff with carrier density are compared to the various mobility law components. As expected, at very low temperature, Coulomb scattering is dominating, whereas, at high temperature, phonon scattering is prevailing. In Fig. 6, we compare, for various temperatures, the variations of the effective mobility µeff with 2D carrier density as obtained from Kubo-Greenwood integral of Eq. ( 2) to those obtained using the Coulomb scattering mobility analytical approximation of Eq. ( 7). As can be seen, such an explicit analytical approximation for the Coulomb scattering mobility provides an overall good description of the global effective mobility, especially at low temperatures. For higher temperatures, a 10 to 20 % discrepancy appears at low carrier density, likely due to the error resulting from the use of the Matthiessen outside the Kubo-Greenwood integral in Eq. (2). However, the Coulomb scattering analytical mobility law of Eq. ( 7) could constitute a reasonably good approximation for a global effective mobility description useful for compact modelling purpose.

Short channel devices

In Figures 7 and8 are displayed the experimental and Kubo-Greenwood-modeled Id(Vg), gm(Vg) and Y(V g ) characteristics for MOS devices with gate length ranging from 30nm up to 1µm and for various temperature, showing again very good agreement. In this case, it should be noted that the best Kubo-Greenwood model fits have been obtained by keeping the same mobility parameters as for long devices (Fig. 3), at the exception of the neutral mobility component µN, and of the reference potential V0, which were adjusted versus gate length for each temperature as shown in Fig. 9. It is found that V0 follows the same trend vs L as the threshold voltage exhibiting a slight roll-off due to short channel effect (not shown). However, as can be seen from Fig. 9b, µN is strongly degraded as the channel length is reduced, regardless of temperature, due to the increased influence of neutral defects close to source and drain in agreement with the data of [START_REF] Shin | Low temperature characterization of mobility in 14nm FD-SOI CMOS devices under interface coupling conditions[END_REF] (square symbols in Fig. 9b) and with the theoretical analysis of [START_REF] Ghibaudo | Mobility characterization in advanced FD-SOI CMOS devices (part III pages[END_REF]. It should also be noted that, when fitting the transfer characteristics for short channel devices, the source/drain series resistance (Rsd) effects were taken care using the Y-function, which is independent of Rsd [START_REF] Ghibaudo | A new method for the extraction of MOSFET parameters[END_REF]. To this end, the mobility parameter µN was first tuned to adjust the experimental Y(Vg) curves. Then, Rsd was adjusted using Eq. ( 14) to fit the experimental Id(Vg) and gm(Vg) characteristics.

(V) Vg (V) Id (A) Id (A) Y (A 1/2 V 1/2 ) gm (A/V) L  L  L  L  Vg (V) Vg (V) T=300K
Typical values for Rsd were found in the range 230-250.µm for temperatures varying from 25K up to 300K.

Conclusions

A comprehensive Kubo-Greenwood modelling of FDSOI MOS devices has been performed down to deep cryogenic temperatures. Interestingly, a single set of mobility parameters was only needed to fit the data versus temperature for long channel devices. Instead, in short channel MOSFETs, the neutral scattering mobility component µN was found to be nearly temperature independent and to significantly decrease at small gate length due to the enhanced influence of neutral defects close to source/drain ends. Therefore, such a Kubo-Greenwood modelling can provide a physical insight into the scattering processes limiting the transport in FDSOI MOSFETs down to very low temperatures.

Besides, a closed-form expression for the Coulomb scattering mobility law as a function of temperature and 2D carrier density has been established and could be very useful for the development of compact models at cryogenic temperature.

Acknowledgment

The authors are grateful to EU H2020 RIA project SEQUENCE (Grant No. 871764) and to ERC Synergy QuCube (Grant No. 810504) for financial support.

Fig. 1 .

 1 Fig. 1. Variations of Coulomb scattering effective mobility µCeff with inversion layer density n for

  Qi(Ef,T)=q.n(Ef,T) is the absolute inversion charge, Vs is the front surface potential, Vfb is the flat band voltage, Vg is the front gate voltage, Vb is the back bias, Cox the front oxide capacitance, Cit (=q.Nit) is the front interface trap capacitance, Nit being the interface trap density, and, 𝐶 𝑏 = (𝐶 𝑠𝑖 𝐶 𝑏𝑜𝑥 ) (𝐶 𝑠𝑖 + 𝐶 𝑏𝑜𝑥 ) ⁄ is the substrate coupling capacitance, Csi being the silicon capacitance and Cbox the BOX oxide capacitance.

  Cgc(Vg) characteristics were measured on large area MOSFETs with W=L=10µm. The inversion charge was obtained by integration of the Ccg(Vg) curves starting from Vg=0 as is usual in split C-V technique. Experimental and modeled Cgc(Vg) and Qi(Vg) characteristics are shown in Fig. 2 for various temperatures. They clearly reveal the adequacy of the single subband approximation of Eq. (1) and of the MOSFET approach of Eq. (10) for the modelling of the gate charge control in FDSOI MOS transistors down to liquid helium temperature. The model parameters used in the simulation are indicated in the caption of Fig. 2. In Fig. 3 are reported the experimental and Kubo-Greenwood modeled Id(Vg), gm(Vg) and Y-function Y(V g ) = I d √g m ⁄ [19] characteristics for such a long channel device and for various temperatures, showing very good agreement. The best Kubo-Greenwood model fits have been achieved after proper calibration of the phonon, Coulomb and Surface roughness mobility law parameters, whose values are those displayed in Eqs (

Fig. 4 .

 4 Fig. 4. a) Variations of µeff with 2D carrier density n for various temperatures T(K)=4.2, 10, 20, 50,

Fig. 5 .

 5 Fig. 5. Variations of µeff (red solid line) and of various scattering mobility component with 2D carrier

Fig. 6 .

 6 Fig. 6. Variations of µeff with 2D carrier density as obtained from Kubo-Greenwood modeling (red

Fig. 7 .

 7 Fig. 7. Experimental (red solid lines) and modeled (blue dashed lines) IdVg), gm(Vg) and Y(Vg)

Fig. 8 .Fig. 9 .

 89 Fig. 8. Experimental (red solid lines) and modeled (blue dashed lines) IdVg) and gm(Vg) characteristics

  𝐹(𝑉 𝑔 , 𝑉 𝑏 ) = 𝑄 𝑖 (𝑉 𝑔 ,𝑉 𝑏 ,𝑇) 2 +𝐶 𝑏 .[𝑉 𝑠 (𝑉 𝑔 ,𝑉 𝑏 )-𝑉 𝑏 ]The drain current is then calculated in linear operation region, i.e. for small drain voltage Vd, as, 𝐼 𝑑 (𝑉 𝑔 , 𝑉 𝑑 , 𝑉 𝑏 , 𝑇, 𝐹) = 𝐼 𝑑 (𝑉 𝑔 , 𝑉 𝑑 , 𝑉 𝑏 , 𝑇, 𝐹) =

			 𝑠𝑖	(12)
	where si is the silicon permittivity.
	𝑊 𝐿	. [𝐸 𝑓 (𝑉 𝑔 , 𝑉 𝑏 ), 𝑇, 𝐹]. 𝑉 𝑑 .	(13)
	In short channel devices, the source/drain series resistance Rsd effect is accounted for through Ohm's
	law as,
		𝑊 𝐿	.[𝐸 𝑓 (𝑉 𝑔 ,𝑉 𝑏 ),𝑇,𝐹].𝑉 𝑑
		1+𝑅 𝑠𝑑 .

𝑊 𝐿

.[𝐸 𝑓 (𝑉 𝑔 ,𝑉 𝑏 ),𝑇,𝐹]

  caption, has been set to a large value for long channel device, indicating that the neutral scattering is negligible in this case.

	Fig. 2. Experimental (red solid lines) and modeled (blue dashed lines) Cgc(Vg) and Qi(Vg)
	characteristics for various temperatures T(K)=4.2, 10, 20, 50, 100, 150, 200, 250 and 300
	(W=L=10µm, model parameters: Cox=2.1 µF/cm 2 , Cbox=0.14 µF/cm 2 , Csi=1.52 µF/cm 2 , Cit=0.16
	µF/cm 2 , V0= 0.5 V).

Fig. 3. Experimental (red solid lines) and modeled (blue dashed lines) IdVg), gm(Vg) and Y(Vg) characteristics for various temperatures T(K)=4.2, 10, 20, 50, 100, 150, 200, 250 and 300 (µN=3000cm 2 /Vs, Vd=50mV, W=L=10µm, Nit=2-810 12 /eVcm 2 ).


Q Zeng

Pierre Navidi

Joseph Zarka

An optimal design of blasting

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Blasting is one of the basic operations in mining and quarrying industry. An efficient blasting for given ge ology conditions, rock volume to break and security conditions, results from a correct choice of explosive types, explosive quantity and the blasting planning.

One can think of diameter, and length of boreholes, drilling pattern, initiation modes, firing sequence, de lay time ... Efficiency is then determined by the qual ity of rock fragmentation and a minimum cost. In this sense it is an optimization problem. At present, al most all of the blast design is based on experiences, at fkrst because the behavior of the explosive and the sur rounding fractured rocks are still not well known, sec ondly because of the complex coupling effects. How ever, the computer aided design techniques, numeri cal blasting models and expert systems can consider ably reduce both the cost and duration of blast design while increasing the quality of prediction. This is the motivation of our research.

The key to a good blast design is a precise predic tion. The first part of our work focus on the predic tion by numerical simulations, in which the modelling is the discussion center.

The second part is on to the optimization of blast design. After introducing . the general methodology of our approach we show an application example for blast design.

Simulation of blasting

Many computational models of rock blasting a. re based on data obtained from experiments or/and em pirical relations derived from experimental data and analytical results for very simple cavities problems [START_REF] Sharpe | The production of elastic wave by explosive pressure, part itheory and empiri cal field observations[END_REF], [7], [START_REF] Digby | Computer simulation of blasking-induced vibration, frature and fragmen tation processes in brittle rock[END_REF].... They are either too much simplified or very speciflc. Recently, more fundamental models based on the material constitutive laws, usually im plemented in a finite element or finite difference code, have been attempted [START_REF] Grady | Continuum mod elling of explosive fracture in oil shale[END_REF],[10], [12] .... From a mechanical point of view, the blasting prob lem couples two difficult fields : the detonation of ex plosives and the rock mechanics. At present, the two fields are still in course of development. This section begins by a brief introduction of an existing 30 com puter code we used for our simulations. An explicit scheme algorithm enable us to handle the non-linear aspects and the fluid-structure coupling. Then a re view of the models for the explosive a. nd the rock is given. After,we introduce a new damaged rock consti tutive relation. Finally some numerical results based on this new model are presented to show its tractabil ity.

Computer code description

The program used can simulate most of the nonlin ear dynamic transient phenomena [START_REF] Mecalog | [END_REF]. This code can handle solid and fruid mechanics problem and the cou pling between them. The main features used for our specifmc problems are :

• explicit central-difference time-integration algo rithm

• Arbitrary Lagrangian-Eulerian (ALE) formula tion [START_REF] Hirt | An arbi trary Jagrangian-eulerian computing method for all fqow speeds[END_REF] for fluid/structure interaction

• artificial viscosity in order to smear shock fronts [START_REF] Richtmayer | A method fsr the numerical calculation of hydrodynamical shock[END_REF] • Incremental form of constitutive relations

Detonation model

Computation of the detonation propagation in explo sives is a complex problem. It results from strong chemico-physics interaction during the detonation process[l4),[6) .... When a coupling with a structure is considered, the numerical modeling of the detona tion is often simplified. This is because the response of the structure is the main interest. Such models are often of continuum type and admit the following hypotheses :

1. reduction of the reactive fluid in a fluid with two components, an initial substance and a final sub stance (explosive reaction products);

2. each of the two substances is separately in thermo-chemical equilibrium;

3. quasi Chapman-Jouguet detonation a. quasisonic detonation velocity).

These hypotheses reduce the complex computation to a simple fluid flow one with discontinuities which a.re then often smeared as follows :

• an artificial viscosity is used for the thermophys ical variables of discontinuity (shock wave);

• a. decomposition law is used for the thermochem ical variables of discontinuity.

Therefore, the behavior of the explosive is repre• sented by an equation of state and a. decomposition law.

Often, in case of the detonation/structure simu lation, the equation of state has a pseudo-potential form p( v, E) (pressure-volume-energy) to avoid the thermo-chemical iterative calculation. The decompo sition law is of �bulk burn» type due to the homo geneity consideration. A variety of explosive behavior ranging from the simple to complex models are avail able, we use the following models[30):

• The Jones-Wilkins-Lee (JWL) equation of state for the explosive products pressure :

p = A(l -�)e-R ,V + B(l -�)e-R2V + � R , v R21 v V where V= ..!!. ""
relative volume explosive products pressure detonation energy per unit volume constant parameters

• if A = B = 0 , we refind the simple �gamma law» with w =, -i (1 = r = �)

• the �C.-J volume burn» chemical decomposi tion law :

m-1-V -1-Vc;
where the burn fraction m is the indicator of de composition progress.

• the current pressure is calculated by :

P=m•p
The detonation velocity, the pressure at Chapman-Jouguet point PcJ and the initial en ergy Eo should also be known.

After having done several simulations (explosion in water, explosion in an elastic medium ), we noticed some interesting results [START_REF] Zeng | Optimisation et /'utilisation des ex p/osifs en Genie civil[END_REF] :

• computation with the equation of state �gamma law» has a similar result compared to the com putation with the equation of state JWL.

• the explosive energy partition (shock en ergy /bubble energy) can be characterized by the coefficient / of the �gamma law», the greater 'Y is, the greater is the shock energy.

• the modeling of blast loading by using defined pressure-time history is only a crude approxima tion specially for the case of the point initiation.

Rock mechanical model

Another essential point for the blasting simulation is a proper representation of the rock behavior. The ap proach which involves the formulation of appropriate material models (constitutive relation) is preferred.

The main phenomena that must be accounted for is crack initiation and propagation. Two alternatives are then possible. One is the discrete model which takes into account the physical presence of the cracks in a direct approach(ll), [START_REF] Munjiza | On a rational approche to rock blasting[END_REF] .... For the computation, this one is often very costly and suitable for the post period of the blasting processus. Another way is to take a continuum damage model. We have chosen the last method which is more consistent with a fi nite element approach. A continuum model [START_REF] Zeng | Optimisation et /'utilisation des ex p/osifs en Genie civil[END_REF] is therefore developed and implemented into the explicit code. This is an elastic anisotropically damaging be havior based on the damage mechanics theory intro duced by Kachanov [START_REF] Mark | Time of the rupture pro cess under creep condition[END_REF]. The constitutive relation is :

u=C(D):£
where C is the elastic effective stiffness tensor and D is the damage tensor.

In an incremental form : this concept is deduced from the impact experi ment results [9)-

2.3.l Damage description

The nature of the rock flaws and the possible complex loadings require an anisotropic damage description.

Generally, an eight order tensor should be introduced to describe the anisotropic damage, but the difficulty to identify all the parameters is obvious. It is thus desirable to reduce the order of that tensor. Based on the third hypothesis above, a symmetric second-order damage tensor D is introduced, reducing the model to an orthotropic one.

We can consider each of the eigenvalues D; as a group of paraUel cracks perpendicular to the associ ated eigenvector. As if :

4 D; = -7r • f 3 ( 2 )
where f = Na•c� is the cracks (Na cracks ofradius Ca) group density . This definition is the same as the one for a homog enization method [START_REF] Nemat-Nasser | Micromechanics overall properties of heterogeneous materials[END_REF], [START_REF] Hoenig | Elastic mduli of non-randomly cracked body[END_REF] so that the results of ho mogenization can be used to determine the effective elastic matrix which is going to be shown lately.

D; can be also defined with the main parameters as introduced by Grady and Kipp [START_REF] Grady | Continuum mod elling of explosive fracture in oil shale[END_REF] :

D; = N v ; with v = �7r • c3 ( c and N are mean values ) 2.3.2 Damage evolution law
Generally, the second thermodynamic principle is used to determine the damage evolution Jaw as in an elastoplastic problem. However the dissipative energy is not easily established. Grady and Kipp have devel oped an evolution law based on the activation and growth of an initial Weibull distribution of fracture producing foaws, ftr an isotropic behavior [START_REF] Grady | Continuum mod elling of explosive fracture in oil shale[END_REF]. It can be extended to an orthotropic material : { I s,,cm+3l, j 1 / 3 ( u�-u. The damage criterion is usually derived from the Grif fith theory [START_REF] Griffith | The phenomenon of rupture dans fnow in solids[END_REF], [START_REF] Bui | Mecanique de la Rupture Fragile[END_REF], with the damage evolution law (3) shown above, we just take the following simple form :

G f >Gu

2.3-4 Effective elastic matrix

Based on the results of the homogenization methods [START_REF] Nemat-Nasser | Micromechanics overall properties of heterogeneous materials[END_REF], [START_REF] Hoenig | Elastic mduli of non-randomly cracked body[END_REF], the symmetric effective elastic matrix takes the form:

_ -[ [ C(l )l 0 ] [c] -o [c< 2 >] C -(1 ) -fi( 2 d d )d C -(1) -E 2 11 -a 1 -I/ 2 3 I 22 -c;(l -I/ d3di )d2 C-(1) -E ( 2d d )d C -(IJ E 33 -o; 1 -v 1 2 3 12 = o;v(l + vd 3 )d1d2 C "(l) E ( -(I) E 1 3 = o; v 1 + vd2)d1rl3 C23 = ;;v(l + vdi)d3d2 C( 2 ) _ 2µ (1-PD2)(!-/iD3) 11 - 2-/JD,-/iD3 C (2) _ 2µ (1-11Di)( l -/iD3) 22 - 2-{JD, -1iD3 C (2) _ 2µ (1-liD2)(l-1iD1)
33 - (i=l,3)

a= 1 -v 2 (d1d 2 + d 2d3 + d1d3 + 2vd 2 d3d 1 ) (J = (1 -v)(2 -v)
v: poisson's ratio The result is in 0(/ 2 ) order to the constant strain or constant stress method and the formulae are easy to calculate.

2.3.5

Incremental algorithm

During each time step the following algorithm is used for the determination of the stress and the damage tensor values : given u(t), e:(t), D(t), i(t + dt) : Grady and Kipp which is strain-rate dependant. For sake of simplicity the formula is now fracture stress dependant :

D(t + dt) = D(t) + D(u(t)) dC = C(D(t + dt))::: C(D(t)) _ u(t + dt) = u(t) + C(t) : i(t + dt)dt + dC : e:(t)
( (!c )-T LM = 1J(k, m) E(m + 3)
where 1) is a constant and qc the tensile fyacture stress.

Simulation results

The damage evolution law based on the Weibull dis• tribution leads to predictions which are in good agree ment with the strain-rate effects discussed by differ ent authors[l6], [START_REF] Margolin | Numerical simula tion of fracture[END_REF][10] .... ). In order to show the ca pability of the model proposed to take into account anisotropic effects, the following numerical experi ment is done : we consider several cases of constant strain rate loading on a single element. The first load ing is carried out up to a value of D1 = 0.4 (largest eigenvalue of the damage tensor). After changing the loading direction, a second loading is applied to the same specimen, the curves in dot lines correspond the results for various directions of the second loading re spectively 0°, 30°, 45°, 90° Fig(!). shown at time t = 0.25ms et t = 0.5ms. This result underlines that :

• the damage begins around the cavity and pro gresses with decreasing value.

• a second damage area appears when the stress wave reaches the free surface and is refpected.

• the eigenvector directions of D tells us that the The results in Fig 7) and (8 are given for the largest eigenvalues of the damage tensor at t = l.5ms where the damage results are convergent.

We observe that the detonation with delay modi• lies the damage distribution and increases the damage volume (10% for this example). The space between the explosive and the surround ing rock is occupied by air. the decoupling ratio (bore hole diameter to explosive diameter) is 1.2, the con finement is perfect (the stemming is also occupied by rock). A point detonator is put at the bottom of the borehole. The rock is the same as in the first example [START_REF] Bui | Mecanique de la Rupture Fragile[END_REF].

In Fig ( 10), the damage result ( largest eigenvalue of the damage tensor) is shown at t = 0.375 ms, I . . . . The blasting simulation based on the advanced mod eling allows us to predict more generally and more precisely the blasting results. It can be used also to better understand fragmentation mechanisms, the roles played by the different design parameters [START_REF] Zeng | Optimisation et /'utilisation des ex p/osifs en Genie civil[END_REF]. However, It can not be yet a perfect tool for the blast design mainly because:

• it is difficult to simulate the entire blasting pro cess with a single model often much simplified.

• the strong non linearity and the local/global problems in blasting make the computation too expensive. Even if it is technically possible to simulate the entire blasting process by more so phisticated modeling (e.g.: a combined finite ele ment/ discrete element model [START_REF] Munjiza | On a rational approche to rock blasting[END_REF]), it is practically impossible to simulate the real sequence blasting problem.

In this context, the use of the automatic learning technics for searching rules (or models) seems to be a rational solution, and opens a realistic way to the blasting optimization.

In this section we will first give a brief introduction to the automatic learning technics. The methodology of blasting optimization is then discussed. Finally an example of automatic learning and optimization are given.

Automatic learning

In many technical domains, experts are aware that, very often, there is no a complete solution for a real world problem. But there is a possibility to build a data base of examples. Each sample of this base is obtained experimentally or numerically with some times some fuzzy or missing informations. The au tomatic learning £xpert Systems generator (l£S), developed at LMS during 1986-1990, can build a set of rules based on the description of such a examples base. Several applications have been shown [START_REF] Schoenauer | Incremental learn ing of rules and meta-rules[END_REF], [START_REF] Navidi | Clever optimal design of materials and structures[END_REF], [START_REF] Michel Te Rrien | Systemes experts par apprentis sage en contrfle non-destructif[END_REF], [START_REF] Hablot | Construction de solutions exactes en elastro-plasticite[END_REF], [START_REF] Bulik | Sur /'Optimisation de la Protec tion Parasismique[END_REF].... The main difficulty is to provide a good description of the examples base. The initial descriptors (let us call descriptor each single field of the data base) appear then as variables for the final model. Suppose that the expert defjnes a set of these ini tial descriptors describing the observed and/ or mea surable phenomena, (x; ) generally in a limited num ber (e.g. the radius of the borehole, explosive energy etc ... ). Indeed only a subset of this initial set is suf ficient to make the describtion of one example but it is not known a priori. These descriptors are then integrated within our actual limited knowledge to de termine a set of more complex ones which are the intef¦gent descriptors, (Xi) (e.g. results of a numeri cal simulation or the ftime-history of the radius of the bubble in the water blast test,etc ... ). These descrip tors are function of the initial descriptors and some of them may be considered among the conclusions (Ck) and then also be learned ( e.g. the damaged rock volume).

In order to build the f{ll knowledge of the problem or the model, the expert must prepare a fi le containing for ea. eh example :

The input descriptors: X1 I ....

I X6 I X1 I •••• • I X40 I • ••••• I
The output descriptors:

C1 I C2 I ...... . The symbolical option of C£S will give a set of rules while its numerical option will generates a polynomial based model which is often easier to handle.

Generating rules

For a simple problem, the data base, with the initial descriptors (x;) and the conclusions C1, can be pre• sented to the system. The system generates then a mathematical expression based on a polynomial ex pression :

C1 =POL (xi,x2,• • •, Xn )
For more complex problems, it is better to use the data base with the intelligent descriptors (Xi) in order to help the system to find a more reliable model fwr the conclusion C 1 :

C1 =POL (X1,X2,•••, X k)•
Once the model is created, new unknown examples may be introduced for evaluating the conclusions.

The building of the set of (X;) requires making them intrinsic in order to be able to reproduce them for various types of problems independently of the nature or type of problem considered.

Optimization

Many optimization technics are available. But we have the following points in mind :

• Design variables can be discrete and also take their value in a finite set.

• From a technological point of view many solu tions can lead to the same type of results. This means that if C1 is the cost function then in the design variable space POL can have flat domains ( with no gradient). Many optima can exist si multaneously that is POL does not define a con vex .

• Constraints generally varies depending on the do main in which the design set is looked for.

Blasting optimization by the traditional methods based on the mathematical programming theory such as quasi-Newton, conjugate gradient ... seems to be impossible. On the other hand, the evolution methods such as the genetic algorithms are more suitable for a general optimization problem.

Contrary to the traditional methods based on the continuum concept, the genetic algorithms make evolve a population of solutions with a manner sug gested by the genetic. Their advantages compared with the traditional methods are as following:

• the probabilistic characteristics can handle local optima.

• there is no other restriction for the functions (cost functions, constraint function), excepted to be calculable.

• the algorithm is extremely general.

• design variables can be of any type.

The main disadvantage is costly computation time.

Optimization by coupling the genetic algorithms with automatic learning is then a good strategy be cause the generation of the specific rules (for the cost function, the constraint functions) by automatic learning allows to compensate the handicap of the ge netic algorithms method.

It is interesting to observe that the fact of express ing all the relations in a simple polynomial way enable us to carry out genetic optimization very effectively.

The other intereseting points of this approach can be summarized as follows:

• It gives us a way of dealing with inverse problems, that cannot be handled by classical schemes.

• It lets us understand the essential descriptors which drives the desired conclusion. This can be a starting point for experts to build new PHYS ICAL or MATHEMATICAL MODELS.

• By a clever choice of the intelligent descriptors, it is possible to extrapolate models resulting from a reduced set of simple samples to more complex situations.

Generation a damaged volume rule by as

In this section we use the above methodolgy in a sim ple blasting design problem. Suppose that the explo sive type and the rock site are given. We search a damaged volume Vd (broken rock) expression in func tion of the drilling pattern parameters in order to carry out next the optimization calculation. This means that their are at least two local so lution. In fact the solution found depends highly on the starting point.

• Genetic Algorithm

With several different starting populations, the following solution is found : {d = 0.04, B = 0.5, Rd= 1.094118} => J = 1187.1654 Now, suppose that borehole diameter d and the ex plosive cartridge diameter can only take some fixed discrete values : One can notice that these points where not in the ini tial database. These values can then be used fvr a new numerical simulation (or experiment). The dam aged volume obtained is then compared with the one predicted by the model. If they agree, the model gen erated by the lt:S is consistent. Otherwise we can add a new example to our data base and build a new model. This work can be done iteratively until the model is consistent.

Conclusion

In this paper, it has been shown that the modelling of micro-craks by a continuum damage model may im prove the quality of the prediction. For the continum rock model considered, the mixed phenomenological -microphysical approach seems to be tractable; the microphysical approach allows to build a consistent model thus limiting the number of parameters, it also allows quantitative results for the fragmentation.The phenomenological approach simplify the model and decrease the computing time.

However, even if the blast prediction by simula tion appears to be a powerfzll prediction tool, it can only be an element of the industrial blasting design. Indeed, many other relevant parameters and techni cal aspects cannot directly enter the mathematical model. The use of the automatic learning lets us overcome this step efficiently, as we showed it on a simple example. The coupling between genetic algo rithms and automatic learning seems to be a powerful potential approach to blasting optimization.

  du= C(D): d£+dC(D): £ (1) The ma.in hypotheses are the following: 1. the dominant mode of failure is mode I as proven by numerous experiments [21],[22),(15) ....

2 .

 2 the material non linearity is totally governed by damage and thus plastic deformations are disre garded 3. The expanding plan directions of the penny shaped cracks are perpendicular to the tensile principal stress directions.

  i. � D; = l ( m+l){m+2) �gn --'--r"" ) 3 D; with i = 1,3 where n(x)=kx m G f >Gu else (3) Cg cr f, uu E constant fracture growth velocity principal stress and a constant depending on the material cohesion Young's modulus k,m constant Weibull parameters 2.3.3 Damage criterion

  2-/JD2-liD1 cU> = o c i ;e j) with d;=l-D;

2. 3 . 6

 36 Dominant fragment size 20.0 .. ..-----�----�----� 15.0 --first loading a-90 The calculation of the dominant fragment size (frag-� ment size corresponding to the largest volume fxaction of material) is also derived from the result given by 5.o

2. 4 . 1

 41 Calculations in plane strain with an imposed pressure function We use the model presented above to carry out the analysis on various geometries. One of them is a cir cular cylindrical cavity problem. Due to symmetry a quarter of the geometry is considered. The configura tion and the loading, a sine function pressure imposed inside the cavity, are given in Fig 2 and 3 : The surrounding rock is an oil shale which has the properties : Fig 4 and 5, the distributions of the largest principal value of the damage tensor are respectively I 0.0020 0.0040 strain in the loading direction

Figure 1 :

 1 Figure 1: Anistropic effect

Figure 2 :

 2 Figure 2: Configuration

Figure 3 :

 3 Figure 3: Pressure function

Figure 4 :

 4 Figure 4: First principal value of D at t = 0.25ms I . . . . ....... ::Jll" 1.1nn '•-'•" .... ..... ,.::m" l . IUU '•"*" ..

Figure 5 :Figure 6 :

 56 Figure 5: First principal value of D at t = 0.5ms

  Figure 7: simultaneous detonation kuPrlu1,1I I . . .. . ..... . O . llU1 . _ ,,. O.UJU O.UU7 Cl.U O.UUJ c.::u1 Cl.HI O.IUU ::um

Figure 8 :

 8 Figure 8: with a delay time r = 0.6ms 2.4.2 3-D calculation coupled with detona tion The configuration is shown in Fig (9) :

Figure 9 :

 9 Figure 9: Configuration in 3-D The explosive is idealized. Its brisance is near the ANFO. The equation of state used is the �gamma law� where the parameters a.re as follows : I Po(�) 2.8 0.9

  Figure 10: damage result in 3-D

  The database is built by numerically simulating blasting in plane strain. The choosen descriptors are borehole diameter d ; burden B;decoupf¥ng ratio Rd and two ratios d/ Rd , d/ B. The result in Ta b(l) is given by eleven simulations d(m) B(m) 0.5066 1.677 0.79336 0.0298 0.0987 0.050 0.600 1.172 1.467875 0.0427 0.0833 0.050 0.600 1.667 0.675639 0.0298 0.0833 0.040 0.500 1.3333 0.94768 0.03 0.08 0.040 0.600 1.333 0.871909 0.03 0.0667 0.100 0.750 1.257 1.520870 0.0796 0.133 0.100 0.750 1.474 1.265204 0.0678 0.133 0.200 1.000 1.2067 1.570654 0.1657 0.2 0.200 0.750 1.2067 2.957848 0.1657 0.267

  •d-2.44B Rd +(24.134-45.6Rd+23.l 7 R�)B (5 ) (and em= 0.00477.)

  The problem can be set as fullowing :For a given explosive type and rock in-site, we search the design parameters values (d,Rd,B) which maxim ize the specific chargeE� = Ve (explos�vevolwnc) •Taking back the result (equation (5) ) divided by (f.) 2 which is in proportion to the explosive volume, the problem becomes : maxi= 4.256B ll;f -88.95B + (90.39 -92.644Rd + 2 R' 25.54Rd):ii'j with the constraints : 1.0 �Rd� 2.0 0.4 � B � 1.0 0.5 < d < 0.2 d/ R � 2'. -0.03 (the explosive cartridge diameter is supposed superior to 0.03m) This "simple" problem is firstly resolved by one of the traditional methods, then by the genetic al gorithm. • Conjugate Gradient With a starting point : {d = 0.04, B = 0.8, Rd= 1.0}, a solution has be found such that : {d = 0.04, B = 0.5, Rd= 1.09358} => J = 1187.1658 With another starting point :{ d = 0.1, B = 0.8, Rd= 2.0} ,a different solution has be found such that : { d = 0.06, B = 0.5, Rd= 1.09358} => I= 995.79

  048, 0.05, 0.066, 0.084, 0.1, 0.15, 0.2} E {0.03, 0.036, 0.04, 0.05, 0.64} The solution found by genetic algorithm and chang ing Rd by d/ Rd is: { d = 0.048, B = 0.5, d/ Rd = 0.045 -Rd = 1.0667} => J = 980.96

Table 1 :

 1 damaged volume database First we ask the [,[S to build a polynomial model. Limiting the system to terms of order 2 with cross terms. the ££S output is :


Francesco Serra

Santa Maria

Christoforos Theodorou

Francis Balestra

Gerard Ghibaudo

Eunjung Cha

Cezar B Zota

In-depth electrical characterization of deca-nanometer InGaAs MOSFET down to deep cryogenic temperatures

Keywords: III-V, MOSFET, cryogenic, Y-function, characterization, InGaAs

This work presents a detailed electrical characterization of planar InGaAs on Insulator MOSFETs from room temperature (namely 300 K) down to deep cryogenic temperatures (10 K). The main electrical parameters of MOSFET operation (threshold voltage Vt, low-field mobility μ0, and subthreshold swing, SS) were extracted in both linear and saturation regions of operation through the consolidated Y-function method, for gate lengths down to 10 nm. The extracted parameters are first analyzed versus temperature and length and then compared against a more mature technology such as Silicon FDSOI MOSFETs. The results reveal competing advantages of the III-V alloy, particularly when going down to cryogenic temperatures.

I. INTRODUCTION

As the research on Quantum Computing (QC) advances , the engineering of the circuitry surrounding the quantum bits themselves requires further study. This necessity is better understood when considering the readout electronics [START_REF] Hornibrook | Cryogenic Control Architecture for Large-Scale Quantum Computing[END_REF] devoted to quantum computing: in order to minimize signal transmission delays and noise amplification due to different temperatures amongst signal stages, readout electronics in QC have to be as close as possible to the actual bits, therefore operating in the temperature range of a few units of Kelvins. Thanks to their enhanced mobility, III-V MOS devices can provide the same ON current at lower power supply voltages, and by turn reduced power consumption and heat dissipation, crucial at QC operation temperature. Thus, the precise identification of their electrical parameters behavior with temperature and channel length is required for reliable modeling and circuit design. On the other hand, as this technology is not as mature as its Silicon-based counterparts (Bulk, FDSOI, and FinFETs), the need of full electrical characterization down to deep cryo-temperatures, becomes challenging and critical, particularly in view of the emerging technology of QC.

To date, Silicon technology has been extensively studied and experimentally characterized [START_REF] Beckers | Theoretical Limit of Low Temperature Subthreshold Swing in Field-Effect Transistors[END_REF][START_REF] Diouf | Y function' method applied to saturation regime: Apparent saturation mobility and saturation velocity extraction[END_REF][START_REF] Liu | On the Interpretation of Ballistic Injection Velocity in Deeply Scaled MOSFETs[END_REF][START_REF] Di Santa Maria | Lambert-W function-based parameter extraction for FDSOI MOSFETs down to deep cryogenic temperatures[END_REF], while on the other hand many modelling and simulation studies have been carried out concerning III-V alloys, particularly highlighting the conduction in satellite valleys [START_REF] O'regan | Modeling the capacitance-voltage response of In0.53Ga0.47As metal-oxidesemiconductor structures: Charge quantization and nonparabolic corrections[END_REF].In this work, we study InGaAs MOS devices so as to compare the effects due to both channel length and temperature reduction with those already known for Si [START_REF] Di Santa Maria | Lambert-W function-based parameter extraction for FDSOI MOSFETs down to deep cryogenic temperatures[END_REF], particularly to see if the same physical limits apply as in the case of subthreshold swing settling [START_REF] Beckers | Theoretical Limit of Low Temperature Subthreshold Swing in Field-Effect Transistors[END_REF]. Moreover, we tried to position our characterization work with respect to results from transport simulation and to see how these impact the ease of parameter extraction.

II. DEVICES UNDER STUDY AND EXTRACTION METHODOLOGY

The devices under test were fabricated by IBM Zurich [START_REF] Zota | Effects of Post Metallization Annealing on InGaAs-on-Insulator MOSFETs on Si[END_REF], based on a III-V on insulator tehnology, incorporating a 20 nm InGaAs film, insulated by a buried oxied and handled through a Silicon substrate (Fig. 1). They share a common gate width value of W = 1 m and channel lengths spanning from L = 300 down to 10 nm. The measurements were carried out at wafer level down to 10K using a HP 4155A parameter analyzer and a Formfactor SussTech 300mm Cryo probe station.

In order to extract the main parameters that define the MOSFET performance, we utilized the Y-function method, expressed through Y(𝑉𝑔) = 𝐼𝑑/√𝑔𝑚  √β(Vg -Vt) [START_REF] Ghibaudo | New method for the extraction of MOSFET parameters[END_REF], where β = Wμ0CoxVd / L. This function suppresses any access resistance effect and thanks to its linear behavior versus Vg at Fig. 1. Schematics of studied III-V MOSFETs. Courtesy of IBM Zurich [START_REF] Zota | Effects of Post Metallization Annealing on InGaAs-on-Insulator MOSFETs on Si[END_REF] strong inversion, the linear fit can provide the threshold voltage, Vt, as well as the low-field mobility, μ0, value through the x-axis intercept and slope, respectively.

III. EXPERIMENTAL RESULTS AND DISCUSSION

A. Preliminary observations

The measured Id -Vg input characteristics in linear regime are plotted in Fig. 2a, for the two temperature limits (10 K and 300 K). Through a first glance, we can already observe how the threshold voltage, Vt, is affected by both channel length and temperature: on one hand, for a shorter channel Vt shifts downwards ("Vt roll-off" short channel effect), while on the other hand, it increases for a decreasing temperature.

Moreover, we can notice some kind of a shoulder/hump in the high Vg region, in both the Id-Vg and the Y-Vg (Fig. 2b) curves: this effect is most likely to be related to the onset of conduction in the L valley [START_REF] O'regan | Modeling the capacitance-voltage response of In0.53Ga0.47As metal-oxidesemiconductor structures: Charge quantization and nonparabolic corrections[END_REF] of III-V. What is most important for our study, however, is that this hump does not impact the overall linear behavior of the Y function in strong inversion, allowing us therefore an easy and reliable extraction of the main operation parameters.

B. Analysis of extracted parameters in ohmic region

When going down to low temperatures (LT), several effects are taking place. First, due to the Boltzmann statistics, as the temperature decreases, fewer and fewer electrons are promoted to the conduction band for the same gate voltage . This is reflected in an increase in threshold voltage as going towards lower temperatures, as shown in the extracted Vt values, plotted versus temperature for all measured gate lengths in Fig. 3. This increase at low temperatures reaches a maximum Vt value, which corresponds to the degenerated behavior of the semiconductor at that range of temperatures .

Moreover, the transition between OFF and ON states becomes sharper (Fig. 4), yielding a consequent decrease of the subtreshold swing, SW, for low temperatures, down to a settling value of about 10 mV/dec, which is attributed to the exponential band tails of states [START_REF] Beckers | Theoretical Limit of Low Temperature Subthreshold Swing in Field-Effect Transistors[END_REF]. This behavior of SW is confirmed for the measured III-V devices, as shown in Fig. 5, where the extracted values of SW are plotted versus temperature for certain channel lengths. Conversely, as the temperature lowers, there is less and less vibrational energy in the lattice, causing a decrease in phonon scattering , which allows an increase in the low-field mobility, μ0, as shown in a, where the extracted values of μ0 are plotted versus temperature for all measured gate channel length. This increase is evidently more pronounced in longer channel devices, whereas, in shorter channels, defect (neutral impurity) scattering is prevailing due to source/drain regions proximity , yielding a generally lower mobility and also a less significant increase at lower temperatures [START_REF] Di Santa Maria | Lambert-W function-based parameter extraction for FDSOI MOSFETs down to deep cryogenic temperatures[END_REF].

C. Behavior in saturation region

As we proceed towards the analysis of the device parameters in saturation region of operation (Vd = 1 V), we observe that both 10 and 300 nm length devices, plotted in Fig. 7, show no significant variation as going to low temperatures.

From the consequent extraction, the variation of both Vt and SW parameters with respect to temperature in the saturation region is very small. On the contrary, we notice that in saturation region, the curvature created by the takeover of L valley is no longer visible at low temperature. This can be explained when considering that the drain current is obtained by integration along the channel from source to drain [START_REF] Taur | Fundaments of Modern VLSI Devices[END_REF]: close to the drain region the influence of high Vd does not allow the quasi Fermi level to fill the L valley, thus attenuating its effect.

When extracting the μ0 values through Y function, plotted in Fig. 8, compared to the linear region, the devices present a μ0 reduction due to velocity saturation (vsat) effect [START_REF] Diouf | Y function' method applied to saturation regime: Apparent saturation mobility and saturation velocity extraction[END_REF]. Moreover, it is worth noticing how the extracted mobility never reaches the ballistic limit [START_REF] Liu | On the Interpretation of Ballistic Injection Velocity in Deeply Scaled MOSFETs[END_REF], revealing that the exponential behaviour of μ0 with respect to channel length is in fact scattering-related. The extracted values of vsat in saturation region are plotted in Fig. 9, showing a good stability with respect to temperature and a slight increase for a decreasing channel length due to overshoot effect .

D. Benchmarking III-V against Si for Cryo

Fig. 10 illustrates the extracted parameters of μ0, SW and vsat along with respective extractions done for Si [START_REF] Di Santa Maria | Lambert-W function-based parameter extraction for FDSOI MOSFETs down to deep cryogenic temperatures[END_REF] channel Fully Depleted (FD) SOI MOSFETs. When comparing the R&D-level III-V devices to this more mature and industrial Si technology like FDSOI (Fig. 10) [START_REF] Di Santa Maria | Lambert-W function-based parameter extraction for FDSOI MOSFETs down to deep cryogenic temperatures[END_REF], we note that although SW is much higher for III-V at 300K, both technologies reach the lowest limit value at 10K. Moreover, despite the higher interface trap density at the Al2O3 interface, III-V shows higher 0 and vsat at all temperatures, revealing a great III-V potential for cryo-related applications with further technology developments and the chance to overperform Si FDSOI in certain cases.

IV. CONCLUSION

We have performed a detailed electrical of planar InGaAs on Insulator MOSFETs from room temperature down to deep cryogenic temperatures (10K). The main MOSFET parameters (threshold voltage Vt, low-field mobility μ0, and subthreshold swing, SS) were extracted in linear region of operation using the consolidated Y-function method for gate lengths down to 10 nm, despite the possible presence of L-valley conduction. The saturation velocity was also extracted and analysed for all lengths and temperatures. The extracted parameters of the III-V devices follow the expected behavior with temperature as in Si, while demonstrating competing advantages as compared to Si MOSFETs, particularly when going down to cryogenic temperatures. and Si [START_REF] Di Santa Maria | Lambert-W function-based parameter extraction for FDSOI MOSFETs down to deep cryogenic temperatures[END_REF]

Fig. 2 .

 2 Fig. 2. Drain current (a) and Y-function (b) versus gate voltage at Vd = 30 mV for L = 10 and 300 nm at T = 10 and 300K. The dashed lines represent the fitting of the plots, based on Y function and linear regression respectively.

Fig. 3 .

 3 Fig. 3. Extracted threshold voltage versus temperature for Vd = 30 mV for various channel lengths

Fig. 4 .

 4 Fig. 4. Drain current versus gate voltage at Vd = 30 mV for L = 10 and 300 nm at T = 10 and 300K

  As it can be shown through Poisson-Schrodinger simulations [ref Oregan I ref list], the effect of L valley carrier population becomes opnly visible below 100K, giving rise to specific structure in Id-Vg and Y-Vg characteristics (Fig. 2).

Fig. 5 .Fig. 6 .

 56 Fig. 5. Extracted subthreshold swing at Vd = 30 mV versus temperature for different lengths from 10 to 300 nm

Fig. 8 .Fig. 9 .

 89 Fig. 8. Extracted low-field mobility versus channel length for linear and saturation regions along with theoretical ballistic limit

ACKNOWLEDGMENT

This work was supported by the European Union H2020 program SEQUENCE (Grant -871764).


Maroun Bsaibes

Yves Quiquempois

Stephane Plus

Adrien Masselot

Guillaume Labroille

Marianne Bigot

Jean-Baptiste Trinel

Pierre Sillard

Laurent Bigot

Light Scattering Mechanisms in Few-Mode Fibers

teaching and research institutions in France or abroad, or from public or private research centers.


Pratik Mullick

Cécile Appert-Rolland

William H Warren

Julien Pettré

Methods of density estimation for pedestrians moving in small groups without a spatial boundary

Keywords: Crowd motion, Density estimation, fundamental diagram

For a group of pedestrians without any spatial boundaries, the methods of density estimation is a wide area of research. Besides, there is a specific difficulty when the density along one given pedestrian trajectory is needed in order to plot an 'individual-based' fundamental diagram. We illustrate why several methods become ill-defined in this case. We then turn to the widely used Voronoi-cell based density estimate. We show that for a typical situation of crossing flows of pedestrians, Voronoi method has to be adapted to the small sample size. We conclude with general remarks about the meaning of density measurements in such context.

Introduction

From the point of view of crowd management, understanding the fundamental relation between crowd density and speed, i.e. fundamental diagram (FD), is extremely useful. Also for simulation based findings, FD acts as an essential tool to evaluate the capacity of simulations to predict realistic pedestrian flow. Thus, to construct a realistic FD, an effective method of density estimation is required.

In the context of self-organizing behaviour of human crowd motion, methods of density estimation are an important topic of research. Existing literature consists of a number of methods, however depending on the crowd situation a debate for the 'best' method has not yet been resolved. Most of the research have been focused on situations where the moving crowd is constrained within a physical boundary. For such cases, density estimation using Voronoi tesselation [START_REF] Duives | Quantification of the level of crowdedness for pedestrian movements[END_REF][START_REF] Steffen | Methods for measuring pedestrian density, flow, speed and direction with minimal scatter[END_REF][START_REF] Nicolas | Mechanical response of dense pedestrian crowds to the crossing of intruders[END_REF][START_REF] Zhang | Transitions in pedestrian fundamental diagrams of straight corridors and t-junctions[END_REF][START_REF] Cao | Fundamental diagrams for multidirectional pedestrian flows[END_REF] and a grid-based measure called XT method [START_REF] Edie | Discussion of traffic stream measurements and definitions[END_REF][START_REF] Duives | Quantification of the level of crowdedness for pedestrian movements[END_REF][START_REF] Bode | The emergence of macroscopic interactions between intersecting pedestrian streams[END_REF] have been reported to work well. However, there lacks a well defined method of density estimation for groups in an unbounded space.

A further difficulty arises when the fundamental diagram, instead of being plotted at a global or meso scale [START_REF] Seyfried | The fundamental diagram of pedestrian movement revisited[END_REF], is related to individual quantities [START_REF] Jelić | Properties of pedestrians walking in line -fundamental diagrams[END_REF]. In one dimension, a proxy for the density can be the inverse of the spacial headway with the predecessor. In 2 dimensions, it is more difficult to find a good estimate of the density along the trajectory of a pedestrian.

For our research, we consider crossing flows of pedestrians without any spatial constraints, where two groups of people have to go across a sporthall from predefined initial positions such that they cross each other at specified values of the crossing angle α. In the context of crossing flows, most of the previous investigations have been concentrated on α = 90 • & 180 • . We performed experimental trials for a wide range of values of α with the goal to study crossing angle dependent properties [START_REF] Mullick | Analysis of emergent patterns in crossing flows of pedestrians reveals an invariant of 'stripe' formation in human data[END_REF]. Currently, we want to investigate how the density-speed relationship for crossing flows is affected by the value of the crossing angle.

As an intermediate result, we remind here why several methods appropriate to determine a density field become inappropriate when it is the density along a pedestrian trajectory which is needed. We discuss also how Voronoi method must be adapted for small groups without boundaries. We shall end with a general discussion on the meaning of defining a density at such scales.

Experimental details

The data for crossing flows of pedestrians used for this research were obtained by experiments [START_REF] Appert-Rolland | Experimental study of collective pedestrian dynamics[END_REF][START_REF] Mullick | Analysis of emergent patterns in crossing flows of pedestrians reveals an invariant of 'stripe' formation in human data[END_REF] using live participants on the campus of University of Rennes, France. Two different sets of volunteered participants (36 on Day 1, 38 on Day 2) were roughly divided into two groups (18 or 19 per group) and were instructed to reach the other side of a sport hall. Initial positions were prescribed so that the groups had to cross each other with seven different crossing angles (0 • to 180 • , at 30 • intervals). During each trial we recorded the head trajectory of each pedestrian as a time series at a frequency of 120 Hz using a motion capture system based on infrared cameras (VICON). In Fig. 1 we show the traces of all the pedestrians for a typical trial using filtered trajectories.

In a first analysis [START_REF] Mullick | Analysis of emergent patterns in crossing flows of pedestrians reveals an invariant of 'stripe' formation in human data[END_REF], we proposed two approaches to detect the stripes formed spontaneously at the crossing. The first one uses the whole dynamics of the crossing and determines individual stripes through an algorithm tracking how individual edges between pedestrian pairs can be cut -or not -during the crossing. The second requires a single snapshot at a given time, that is matched with a regular stripe pattern.

The analysis of the experimental data based on these two stripe detection methods has allowed us to show, first, that the average number of stripes decreases when the crossing angle increases. Besides, we found a squeezing effect, meaning that during the crossing, stripes become thinner and thinner until their width reaches a minimum and expends again when pedestrian are close to the exit of the crossing region.

For further analysis, we would like to compute the fundamental diagram as measured during the crossing manoeuvre. However, the fact that groups move in an unbounded space and are relatively small raises specific difficulties as far as density estimation is concerned.

Estimating pedestrian density ρ along a pedestrian trajectory

A large number of methods have been developed to measure the density field of pedestrian flows. In order to determine the fundamental diagram at an individual scale, one has to associate a density to a specific pedestrian location. We shall see in this section that, while some methods do it without any difficulty, others completely fail once they are used along a trajectory. We illustrate how this failure occurs using a classical grid-based method. The classical grid-based method to estimate the density of pedestrians follows an Eulerian approach. We divide the entire tracking region into square-shaped grids. For d g to be the length of these grids, the classical density ρ g is given by ρ g = n/d g 2 , where n is the number of pedestrians that are located inside the grid. Note that sometimes, a single cell is used to measure density in a region of interest, for example in the crossing area of crossing flows [START_REF] Guo | A microscopic pedestrian-simulation model and its application to intersecting flows[END_REF][START_REF] Zanlungo | Macroscopic and microscopic dynamics of a pedestrian cross-flow: Part I, experimental analysis[END_REF][START_REF] Zanlungo | Macroscopic and microscopic dynamics of a pedestrian cross-flow: Part II, modelling[END_REF].

As we are interested in small systems of pedestrians, d g must be small enough to capture the immediate vicinity of a given pedestrian. In the limit of small cells, at most one pedestrian can be contained in the cell. The corresponding density is

ρ max ≡ 1/d g 2 .
If we plot the whole density field, we get typically a figure as in Fig. 2(a) with some empty cells and some cells with density ρ max . The density field is not smooth, but using a spatial average, we get a reasonable density estimate -provided the area on which we average is homogeneous. For stationary flows, it is also possible to average in time even for a rather small measurement area [START_REF] Zhang | Transitions in pedestrian fundamental diagrams of straight corridors and t-junctions[END_REF] and to get average density values with rather low fluctuations.

If instead, we decide to store the density along the trajectory of a given pedestrian, the density in the cell where this pedestrian is located will always be ρ max -a value which can be arbitrarily high. Any kind of averaging will always give ρ max (red line in Fig. 2(b)), very far from the real density average value found previously (black line in Fig. 2(b)). Actually density values are completely biased by the fact that we have selected cells conditionally, along the pedestrian trajectory. Another variant could be to measure the density along the trajectory of a pedestrian without counting the pedestrian himself. But then again there is a bias: for small cells, the density along the trajectory would always be zero, far from the real density value.

This grid-based method could be used to estimate the local density around the pedestrian trajectory only if the cell size could be chosen large compared to the pedestrian interdistance, but small compared to density gradients. This is obviously not possible for the small pedestrian groups involved in our experiment.

A similar failure for density measurements along trajectories can be observed for other methods. The systematic bias and the fluctuations around this biased value have been studied in detail for 1 dimensional systems in [START_REF] Tordeux | Quantitative comparison of estimations for the density within pedestrian streams[END_REF]. Here we illustrate it for two other methods in bidimensional space:

-The XT method originally proposed by Edie [START_REF] Edie | Discussion of traffic stream measurements and definitions[END_REF] and modified by Duives et al [START_REF] Duives | Quantification of the level of crowdedness for pedestrian movements[END_REF]: The density ρ xt (c, t) estimated in a cell c of length d x and centered on pedestrian p at time t is given by

ρ xt (c, t) = q (T q end -T q begin ) d x 2 × T , (1)
where the summation is performed over all the pedestrians q placed inside cell c at time t. T q begin and T q end denote respectively the entrance and exit time in the cell of pedestrian q with values bounded by t -0.5T and t + 0.5T .

-The Kernel method: Each pedestrian contributes to the density field through a kernel K h such that

ρ h k (X) = N i=1 K h (X -X i ). (2)
An often used kernel is the Gaussian distribution function K h = 1 2πh e -1 2h 2 X T X . Along a pedestrian trajectory, the density value which is measured at the pedestrian location directly depends on h through the normalization: ρ max = 1 2πh and can thus be arbitrarily high for small h.

In these methods, the parameters d g , d x , or h determine the typical spatial length on which the presence of a pedestrian has an effect on the density when measuring the density field. It also determines the density value at the pedestrian location when density is measured along a trajectory.

In Fig. 3 we show how the density along the trajectory of a given pedestrian measured by the 3 aforementioned methods varies as the spatial parameters d g , d x , or h decrease. In the three cases, the effect is twofold: fluctuations increase, and the average value increases towards arbitrarily high values. While the increase of fluctuations could possibly be accounted for by more averaging, the bias introduced in the average value is more problematic and makes these methods ill-defined for a measurement along a pedestrian trajectory.

Voronoi method

It is still possible to find methods that do not have this bias along pedestrian trajectories. These are methods which, by definition, scan the surroundings of a pedestrian on a scale such that the nearest neighbors' positions can be taken into account.

Among the various methods having this property, we shall here focus on Voronoi method -a method which is very popular due precisely to its good behavior in a large range of situations [START_REF] Duives | Quantification of the level of crowdedness for pedestrian movements[END_REF][START_REF] Steffen | Methods for measuring pedestrian density, flow, speed and direction with minimal scatter[END_REF][START_REF] Nicolas | Mechanical response of dense pedestrian crowds to the crossing of intruders[END_REF][START_REF] Cao | Fundamental diagrams for multidirectional pedestrian flows[END_REF][START_REF] Zhang | Transitions in pedestrian fundamental diagrams of straight corridors and t-junctions[END_REF], in spite of the fact that some fluctuations are introduced by the piecewise constant nature of the density field. We shall now describe how this method must be adapted to account for small system sizes.

The Voronoi cell of a pedestrian is the area A of the surface within which all the points are closer to that pedestrian than to any other one. The density estimate is then

ρ v = 1/A.
If we directly use this definition to construct the Voronoi diagram of the pedestrians at every instant of the trial, we notice that pedestrians on the edge of the groups may have a large and possibly infinite Voronoi cell. While this indicated that these pedestrians are not constrained by neighbors at least on one side, it may be more appropriate to associate with them the density felt on the group side.

Note that if the group was confined in a space limited by walls, the latter could be used to bound the Voronoi cells. When there is no wall, a possibility is to bound Voronoi cells with an arbitrary limit. For example a restriction to 2 m 2 was used in [START_REF] Steffen | Methods for measuring pedestrian density, flow, speed and direction with minimal scatter[END_REF]. Here we prefer a method that concentrates on the density felt on the group side.

To do this, a first correction already proposed in [START_REF] Nicolas | Mechanical response of dense pedestrian crowds to the crossing of intruders[END_REF] is to restrict the Voronoi cells to the angular sectors in which the Voronoi cell lies inside the convex hull encompassing the whole set of pedestrians. Then a correction in the density accounting for the suppressed angular sectors 2π -α must be performed as

ρv = α 2π 1 A (3)
where α is the angular sector on which the Voronoi cell is defined. This modification takes care of the angular adjustments for the agents located 'on' the convex hull and the agents whose Voronoi cell extend beyond the convex hull.

For very small samples as the ones in our experiments, specific cases have to be taken into account. Indeed a cell can extend on several sides of the group, requiring to suppress multiple sectors. Some examples are illustrated in Fig. 4(a).

In Figure 4(b) the variation of ρ v (considering Voronoi cells clipped within the convex hull but without angular modifications) and ρv (with angular modifications) as a function of time has been plotted for a typical pedestrian. Clearly ρ v produces an overestimation, which might result in an inaccurate fundamental diagram. Hence the improvement in density estimation is brought about by the angular modification in ρv .

Discussion and conclusion

We have illustrated that not all methods are appropriate for density measurements along a pedestrian trajectory, and that Voronoi method, though well defined, has to be adapted in the case of small pedestrian groups.

Several other methods can also allow to measure density along a trajectory, for example using the harmonic sum of the inter-pedestrian distances from the pedestrians that are within the field of view (FOV) [START_REF] Duives | Quantification of the level of crowdedness for pedestrian movements[END_REF], or considering the region where the convex hulls of the two crossing groups have an intersection, etc.

More generally, one must keep in mind that at small scales, a density cannot be defined in a strict way as for example for a fluid. Indeed, a proper definition of density requires that there is a scale separation between, on the one hand, the individual scale and the dimension of the area in which density is computed by a proper averaging, and, on the other hand between the latter and the scale on which density gradients occur, or on which the boundary effects become important. This is rarely the case for pedestrians, and even less for the small samples that we consider here.

What we call density is thus rather an observable that we think relevant for the navigation of pedestrians, and that we could call perceived density, following [START_REF] Steffen | Methods for measuring pedestrian density, flow, speed and direction with minimal scatter[END_REF]. It certainly depends on neighboring pedestrians, but could depend more specifically on their relative positions (rear or front), on possible visual occlusions [START_REF] Dachner | The visual coupling between neighbours explains local interactions underlying human 'flocking[END_REF], on interpersonal distances [START_REF] Geoerg | How people with disabilities influence crowd dynamics of pedestrian movement through bottlenecks[END_REF], etc. One must thus keep in mind that the choice of the perceived density plotted in the fundamental diagram actually is a hypothesis on the driving cue determining pedestrian dynamics. We shall explore various choices for this observable in our future research.

Fig. 1 .

 1 Fig.1. Traces of all the pedestrians involved in a typical trial whose crossing angle is expected to be 60 • . The dots represent the pedestrians and the tails behind each dot are the distances travelled by the pedestrian in previous 1.25 sec.

Fig. 2 .

 2 Fig. 2. Sketch representing the density variations in a grid-based method (a) when the whole density field is plotted, (b) for a time evolution in a fixed cell (black triangles) or along a pedestrian trajectory (red circles). Solid lines indicate the average values.

Fig. 3 .

 3 Fig. 3. Temporal variations of density as a function of time, along the trajectory of a given pedestrian. Density is determined by (a) the classical grid-based method, (b) the XT-method, (c) the Gaussian kernel method. Time sequences are shown for several values of the spatial scale dg, dx, or h and in (b), of the time window T .

Fig. 4 .

 4 Fig. 4. (a) Modified Voronoi cells for an instant from a typical trial of α = 90 • . The solid black line denotes the current convex hull. The dots and triangles indicate the pedestrians from two groups, which move along the x-axis and y-axis respectively. The color bar indicates the density ρv experienced by each pedestrian. (b) Time sequences of the density ρv (considering Voronoi cells clipped within the convex hull but without angular modifications) and ρv (considering modified Voronoi cells clipped within the convex hull) for a typical pedestrian, denoted by a bigger black circle in (a).


Alex Chedid

Esben Yves Quiquempois

Ravn Andresen

Laurent Bigot

Experimental characterization of mode content in fibers using S2 imaging without a DFT

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.