RED: floating-gate nMOS transistor

This example shows the distribution of the electric potential V in a nanoscale floating-gate nMOS (Metal-Oxide-Semiconductor) transistor working in inversion conditions. The prefix "n" indicates that the electric current is due to negatively charged electrons. This simulation involves a non-homogeneous problem with internal interfaces and it has been solved with the HDG method applied to the Darcy problem. A more detailed description of this test case is illustrated in [HDG2018].

1. Running the case

The command line to run this case is

mpirun -np 4 feelpp_toolbox_mixed-poisson-model_3DP1 --case "github:{repo:toolbox,path:examples/modules/poisson/examples/red}"
sh

2. Data files

The case data files are available in Github here

3. Geometry

A scheme of a realistic floating-gate nMOS transistor used as a nonvolatile memory device is shown with the schematic representation of the two-dimensional cross-section in which the control gate is now included for sake of simplicity in the following Figure:

nMOS

The device is composed of a pair of n-doped source and drain regions, a p-doped substrate region and a silicon dioxide (SiO2) region in which a control gate and a floating gate are buried. The floating gate is isolated from the control gate by an inter-polysilicon dielectric (IPD) layer beneath the control gate.

The geometry has been generated with GMSH:

h = 0.1;

Lx = DefineNumber[480*10^(-9), Name "x Lenght"];
Ly = DefineNumber[320*10^(-9), Name "y Lenght"];
Lz = DefineNumber[320*10^(-9), Name "z Lenght"];
Lj = DefineNumber[160*10^(-9), Name "j Lenght"];
t0x = DefineNumber[10*10^(-9), Name "Gate Length"];

// Discretization
npt = DefineNumber[20, Name "Nb point interp"];
l = DefineNumber[3.5, Name "Interv for tanh"];

A = newp; Point(A) = {0,0,0,h};
B = newp; Point(B) = {Lx,0,0,h};
C = newp; Point(C) = {Lx,Ly,0,h};
D = newp; Point(D) = {2./3.*Lx,Ly,0,h/3};
E = newp; Point(E) = {2./3.*Lx-Lx/8.,Ly,0,h/3};
F = newp; Point(F) = {1./3.*Lx+Lx/8.,Ly,0,h/3};
G = newp; Point(G) = {1./3.*Lx,Ly,0,h/3};
H = newp; Point(H) = {0,Ly,0,h};
K = newp; Point(K) = {2./3.*Lx-Lx/8.,Ly+t0x,0,h/3};
L = newp; Point(L) = {1./3.*Lx+Lx/8.,Ly+t0x,0,h/3};
M = newp; Point(M) = {0,Ly-Lj,0,h};
N = newp; Point(N) = {1./3.*Lx,Ly-Lj,0,h};
P = newp; Point(P) = {Lx,Ly-Lj,0,h};
Q = newp; Point(Q) = {2./3.*Lx,Ly-Lj,0,h};

AB = newl; Line(AB) = {A,B};
BP = newl; Line(BP) = {B,P};
PC = newl; Line(PC) = {P,C};
CD = newl; Line(CD) = {C,D};
DE = newl; Line(DE) = {D,E};
EF = newl; Line(EF) = {E,F};
FG = newl; Line(FG) = {F,G};
KL = newl; Line(KL) = {K,L};
GH = newl; Line(GH) = {G,H};
HM = newl; Line(HM) = {H,M};
MA = newl; Line(MA) = {M,A};
MN = newl; Line(MN) = {M,N};
PQ = newl; Line(PQ) = {P,Q};

pListNF[0] = N;
pListQE[0] = Q;
pListGL[0] = G;
pListDK[0] = D;
For i In {1:npt-1}
    xhi = i*Lx/8./npt;
    eta = 64*Lj/(Lx*Lx)*xhi*xhi;
    xNF = 1./3.*Lx + xhi;
    xQE = 2./3.*Lx - xhi;
    y = eta + Ly - Lj;
    pListNF[i] = newp;
    Point(pListNF[i]) = {xNF,y,0,h};
    pListQE[i] = newp;
    Point(pListQE[i]) = {xQE,y,0,h};
    eta = t0x/2.*(Tanh(16.*l/Lx*xhi-l)+1);
    xGL = xNF;
    xDK = xQE;
    y = Ly + eta;
    pListGL[i] = newp;
    Point(pListGL[i]) = {xGL,y,0,h/3};
    pListDK[i] = newp;
    Point(pListDK[i]) = {xDK,y,0,h/3};
EndFor
pListNF[npt] = F;
pListQE[npt] = E;
pListGL[npt] = L;
pListDK[npt] = K;

