Skip to content
Publicly Available Published by De Gruyter February 12, 2019

Approximate Solution of Coupled Schrödinger and Poisson Equation in Inversion Layer Problem: An Approach Based on Homotopy Perturbations

  • Tijana Kevkić and Vladica Stojanović ORCID logo EMAIL logo

Abstract

In this paper, the homotopy perturbation method (HPM) is applied to the coupled set of Schrödinger–Poisson (SP) equations in inversion layer problem for obtaining the approximate analytical solution. Inversion layer of n-type is considered, and the electric quantum limit is assumed. By introducing some dimensionless quantities, the SP system has been turned into one which can be solved along the infinite interval. After some appropriate transformations, the infinite interval has been reduced to finite one (0,1), and recurrence series of the HPM approximate solutions of the coupled SP system have been obtained. The existence and convergence of obtained HPM approximate solutions have been formally proved. Moreover, these solutions show relative simple mathematical form, as well as high degree of accuracy what is desirable for semiconductor device modelling.

1 Introduction

Application of an electric field normal to the surface of a semiconductor results in the redistribution of charge carriers and creation of a space-charge layer at the surface, causing a band bending there [1]. This leads to the carrier confinement within a narrow potential well close to the semiconductor surface and, consequently, a quantized motion in the direction perpendicular to the surface [1], [2]. The energy levels of these carriers are grouped in discrete subbands, each of which corresponds to a quantized level for motion in the direction perpendicular to the surface, with a continuum for motion in the plane parallel to the semiconductor surface [3]. In the case of inversion surface layer, most of the carriers responsible for current transport reside in the lowest subband [4], and their motion is described by Schrödinger wave equation. On the other side, the charge distribution of these carriers must satisfy the Poisson equation. Thus, Schrödinger and Poisson (SP) equations are coupled here, and these need to be solved simultaneously and self-consistently in order to find the quantized energy levels and carrier concentration present in the inversion layer.

Besides the quantum charge transport in semiconductor devices, e.g. metal–oxide–semiconductor transistors (MOSTs), the SP equations also govern many other important physical systems, such as the quantum many-body system and time-dependent density functional theory [5]. Therefore, a lot of attempts to solving the considered SP system have been made. Unfortunately, as the coupled SP system cannot be solved analytically, the solution procedures are exclusively numerical [6], [7], [8]. However, applications of the numerical methods in semiconductor device modelling are mostly limited to the single device simulation or restricted only at finite approximation intervals [9], [10], [11] without giving the analytical expressions that would give an insight into overall behaviour of the obtained solutions. On the other side, the analytical solutions are necessary for understanding the behaviour of charge carriers in surface layers and preferable from a circuit modelling point of view due to simplicity in their format and fast computational speed. In order to obtain an analytical solution of the coupled SP equations in the case of inversion layer of MOST, different approximations have been introduced [12]. Commonly used one is that the barrier height at the oxide–semiconductor interface is infinite, and as a consequence, the wave function at the interface is forced to zero [13]. The other is the triangular potential well approximation within the semiconductor, which results in a constant electric field throughout the semiconductor, equal to the surface electric field [14]. Further, under the assumption that electrons occupy only the lowest subband, the variational approach gives a good estimate of the eigen energy of this subband [15].

Unlike the previous approximate methods relating mainly to the finite intervals, here we started from approximate solving of the SP system along the whole semiconductor, i.e. along the infinite interval. In this purpose, we introduce the homotopy perturbation method (HPM), similarly as in [16], where the coupled SP system was considered in the case of the accumulation layer. However, unlike the accumulation layer, here the so-called effective potential becomes unbounded making the inverse layer problem considerably more complicated (see, for more detail, Section 3). To ensure the convergence on infinite (0,) interval, after some appropriate transformations, we have reduced this interval on the finite one (0,1). In that way, the system of dimensionless coupled equations that can be solved by using the HPM has been obtained. The HPM, proposed by He [17], [18], [19], [20], [21], found significant applications in approximate solving of the nonlinear equations of various types [22], [23], [24], mainly in the physical sciences [25], [26], [27], [28], [29]. Here, starting from the coupled system of homotopy equations, a procedure that provides convergence with exact solutions of the coupled SP equations is developed. That convergence is formally confirmed, and this is one of the most important advantages of HPM technique in comparison to some previously used approximate methods. Furthermore, the obtained HPM approximate solutions can be computed recursively with an arbitrary numerical precision, and the whole approximation procedure can be implemented by using the appropriate software.

2 Self-Consistent Field SP Equations

The electrons in an inversion layer of NMOS transistor are confined in potential well whose width in the ξ direction, the direction perpendicular to the semiconductor surface, is small compared to the wavelengths of the electrons. The energy of these electrons is quantized, and the semiconductor conduction band is split into discrete subbands. In the effective mass approximation, under the assumption that only lowest subband E0 is occupied, the electronic wave function ψ0(ξ) for the lowest subband is determined by the Schrödinger equation:

(1)22md2ψ0(ξ)dξ2eϕe(ξ)ψ0(ξ)=E0ψ0(ξ).

Here, m is the effective mass of electrons, e is the electron charge, ℏ is the reduced Planck constant, and ϕe(ξ) is electrostatic potential in the semiconductor determined by Poisson equation:

d2ϕe(ξ)dξ2=ρεsi.

Here, ρ is space charge density, and εsi is permittivity of semiconductor. In the so-called depletion layer approximation, ϕe(ξ) is divided into two parts:

ϕe(ξ)=ϕdep(ξ)+ϕI(ξ),

where ϕdep(ξ) is only determined by the space charge density due to the ionized acceptors in the depletion layer, and ϕI(ξ) is determined by the electron concentration in the inversion layer. Because the inversion layer is found to be very thin compared with the depletion layer thickness, for most practical cases it is possible to approximate ϕdep(ξ) by the linear term giving ϕdep(ξ)=Kξ, where K is the electric field at the surface due to ionized acceptors charge in the depletion layer.

By putting this expression in (1), we get reduced Schrödinger equation for ψ0(ξ):

[22md2dξ2+eKξeϕI(ξ)]ψ0(ξ)=E0ψ0(ξ).

The distance between the conduction band edges of oxide and semiconductor presents the potential barrier at the semiconductor–oxide interface, placed at ξ = 0. For commonly used SiSiO2 system, the height of that barrier is 3.1 eV [30], what is about hundred times larger than the energy-level separation of the quantized motion perpendicular to the interface [4]. That is why it can be assumed that the potential barrier at interface is infinite, and consequently, the electron wave function vanishes there, i.e.

ψ0(ξ0)=0.

Moreover, the inversion charge is distributed near the surface of the semiconductor, what indicates that both the wave function and its derivative should vanish in the infinite distance. This results in the following boundary condition:

ψ0(ξ)=0.

At last, the wave function ψ0(ξ) must obey the normalization condition:

0|ψ0(ξ)|2dξ=1.

On the other hand, the electrostatic potential ϕI(ξ) in Schrödinger equation is a solution of Poisson equation:

(2)d2ϕI(ξ)dξ2=eN0εsi|ψ0(ξ)|2,

where N0 is the total number of electrons in inversion layer as only the lowest level E0 is occupied. In the bulk of the semiconductor (ξ), the space charge density is zero, and the electrostatic potential ϕI(ξ) stays constant. For convenience, it is chosen to be the reference potential. Thus, the potential boundary conditions at infinity is given by

ϕI(ξ)=dϕI(ξ)dξ=0,whenξ.

In the following, we can get a set of equations that are independent of the parameters e, m, ℏ, and K, by introducing the following dimensionless quantities:

z=ξ(emK2)1/3,ϕ(z)=ϕI(ξ)(em2K2)1/3,ε0=E0(me22K2)1/3,θ(z)=ψ0(ξ)(2emK)1/6

and α=(eN0)/(εsK). With these new terms, SP system (1)–(2) change into

(3)12d2θ(z)dz2+[zϕ(z)]θ(z)=ε0θ(z)
(4)d2ϕ(z)dz2=αθ2(z),

with the boundary conditions:

(5)θ(z)=0,whenz0andz,ϕ(z)=dϕ(z)dz=0,whenz,

and the normalized condition 0|θ(z)|2dz=1.

Figure 1: Left: Graph of the transformation z(x)=ln(1/x)$z(x)=\ln(1/x)$. Right, above: Graphs of the potential ϕ(z)$\phi(z)$ and effective potential g~(z)$\tilde{g}(z)$, when z∈(0,+∞)$z\in(0,+\infty)$. Right, below: Graph of the effective potential g(x)$g(x)$, when x∈(0,1)$x\in(0,1)$.
Figure 1:

Left: Graph of the transformation z(x)=ln(1/x). Right, above: Graphs of the potential ϕ(z) and effective potential g~(z), when z(0,+). Right, below: Graph of the effective potential g(x), when x(0,1).

3 Solving the Coupled SP Equations with HPM

In general, as it is already mentioned, the coupled system of (3)–(4) can be solved numerically or by using some approximate analytical methods mostly based on popular triangular potential well approximation or a variation calculus. However, first of these approximate methods fails when the inversion layer charge density is not sufficiently low, while the second is not suitable for a device model aimed at circuit simulation [4]. In addition, there is no general proof on convergence of solutions obtained using mentioned approximate analytical methods [31]. In order to obtain more exact approximate solution, which convergence can be formally confirmed, here we propose the application of the HPM for solving the coupled SP equations.