NF = newl; Spline(NF) = pListNF[];
QE = newl; Spline(QE) = pListQE[];
GL = newl; Spline(GL) = pListGL[];
DK = newl; Spline(DK) = pListDK[];

Bll = newll; Line Loop(Bll) = {AB,BP,PQ,QE,EF,-NF,-MN,MA};
OB = news; Plane Surface(OB) = {Bll};

Dll = newll; Line Loop(Dll) = {PC,CD,DE,-QE,-PQ};
OD = news; Plane Surface(OD) = {Dll};

Sll = newll; Line Loop(Sll) = {FG,GH,HM,MN,NF};
OS = news; Plane Surface(OS) = {Sll};

Gll = newll; Line Loop(Gll) = {DE,EF,FG,GL,-KL,-DK};
OG = news; Plane Surface(OG) = {Gll};

// 3D
out[] = Extrude{0,0,Lz} {Surface{OB,OD,OS,OG};};

// For i In {0:#out[]-1}
//     Printf("out[%g] : %g",i, out[i]);
// EndFor

Physical Surface("Bulk") = {out[2]};
Physical Surface("Drain") = {out[13]};
Physical Surface("Source") = {out[20]};
Physical Surface("Gate") = {out[30]};
Physical Surface("IntBD") = {out[4],out[5]};
Physical Surface("IntBG") = {out[6]};
Physical Surface("IntBS") = {out[7],out[8]};
Physical Surface("IntDG") = {out[14]};
Physical Surface("IntSG") = {out[19]};
Physical Surface("WallB") = {OB,out[0],out[3],out[9]};
Physical Surface("WallD") = {OD,out[10],out[12]};
Physical Surface("WallS") = {OS,out[17],out[21]};
Physical Surface("WallG") = {OG,out[24],out[29],out[31]};

Physical Volume("OmegaB") = {out[1]};
Physical Volume("OmegaD") = {out[11]};
Physical Volume("OmegaS") = {out[18]};
Physical Volume("OmegaG") = {out[25]};
cpp

N.B.: the unit of measure used in the .geo file are in mm.

4. Input parameters

Physical parameters:

Name Symbol Value Unit

electron charge

q

1.602e-19

C

permittivity of vacuum

ε0

8.854e-12

Fm1

relative permittivity of silicon

εsir

11.7

stem:[ - ]

relative permittivity of silicon dioxide

εoxr

3.9

stem:[ - ]

thermal voltage (at 300K)

Vth

0.02589

V

intrinsic concentration (at 300K)

ni

1.45e16

m3

Device parameters:

Name Symbol Value Unit

horizontal length

Lx

480e-9

m

vertical length

Ly

320e-9

m

width

Lz

320e-9

m

oxide thickness

tox

10e-9

m

source and drain lengths

|ΓS|,|ΓD|

160e-9

m

channel length

Lch

40e-9

m

junction depth

tJ

106e-9

m

source potential

ˉVS

0

V

bulk potential

ˉVB

0

V

drain potential

ˉVD

0

V

source doping concentration

ˉNS

1e26

m3

bulk doping concentration

ˉNB

5e22

m3

drain doping concentration

ˉND

1e26

m3

5. Model & Toolbox

  • The main goal of the simulations is to determine, via the integral boundary condition, the value attained on ΓG at inversion conditions by the electric potential. This value is the threshold voltage of the device, denoted henceforth by VT , and is a fundamental design parameter in integrated circuit nanoelectronics. The accurate prediction of the threshold voltage may be significantly influenced by the numerical treatment of the heterogeneous material device structure across the interface Γint. With this respect the HDG framework appears to be the ideal computational approach because it allows to enforce, at the same time and in a strong manner, (a) the jump of the normal component of the electric displacement vector across the interface and (b) the continuity of the electric potential, at suitable quadrature nodes on each face belonging to Γint, unlike standard displacement-based finite element discretization schemes which satisfy property (b) but fail at satisfying property (a).

  • Config file:

directory=toolboxes/hdg/mixed-poisson/red

[picard]
itol=1e-15
itmax=5

[exporter]
element-spaces=P0

[gmsh]
filename=$cfgdir/new_red.geo
hsize=5e-5

[mixedpoisson]
// pc-type=gasm
// sub-pc-factor-mat-solver-package-type=umfpack
// sub-pc-type=lu
ksp-monitor=1
ksp-rtol=1e-14
conductivity_json=k
filename=$cfgdir/new_red.json

[bdf]
order=1
[ts]
time-initial=0.0
time-step=10
time-final=150
steady=1 #false #true
ini

5.1. Materials

    "Materials":
    {
        "OmegaG":
        {
            "name":"oxide",
	    	"k":"3.9 *  8.854*10^(-15) * 10^15", // k = epsR * eps0 / q
            "scale_flux":"1e-15 / (3.9 * 8.854^(-15))" // E = D/(epsR*eps0)
        },
        "OmegaS":
        {
            "name":"silicium",
	    	"k":"11.7 *  8.854*10^(-15) * 10^15", // k = epsR * eps0 / q
            "scale_flux":"1e-15 / (11.7 * 8.854^(-15))" // E = D/(epsR*eps0)
        },
        "OmegaD":
        {
            "name":"silicium",
	    	"k":"11.7 *  8.854*10^(-15) * 10^15", // k = epsR * eps0 / q
            "scale_flux":"1e-15 / (11.7 * 8.854^(-15))" // E = D/(epsR*eps0)
        },
        "OmegaB":
        {
            "name":"silicium",
	    	"k":"11.7 *  8.854*10^(-15) * 10^15", // k = epsR * eps0 / q
            "scale_flux":"1e-15 / (11.7 * 8.854^(-15))" // E = D/(epsR*eps0)
        }
    },
json

5.2. Boundary conditions

Boundary conditions
V=ˉVj+ˉV|bi,j on Γj,j=S,D,BD_n_=0 on ΓlatsiΓlatoxΓGD_n_dΣ=qˉNBδLchLz on ΓG
Interface conditions
[V]Γint=0[D_]Γint=qˉNBδ

The quantities ˉV|bi,j are the built-in potentials associated with the subdomain regions Ωj,j=S,D,B, while ˉNB is the concentration of ionized dopants in the bulk region and δ is the width of the accumulation region in the y direction of the channel.

    "BoundaryConditions":
    {
        "potential":
        {
            "InitialSolution":
            {
                "OmegaG":
                {
                    "expr":"0.0"
                },
                "OmegaS":
                {
                    "expr":"0.0"
                },
                "OmegaB":
                {
                    "expr":"0.0"
                },
                "OmegaD":
                {
                    "expr":"0.0"
                }
            },
            "SourceTerm":
            {
                "OmegaG":
                {
                    "expr":"0.0:x:y:z"
                },
                "OmegaB":
                {
                    "expr":"0.0:x:y:z"
                },
                "OmegaD":
                {
                    "expr":"0.0:x:y:z"
                },
                "OmegaS":
                {
                    "expr":"0.0:x:y:z"
                }
            },
            "Dirichlet":
	    	{
	    		"Source":
        		{
                	"expr":"0.5865"
               	},
	    		"Drain":
        		{
                	"expr":"0.5865"
               	},
	    		"Bulk":
        		{
                	"expr":"-0.3897"
               	}
	    	},
            "Neumann":
            {
	    		"WallG":
        		{
                	"expr":"0.0"
               	},
	    		"WallS":
        		{
                	"expr":"0.0"
               	},
	    		"WallD":
        		{
                	"expr":"0.0"
               	},
	    		"WallB":
        		{
                	"expr":"0.0"
               	}
            }
		},
		"flux":
		{
	    	"Integral":
	    	{
        		"Gate":
                {
					"expr": "-1.1858*10^(-17) * 10^15"
                }
 	    	},
            "InterfaceCondition":
            {
                "IntSG":
                {
                    "expr":"11*10^(-9) * 10^15"
                },
                "IntDG":
                {
                    "expr":"11*10^(-9) * 10^15"
                },
                "IntBG":
                {
                    "expr":"11*10^(-9) * 10^15"
                }
            }
    	}
     },
json

6. Outputs

6.1. Potential and current

Spatial distribution of the potential and direction of the current in the electronic device at steady state in inversion condition.

potential+arrow

Spatial distribution of the electric field and its direction in the electronic device at steady state in inversion condition.

electricfield with arrow

6.2. 3D model

In the window below, you can manipulate the 3D model at steady state.

Click top left button on opengl window to change basic visualisation features

6.3. Measures

The predicted value of VT is positive and given by VcomputedT=0.854926 V.

The sign of VcomputedT is in physical agreement with the fact that the transistor is of n type, so that electron charges (negative) are attracted from the bulk region up towards the gate contact.

To quantitatively assess the accuracy of the estimate of the threshold voltage we adopt the analytical theory for an ideal MOS system developed in [Muller1986][Chapter 8] and use formula (8.3.18) of [Muller1986] to obtain VidealT=0.8591 V, which agrees with the value predicted by the numerical simulation.

7. Bibliography

References for this example
  • [] Muller RS, Kamins TI, Chan M, Ko PK. Device electronics for integrated circuits. 1986.

  • [] Bertoluzza S, Guidoboni G, Hild R, Prada D, Prud’homme C, Sacco R, Sala L, Szopos M. An implementation of HDG methods with Feel++. Application to problems with integral boundary conditions. In preparation.