In the following, by transformation z(x):=ln(1/x) we make the change of z variable in SP system, i.e. introduce the new variable x(z):=exp(z), where x(0)=1 and x(+)=0. Accordingly, the introduced transformation is a monotonic, nonnegative function that translates the infinite z interval (0,+) to the unit x interval (0,1) (Fig. 1, left panel). In this way, as already mentioned, the simplified convergence of the approximate solutions of the SP system to the exact solution is enabled.

Firstly, after transformation z=z(x), the coupled (3)–(4) can be additionally rearranged in the following dimensionless form:

(6)x2d2f(x)dx2+xdf(x)dxf(x)g(x)=0
(7)x2d2g(x)dx2+xdg(x)dx+Af2(x)=0.

Here, we denoted A=2α, f(x)=θ(z(x)), and g(x)=2(z(x)ϕ(z(x))ε0) are effective potential. Let us notice that effective potential, observed as a function of z variable, e.g. g~(z):=2(zϕ(z)ε0), in accordance to (3)–(5) satisfies following conditions:

g~(z)2z,whenz,d2g~(z)dz2=Aθ2(z)<0,whenz>0.

From these conditions, it is apparent that g~(z) is a concave function, unbounded at the end of the (0,+) interval. On the other hand, effective potential g(x), even unlimited at x = 0, is defined at the finite unit interval and, similarly as potential ϕ(z), is a convex function on x(0,1) (Fig. 1, right panels). Among the others, this is one of the main motivation for introducing the aforementioned transformation of the z to the x variable. Finally, in accordance to the above transformation, the boundary conditions given by (5) now become

(8)f(0)=f(1)=0,andg(x)lnxxdg(x)dx2,

when x → 0. Also, for function f(x), the normalized condition 01x1|f(x)|2dx=1 holds.

In order to apply HPM technique on the coupled (6)–(7) with boundary conditions (8), we have constructed the following homotopy equations:

(9)(1p)[x22F(x,p)x2x2d2f0(x)dx2]+ph1[x22F(x,p)x2+xF(x,p)xF(x,p)G(x,p)]=0
(10)(1p)[x22G(x,p)x2x2d2g0(x)dx2]+ph2[x22G(x,p)x2+xG(x,p)x.+AF2(x,p)]=0.

Here, p[0,1] is the embedding parameter, h1,h2(0,1] are the auxiliary parameters, and f0(x), g0(x) are the so-called initial solutions. The basic assumption of the HPM is that solutions of coupled homotopy (9)–(10) can be expressed as power series in p:

(11a)F(x,p)=j=0pjfj(x),
(11b)G(x,p)=k=0pkgk(x).

Owing to these series is obvious why f0(x) and g0(x) are the initial solutions of (9)–(10), respectively. Namely, from (9)–(10), as well as from (11a)–(11b), when p = 0 follows:

F(x,0)=f0(x)andG(x,0)=g0(x).

In that way, f0(x) and g0(x) can be chosen arbitrarily, but also, they should satisfy the boundary conditions:

(12)f0(0)=f0(1)=0,andg0(x)lnxxdg0(x)dx2,

when x → 0, as well as the normalized condition:

01x1|f0(x)|2dx=1.

On the other hand, when p = 1, (9)–(10) turn into (6)–(7), with the same boundary conditions (8). According to (11a)–(11b), solutions of (6)–(7) will be

(13a)f(x)=limp1F(x,p)=j=0fj(x),
(13b)g(x)=limp1G(x,p)=k=0gk(x),

under the condition of convergence of both series in (13a)–(13b).

Now, rearranging terms in homotopy (9)–(10), and by substituting (11a)–(11b) in them, it follows:

(14){x2j=1pjd2fj(x)dx2+p[x2h1j=0pjd2fj(x)dx2x2j=1pjd2fj(x)dx2+xh1j=0pjdfj(x)dxh1(j=0pjfj(x))(k=0pkgk(x))]=0,x2k=1pkd2gk(x)dx2+p[x2h2k=0pkd2gk(x)dx2x2k=1pkd2gk(x)dx2+xh2k=0pkdgk(x)dx+Ah2(j=0pjfj(x))2]=0.

Equalization of the expressions with identical powers pj, pk in (14) leads to the following differential equations:

(15){x2d2f1(x)dx2+h1[x2d2f0(x)dx2+xdf0(x)dxf0(x)g0(x)]=0x2d2g1(x)dx2+h2[x2d2g0(x)dx2+xdg0(x)dx+Af02(x)]=0

and, for k ⩾ 2:

(16){x2d2fk(x)dx2+x2(h11)d2fk1(x)dx2+xh1dfk1(x)dxh1j=0k1fj(x)gkj1(x)=0x2d2gk(x)dx2+x2(h21)d2gk1(x)dx2+xh2dgk1(x)dx+Ah2j=0k1fj(x)fkj1(x)=0.

Furthermore, we suppose that for arbitrary k=1,2, the following boundary conditions are fulfilled:

fk(0)=fk(1)=0,andgk(x)lnxxdgk(x)dx0,

when x → 0. These conditions, along with (12), enable that functions f(x) and g(x), defined by (13a)–(13b), satisfy boundary conditions (8). It is clear that (15)–(16) can be solved recursively on fk(x) and gk(x), using the double integration. The solutions of above differential equations, for any k=0,1,2,, give the so-called HPM approximations:

(17)f^k(x):=j=0kfj(x),g^k(x):=j=0kgj(x)

as the estimates of unknown functions f(x) and g(x), respectively. Notice that HPM approximations, at the same time, satisfy the boundary conditions (8). On the other hand, for any k=0,1,2,, HPM approximations f^k(x) of the wave function f(x) must obey normalization condition:

(18)01x1|f^k(x)|2dx=1.

In accordance to this, the multiplicative factor h1 can be chosen such that the condition (18) provides, while h2 enables a faster convergence of the HPM approximations (17). Finally, we give some sufficient conditions for the convergence of these series.

Theorem 1

Let {fk(x)}k=0 and {gk(x)}k=0 be the series of continuous functions, defined on (0,1) by recurrence relations (15)–(16), such that for some C > 0 and arbitrary x(0,1) following inequalities

(19)max{|fk(x)|,|dgk(x)dx|}Cxk+1

hold for each k=0,1,2. Then, series of HPM approximations {f^k(x)}k=0 and {g^k(x)}k=0 defined by (17) converge on (0,1), and their limits are given by (13a)–(13b); i.e. they are solutions of the (6)–(7), respectively.

Proof.

According to assumption given in (19), for fixed but an arbitrary x(0,1) and k ⩾ 1, we have

|fk(x)|1/kC1/kx1+1/k.

On the other side, according to the second equations in (15)–(16), for the sequence {gk(x)}, we have

(20)|g1(x)|h2(|g0(x)|+dx|dg0(x)x|+Adxf02(x)x2dx)h2(|g0(x)|+C(1+AC)x22)

and, when k ⩾ 2:

(21)|gk(x)|(1h2)|gk1(x)|+h2dx|dgk1(x)x|+Ah2j=0k1dx|fj(x)x||fkj1(x)x|dx(1h2)|gk1(x)|+h2C(1k+AC)xk+1k+1.

Applying the Cauchy–Hadamard theorem, we find that radius of convergence of the power series in (11a), denoted as r1(x), satisfies inequalities:

r1(x)=(lim supk|fk(x)|1/k)1limk(C1/kx(1+1/k))=1x>1.

On the other hand, for the radius of convergence r2(x) of the power series in (11b), according to (21), is valid:

(22)r2(x)=limk|gk1(x)gk(x)|limk|gk1(x)|(1h2)|gk1(x)|+h2C(1k+AC)xk+1k+1=11h2>1.

Thus, both radius of convergence of these power series are greater than 1.

Further on, applying the induction method in (20)–(21), for an arbitrary k ⩾ 1 follows:

|gk(x)|h2(1h2)k1|g0(x)|+h2C(1+AC)j=0k1(1h2)jxkj+1kj+1h2(1h2)k1|g0(x)|+h2C(1+AC)xk+1j=0k1(1h2x)jh2(1h2)k1|g0(x)|+h2C(1+AC)xk+2x1+h2.

Now, according to the definition of HPM approximations f^k(x) and g^k(x), given by (17), for an arbitrary x(0,1) and k ⩾ 1, we have

|f^k(x)|j=0k|fj(x)|Cj=0kxj+1,|g^k(x)||g0(x)|+h2|g0(x)|j=0k1(1h2)j+h2C(1+AC)x1+h2j=0kxj+2.

Finally, when k → ∞, the inequalities above imply

|j=0fj(x)|=limk|f^k(x)|Cj=0xj+1=Cx1x<+,|j=0gj(x)|=limk|g^k(x)||g0(x)|[1+h2j=0(1h2)j]+h2C(1+AC)x1+h2j=0xj+22|g0(x)|+h2C(1+AC)x2(1x)(x1+h2)<+.

In that way, the convergence of the sequences {f^k(x)}k=0 and {g^k(x)}k=0 is proved; i.e. both series converge on (0,1). According to Abel’s theorem, it follows that functions F(x,p) and G(x,p), defined by (11a)–(11b), are continuous from the left at p = 1. Therefore, (13a)–(13b) hold, i.e. the series j=0fj(x) and j=0gj(x) are solutions of (6)–(7), respectively.

Remark 1

Thanks to the appropriate choices of the electronic wave function and the electrostatic potential, the conditions of the previous theorem can be easily fulfilled in solving of the inversion layer problem. In addition, (22) indicates to the fact that “the most preferable” value of the second auxiliary parameter is h2 = 1. In that case, the power series in (11b) has a maximum radius of convergence:

r2(x)=(1h2)1=+.

The practical interpretation (and validation) of these conditions will be seen immediately, further in this article.

4 Validation of the HPM Procedure

The practical application and the quality of approximations of the aforementioned HPM procedure have been investigated and described in more detail in this section. The primary question in HPM is choice of suitable initial solutions that will ensure (fast) convergent series in (17). Additionally, these initial solutions must satisfy the boundary conditions (8) and should be as simple as possible for integration.

Having in mind the form of solution of the Schrödinger equation (3) in the electric quantum limit [6], here for electron wave function of the lowest subband we propose the following initial approximation:

θ0(z)=P(z)exp(az),a>0.

Here, P(z) is a polynomial on z ⩾ 0, which satisfies condition P(0)=0, and as such, in addition, the normalized condition:

0|θ0(z)|2dz=1

is fulfilled. According to the transformation z(x)=ln(1/x), as well as the aforementioned facts, we obtain as the initial solution of the homotopy (9) the following function:

f0(x):=θ0(z(x))=xaP(lnx1).

Similarly as in [4], we have taken P(z)=2z and a = 1; i.e. the first sequence of HPM approximations is based on the function f0(x)=2xlnx, where x(0,1). Note that, in this case is

01x1|f0(x)|2dx=401xln2xdx=1,

i.e. the normalized condition obviously holds.

On the other side, for the initial approximation ϕ0(z) of the potential ϕ(z), in general, any function that satisfies the second boundary condition in (5) can be used; i.e. such that

ϕ0(z)0,dϕ0(z)/dz0,whenz.

With regard to this, as the initial approximation ϕ0(z) of the potential ϕ(z), we have used two specific functions:

  1. Firstly, we have taken the trivial function ϕ0(z)0; i.e. the initial solution of the homotopy (10) will be

    g0(x):=2(lnx1ϕ0(lnx1)ε0)=2(lnx+ε0).
  2. The second initial approximation for the potential ϕ(z) has been obtained by using Poisson equation (4), i.e. by double integration:

    ϕ0(z)=αdzθ02(z)dz=αe2z(z2+2z+32).

    Then, as in the previous case, using transformation z(x)=lnx1, initial solution of (10) becomes

    g0(x)=2(αx2(ln2x2lnx+32)+lnx+ε0).

Therefore, by using the recurrence relations in (15)–(16), the functions {fk(x)} and {gk(x)} can be obtained for an arbitrary k ⩾ 1. Hereinafter, these two series will be labelled as A and B, respectively, and their first few components are given as below:

Series A:

f1(x)=43xln3x+3xln2x2xlnx,g1(x)=x2(2ln2x6lnx+7);f2(x)=281(18x3ln3x90x3ln2x+171x3lnx+108xln4x783xln3x+2835xln2x6156xlnx),g2(x)=x26(8ln4x50ln3x+138ln2x201lnx+132);f3(x)=22187(648x3ln5x5994x3ln4x+24156x3ln3x51876x3ln2x+57033x3lnx5832xln5x+83106xln4x574452xln3x+2322594xln2x4846392xlnx),g3(x)=x2288(189x2+32x2ln4x224x2ln3x+616x2ln2x532x2lnx16200ln4x+75264ln3x218988128ln6x+2112ln5x221256ln2x+329616lnx);

etc.

Series B:

f1(x)=19(6x3ln3x27x3ln2x+48x3lnx+12xln3x+27xln2x18xlnx),g1(x)=x2(2ln2x6lnx+7);f2(x)=111390625(607500x5ln5x+303750x5ln4x10388250x5ln3x+29279475x5ln2x34497900x5lnx+3375000x3ln5x23906250x3ln4x+89156250x3ln3x152203125x3ln2x+119031250x3lnx+30375000xln4x220218750xln3x+797343750xln2x1809843750xlnx),g2(x)=x2576(653x296x2ln4x240x2ln3x+1200x2ln2x1606x2lnx14656768ln4x+4800ln3x13248ln2x+21280lnx);

etc.

Notice that the both series of HPM approximations meet conditions of Theorem 1. For both series have been used following values of parameters α = 0.5 and ε0 = 3. On the other hand, auxiliary parameter h1 is chosen such that the normalized condition (18) is satisfied, while for the other parameter the constant value h2 = 1 is taken. According to these, and by using (17), the appropriate HPM approximations f^k(x) and g^k(x) can be easily obtained. The whole procedure of computation fk(x) and gk(x), i.e. the appropriate HPM approximations f^k(x) and g^k(x), when k=1,2,,10, has been implemented in the software package MATHEMATICA 11.0. In the same procedure, three different measures of convergence of the obtained HPM approximations have been computed, i.e. corresponding approximation errors:

  1. The maximum iteration differences (MID):

    MID(f^k):=f^kf^k1:=max0<x<1|fk(x)fk1(x)|,MID(g^k):=g^kg^k1:=max0<x<1|gk(x)gk1(x)|,

    represent the measures of convergence of the HPM approximations in relation to their approximation’s order k ⩾ 1.

  2. The maximum absolute errors (MAE):

    MAE(f^k):=max0<x<1|x2d2f^k(x)dx2+xdf^k(x)dxf^k(x)g^k(x)|,MAE(g^k):=max0<x<1|x2d2g^k(x)dx2+xdg^k(x)dx+Af^k2(x)|,

    are quantitative indicators of the accuracy of HPM approximations errors in each iteration k ⩾ 1. Note that for exact solutions f(x), g(x) of the SP system, given by (6)–(7), will MAE(f)=MAE(g)=0. As the values of f(x), g(x) cannot be calculated accurately, the MAE are equivalent to the maximum standard absolute errors (MSAE):

    MSAE(f^k):=max0<x<1|f(x)f^k(x)|=max0<x<1|j=k+1fj(x)|,MSAE(g^k):=max0<x<1|g(x)g^k(x)|=max0<x<1|j=k+1gj(x)|.
  3. The maximum fractional errors (MFE):

    MFE(f^k)=max0<x<1|x2d2f^k(x)/dx2+xdf^k(x)/dxf^k(x)g^k(x)f^k(x)|×(100%),MFE(g^k)=max0<x<1|x2d2g^k(x)/dx2+xdg^k(x)/dx+Af^k2(x)g^k(x)|×(100%),

    represent the maximum values of the relative errors, defined as ratio of the absolute errors and the values of the HPM approximations in kth iteration.

Table 1:

Maximum iteration differences, absolute and fractional errors, as well as CPU time of the HPM approximations (Series A).

Order (k)MIDMAEMFE (%)CPU(s)
MID(f^k)MID(g^k)MAE(f^k)MAE(g^k)MFE(f^k)MFE(g^k)
13.75E−19.64E−12.11E−18.54E−16.22E+11.68E+00.125
26.61E−21.17E−19.03E−23.44E−13.92E+09.04E−10.234
37.71E−31.36E+03.07E−21.23E−11.53E+02.22E−10.390
41.21E−21.30E−18.80E−33.65E−24.67E−13.90E−20.765
51.05E−21.19E−11.97E−37.29E−31.24E−19.05E−31.203
69.14E−31.06E−14.12E−41.57E−32.43E−21.37E−31.906
78.03E−39.28E−27.84E−53.06E−45.54E−33.53E−42.718
87.15E−38.12E−21.34E−55.36E−51.04E−35.26E−53.796
96.43E−37.14E−22.08E−68.49E−61.19E−48.37E−65.061
105.83E−36.31E−29.42E−71.20E−62.35E−51.43E−66.889
Table 2:

Maximum iteration differences, absolute and fractional errors, as well as CPU time of the HPM approximations (Series B).

Order (k)MIDMAEMFE (%)CPU(s)
MID(f^k)MID(g^k)MAE(f^k)MAE(g^k)MFE(f^k)MFE(g^k)
13.24E−18.81E−12.48E−11.53E−16.84E+11.07E+00.265
21.40E−21.75E−11.06E−26.61E−22.30E+03.80E−10.530
36.66E−31.20E−19.60E−32.35E−28.03E−11.40E−11.077
46.49E−31.09E−11.03E−36.96E−32.34E−13.89E−21.733
55.41E−31.02E−17.36E−41.50E−35.79E−29.03E−32.826
64.98E−39.18E−29.94E−54.76E−44.52E−33.84E−44.529
74.49E−38.19E−21.39E−55.33E−51.25E−31.98E−46.388
84.07E−37.30E−23.60E−69.27E−64.17E−44.67E−57.688
93.29E−36.53E−24.48E−71.46E−69.35E−57.51E−68.935
102.75E−35.86E−21.55E−79.09E−71.53E−51.05E−612.497
Figure 2: Left: Graphs of the HPM approximations f^k(x)${\widehat{f}_{k}}(x)$ and g^k(x)${\widehat{g}_{k}}(x)$ obtained in Series A. Right: Appropriate approximations of the wave function θ(z)$\theta(z)$ and the potential ϕ(z)$\phi(z)$.
Figure 2:

Left: Graphs of the HPM approximations f^k(x) and g^k(x) obtained in Series A. Right: Appropriate approximations of the wave function θ(z) and the potential ϕ(z).

Figure 3: Left: Graphs of the HPM approximations f^k(x)${\widehat{f}_{k}}(x)$ and g^k(x)${\widehat{g}_{k}}(x)$ obtained in Series B. Right: Appropriate approximations of the wave function θ(z)$\theta(z)$ and the potential ϕ(z)$\phi(z)$.
Figure 3:

Left: Graphs of the HPM approximations f^k(x) and g^k(x) obtained in Series B. Right: Appropriate approximations of the wave function θ(z) and the potential ϕ(z).

The values of all approximation errors, computed for both series of HPM approximations up to 10th order, are shown in Tables 1 and 2. As can be seen, the accuracy of approximations in the last, 10th iteration is the order of 107, which is similar to HPM approximations obtained in solving some other problems, such as Fredholm and Volterra integral equations [32], [33]. Additionally, the usage of HPM technique increases the accuracy of approximations, in comparison to some other known solver techniques, applicable for the coupled SP system. For instance, the MFE in the last few iterations here are significantly less than 102%, i.e. the fractional error’s order of the numerical solution of the SP system based on two-directional fourth-order Runge–Kutta (RK4) algorithm [34].

Further, the last columns of the both Tables show the central processing unit (CPU) time for each iteration (calculated on a Core2 duo 2.16 GHz PC processor). Let us emphasize here that computation time is significantly reduced in comparison to some other methods existing in literature [35], [36]. For instance, the longest CPU time for HPM approximations (the 10th iteration of Series B) is some less than 12.5 s. On the other hand, the shortest computational time of the one-dimensional SP self-consistent solver, using the finite element technique by means of FEMLAB with MATLAB [36], is around 25 min on a Core2 duo 1.86 GHz PC processor. This indicates that CPU time for HPM is more than 120 times faster, what can be understood as an additional motivation for applying the HPM in solving the coupled SP system for MOS devices.

Finally, it can be easily seen that, in the case of A-series HPM approximations have a relative short computation time, but slightly slower convergence to the solutions of the SP-coupled equations. This is a consequence of distant values of initial solution g0(x) in comparison to solutions of (6)–(7). Conversely, HPM approximations of the B-series, with “closer” initial solution g0(x), converge somewhat faster, although the computation time is significantly longer. Convergences of both series are shown also in Figures 2 and 3. Here, the graphs of the HPM approximations f^k(x) and g^k(x), as well as the appropriate approximations of the wave function θ^k(z)=f^k(x(z)) and the potential ϕ^k(z)=x(z)g^k(x)/2ε0 are plotted, respectively.

5 Conclusion

The coupled system of SP equations that describes the behaviour of electrons in the n-type inversion layer has been solved by application of the HPM in this article. The case of electric quantum limit, with only the lowest energy subband occupied, has been considered. The SP system firstly has been simplified by introducing the appropriate dimensionless variables. Thereafter, by reducing (0,) interval to the finite one, the HPM approximate solutions that satisfy the requested boundary conditions have been obtained. The convergence of HPM approximations to exact solutions of SP system demonstrates the fact that the solutions of HPM are very accurate and are in excellent agreement with the exact solutions. Besides the efficiency and suitability, one more advantage of the HPM approach is that the solutions can be computed recursively with an arbitrary numerical precision, and the whole approximation procedure can be implemented by using the appropriate software. Moreover, the CPU time is significantly shorter than for some other methods that can be seen in literature.

References

[1] B. P. K. Yadav and A. K. Dutta, J. Semicond. Technol. Sci. 10, 203 (2010).10.5573/JSTS.2010.10.3.203Search in Google Scholar

[2] J. He, X. Xi, H. Wan, M. Dunga, M. Chan, et al. Solid-State Electron. 51, 433 (2007).10.1016/j.sse.2006.12.006Search in Google Scholar

[3] A. Chaudhry and J. N. Roy, Electron. 14, 86 (2010).Search in Google Scholar

[4] J. A. Pals, Quantization Effects in Semiconductor Inversion and Accumulation Layer, dissertation, Technical School of Eindhoven, Eindhoven, The Netherlands, 1972.Search in Google Scholar

[5] L. Bian G. Pang, S. Tang, and A. Arnold, J. Comput. Phys. 313, 233 (2016).10.1016/j.jcp.2016.02.025Search in Google Scholar

[6] F. Stern J. Comput. Phys. 6, 56 (1970).10.1016/0021-9991(70)90004-5Search in Google Scholar

[7] A. M. C. Serra and H. A. Santos, J. Appl. Phys. 70, 2734 (1991).10.1063/1.349389Search in Google Scholar

[8] I.-H. Tan, G. L. Snider, L. D. Chang, and E. L. Hu, J. Appl. Phys. 68, 4071 (1991).10.1063/1.346245Search in Google Scholar

[9] J. Heyl, M. W. Choptuik, and D. Shinkaruk, Phys. Rev. D 96, (2017). DOI: https://doi.org/10.1103/PhysRevD.96.103010.10.1103/PhysRevD.96.103010Search in Google Scholar

[10] L. Wang, D. Wang, and P. M. Asbeck, Solid State Electron. 50, 1732 (2006).10.1016/j.sse.2006.09.013Search in Google Scholar

[11] K. A. Berland, Superlattices Microstruct. 50, 411 (2014).10.1016/j.spmi.2011.08.003Search in Google Scholar

[12] M. Claus, S. Mothes, S. Blawid, and M. Schröter, J. Comput. Electron. 13, 689 (2014).10.1007/s10825-014-0588-6Search in Google Scholar

[13] M. J. van Dort, P. H. Woerlee, and A. J. Walker, Solid State Electron. 37, 411 (1994).10.1016/0038-1101(94)90005-1Search in Google Scholar

[14] Y. Ma, L. Liu, Z. Yu, and Z. Li, Microelectron. J. 31, 913 (2000).10.1016/S0026-2692(00)00097-5Search in Google Scholar

[15] F. Pregaldiny, C. Lallement, and D. Mathiot, Solid State Electron. 48, 781 (2004).10.1016/j.sse.2003.12.010Search in Google Scholar

[16] T. Kevkić, V. Stojanović, and D. Randjelović, Rom. J. Phys. 62, 122 (2017).Search in Google Scholar

[17] J.-H. He, Comput. Methods Appl. Mech. Engrg. 178, 257 (1999).10.1016/S0045-7825(99)00018-3Search in Google Scholar

[18] J.-H. He, Int. J. Non-Linear Mech. 35, 37 (2000).10.1016/S0020-7462(98)00085-7Search in Google Scholar

[19] J.-H. He, Appl. Math. Comput. 135, 73 (2003).10.1016/S0096-3003(01)00312-5Search in Google Scholar

[20] J.-H. He, Appl. Math. Comput. 156, 591 (2004).10.1016/j.amc.2003.08.011Search in Google Scholar

[21] J.-H. He, Comput. Math. Appl. 57, 410 (2009).10.1016/j.camwa.2008.06.003Search in Google Scholar

[22] A. M. A. El-Sayed, A. Elsaid, I. L. El-Kalla, and D. Hammad, Appl. Math. Comput. 218, 8329 (2012).10.1016/j.amc.2012.01.057Search in Google Scholar

[23] A. A. Hemeda, Appl. Math. Sci. 96, 4787 (2012).Search in Google Scholar

[24] M.-F. Zhang, Y.-Q. Liu, and X.-S. Zhou, Therm. Sci. 19, 1167 (2015).10.2298/TSCI1504167ZSearch in Google Scholar

[25] M. Zeb, T. Haroon, and A. M. Siddiqui, U.P.B. Sci. Bull. Series A 76, 179 (2014).Search in Google Scholar

[26] P. K. Roy and A. Mallick, Alexandria Engrg. J. 55, 2269 (2016).10.1016/j.aej.2016.05.020Search in Google Scholar

[27] K. Grysa and A. Maciag, Int. J. Heat Mass Tran. 100, 627 (2016).10.1016/j.ijheatmasstransfer.2016.04.103Search in Google Scholar

[28] V. Stojanović, T. Kevkić, G. Jelić, and D. Randjelović, U.P.B. Sci. Bull. Series A 80, 119 (2018).Search in Google Scholar

[29] T. Kevkić, V. Stojanović, and D. Petković, Rom. Rep. Phys. (to appear) http://rrp.infim.ro/IP/A352.pdf.Search in Google Scholar

[30] D. R. Islamov, V. A. Gritsenko, T. V. Perevalov, O. M. Orlov and G. Ya, Krasnikov. Appl. Phys. Lett. 109, 052901 (2016).10.1063/1.4960156Search in Google Scholar

[31] J. A. Barrett, Computational and Analytical Methods for the Simulation of Electronic States and Transport in Semiconductor Systems, dissertation, Anglia Ruskin University, Cambridge, England 2014.Search in Google Scholar

[32] E. Hetmaniok, I. Nowak, D. Slota, and R. Witula, Appl. Math. Comput. 218, 10717 (2012).10.1016/j.amc.2012.04.041Search in Google Scholar

[33] E. Hetmaniok, D. Slota, and R. Witula, Appl. Math. Lett. 26, 165 (2013).10.1016/j.aml.2012.08.005Search in Google Scholar

[34] M. J. Hargrove, A. K. Henning, J. A. Slinkman, and C. E. Hembree, Self-Consistent Solution of the Schrödinger and Poisson Equations Applied to Quantum Well Heterostructures, Technical Report, Thayer School of Engineering, Dartmouth College, Hanover, MS, USA, 1998.Search in Google Scholar

[35] H.-Y. Jiang and P.-W. Zhang, J. Comput. Math. 24, 401 (2006).Search in Google Scholar

[36] M. K. Alam, A. Alam, G. Rabbani, and Q. D. M. Khosru, An Accurate and Fast Schrödinger-Poisson Solver using Finite Element Method, in: Proceedings of the 18th IASTED International Conference on Modelling and Simulation (MS 2007), Montreal, Canada, pp. 246–249, 2007.Search in Google Scholar

Received: 2018-11-03
Accepted: 2019-01-11
Published Online: 2019-02-12
Published in Print: 2019-06-26

©2019 Walter de Gruyter GmbH, Berlin/Boston

Downloaded on 28.3.2024 from https://www.degruyter.com/document/doi/10.1515/zna-2018-0495/html
Scroll to top button