Relocation dynamics and multistable switching in semiconductor superlattices

Journal of Computational and Applied Mathematics 204 (2007) 18 – 24 www.elsevier.com/locate/cam

Relocation dynamics and multistable switching in semiconductor superlattices Guido Dell’Acqua∗ , Luis L. Bonilla, Ramón Escobedo Grupo de Modelización y Simulación Numérica, Escuela Politécnica Superior, Universidad Carlos III de Madrid, Av. de la Universidad 30, 28911 Leganés, Madrid, Spain Received 15 July 2005; received in revised form 10 December 2005

Abstract A numerical study of electric field domain relocation during slow voltage switching is presented for a spatially discrete model of doped semiconductor superlattices. The model is derived from the Poisson’s equation and the charge continuity equation. It consists of an Ampère equation for the current density and a global summatory condition for the electric field and it has been particularly effective in the prediction and reproduction of experimental results. We have designed a fast numerical scheme based on the use of an explicit expression for the current density. The scheme reproduces both previous numerical and experimental results with high accuracy, yielding new explanations of already known behaviors and new features that we present here. © 2006 Elsevier B.V. All rights reserved. MSC: 65P40; 65L07; 81T10; 81T80 Keywords: Numerical simulations; Quantum transport; Discrete model; Multistability; Weakly coupled superlattices; Domain relocation

1. Introduction Semiconductor superlattices (SLs) are essential ingredients in fast nanoscale oscillators, quantum cascade lasers and infrared detectors. Quantum cascade lasers are used to monitor environmental pollution in gas emissions, to analyze breath in hospitals and in many other industrial applications. A SL is formed by growing a large number of periods with each period consisting of two layers, which are semiconductors with different energy gaps but having similar lattice constants, such as GaAs and AlAs. The conduction band edge of an infinitely long ideal SL is modulated so that it looks like a one-dimensional (1D) crystal consisting of a periodic succession of a quantum well (GaAs) and a barrier (AlAs). Vertical charge transport in a SL subject to strong electric fields exhibits many interesting features, and it is realized experimentally by placing a doped SL of finite length in the central part of a diode (forming a n+ –n–n+ structure) with contacts at its ends. Depending on the bias condition, the SL configuration, the doping density, the temperature and other control parameters, the current through the SL and the electric field distribution inside the SL display a great variety of nonlinear phenomena such as pattern formation, current self-oscillations and chaotic behavior [2]. In 1998, Luo et al. reported experimental results on how a domain wall relocates if the voltage across the SL is suddenly changed [3]. These experiments have been explained recently by numerical simulating a discrete model with ∗ Corresponding author.

E-mail address: [email protected] (G. Dell’Acqua). 0377-0427/$ - see front matter © 2006 Elsevier B.V. All rights reserved. doi:10.1016/j.cam.2006.04.025

G. Dell’Acqua et al. / Journal of Computational and Applied Mathematics 204 (2007) 18 – 24

19

the tunneling current density given by a constitutive relation in terms of the local electric field and the electron densities at adjacent wells [1]. Here we reproduce the model, we derive an explicit expression for the total current density which leads to an equivalent model, and we solve it with an effective numerical algorithm. We present also new features of the model which can help us to explain already known behaviors of the physical device. 2. The model The model consists of the following Poisson and charge continuity equations: e Fi − Fi−1 = (ni − ND ), ε

(1)

dni = Ji−1→i − Ji→i+1 , dt

(2)

for the average electric field −Fi and the two-dimensional (2D) electron density ni at the ith SL period (which starts at the right end of the (i − 1)th barrier and finishes at the right end of the ith barrier), with i = 1, . . . , N. Here ND , ε, −e and eJ i→i+1 are the 2D doping density at the ith well, the average permittivity, the electron charge and the tunneling current density across the ith barrier, respectively. The SL period is l = d + w, where d and w are the barrier and well widths, respectively. Time-differencing Eq. (1) and inserting the result in Eq. (2), we obtain the following form of Ampere’s law: ε dFi + Ji→i+1 = J (t), e dt

i = 0, . . . , N

(3)

which may be solved with the bias condition for the applied voltage V (t): N V (t) 1  . Fi = N +1 (N + 1)l

(4)

i=0

The space-independent unknown function eJ (t) is the total current density through the SL. We use a simplified constitutive relation for the tunneling current density across barriers in terms of the local electric field and the electron densities at adjacent wells [2]: J0→1 = F0 ,

(5)

   2   v (f ) (Fi ) m ∗ kB T eF i l h¯ ni+1 Ji→i+1 = ln 1 + exp − ni − exp −1 , l kB T m∗ k B T h¯ 2 

JN→N+1 = FN



nN . ND

i = 1, . . . , N − 1, (6) (7)

Here  is the contact conductivity (assumed to be the same at both contacts for simplicity), m∗ the effective mass, T the temperature, kB and h¯ are the Boltzmann and Planck constants, respectively, and v (f ) is the drift velocity:

v (f ) (Fi ) =

n  j =1

Ti () = √

h¯ 3 l(C1 + Cj ) Ti (EC1 ) 2m∗2 , (EC1 − ECj + eF i l)2 + (C1 + Cj )2

2 2 (k 2 + 2 )−1 (k 2 + 2 )−1 16ki2 ki+1 i i i i+1 i

−1 −1 −1 2i d (w + −1 i−1 + i )(w + i+1 + i )e

2m∗ , h¯ ki+1 = 2m∗ [ + e(d + w)Fi ], hk ¯ i=

,

(8)

(9) (10) (11)

20

G. Dell’Acqua et al. / Journal of Computational and Applied Mathematics 204 (2007) 18 – 24

Table 1 Parameters of the SL N 40

ND (cm−2 )

w/d (nm/nm)

m∗ (10−32 kg)

EC 1

E C2

EC 3

(meV)

(meV)

(meV)

(meV)

Vb (V)

1.5 × 1011

9.0/4.0

8

8.43

44

180

410

0.982



60 5 4

(FM,JM)

3

40 JM (A/cm2)

2 1 0 0

2

4

6

8

10

20

0

0

100

200 FM (kV/cm)

300

400

Fig. 1. Constitutive relation of the tunneling current density for doping density, showing the calculation of FM ≈ 3.945 kV/cm and JM ≈ 3.1269 A/cm2 for T = 5 K and the SL values of Table 1.

h ¯ i−1 =  h ¯ i=

(12)

  ewF i 2m∗ eV b − − , 2 

h¯ i+1 =

  w 2m∗ eV b + e d + Fi −  , 2

(13)

    3w 2m∗ eV b − e d + Fi −  . 2

(14)

Here Cj indicates the jth subband in a well, ECj is the energy level, Cj is the scattering width, Ti is the dimensionless transmission probability across the ith barrier, and eV b is the barrier height in absence of potential drops. Typical values of these parameters are shown in Table 1 (Fig. 1). These formulae have the advantage over pure numerical computations of being analytical and they can be used to develop the theory further. Given a known configuration of a sample used in experiments, our formulae allow calculation of constitutive relations that can be easily used to determine the dynamical behavior of the SL. For numerical treatment, it is convenient to render the equations dimensionless. We have used the following notation, introducing the typical scales of each physical magnitude: Fi = FM Ei ,

ni = ND n˜ i ,

v(Fi ) = vM v(E ˜ i ),

Ji→i+1 = JM J˜i→i+1 ,

V = V0 ,

 = c , ˜

J = JM J˜,

x = x0 x, ˜

t = t0 t˜.

(15) (16)

G. Dell’Acqua et al. / Journal of Computational and Applied Mathematics 204 (2007) 18 – 24

21

Table 2 Typical scales for T = 5 K FM (kV/cm)

eJ M (A/cm2 )

vM (m/s)

x0 (nm)



0

(–)

(–)

–

–

JM l ND

εFM l eN D

eN D εFM

3.945

3.127

1.691

2.494

5.212

m∗ kB T h¯ 2 ND 0.111

c ( m)

V0 (V)

lF M ev M ND

FM Nl

12.62

0.205

The values of FM and JM are calculated as the coordinates of the first relative maximum of the function Ji→i+1 (Fi , ni , ni+1 ) = Ji→i+1 (Fi , ND , ND ). We have also defined the following dimensionless parameters: =

eN D , εFM

0 =

m ∗ kB T , h¯ 2 ND

1 = e1/0 − 1,

a=

elF M , kB T

(17)

where is the dimensionless doping and 0 ?1 (0 >1) denotes the high (low) temperature limit behavior of the system. The two other parameters are just for notation. The values of these parameters corresponding to the SL described in Table 1 are given in Table 2. Then the dimensionless model can be written (dropping tildes) as dEi (t) + Ji→i+1 (Ei , ni , ni+1 ) = J (t), dt Ei − Ei−1 + 1, i = 1, . . . , N, ni = N 

Ei (t) = (N + 1)(t),

i = 0, . . . , N,

(18) (19)

t 0,

(20)

i=0

Ji→i+1 = v(Ei ){ni − 0 ln[1 + e−aE i (eni+1 /0 − 1)]},

i = 1, . . . , N − 1,

(21)

J0→1 = E0 ,

(22)

JN→N+1 = EN nN .

(23)

3. Numerical solution The efficiency of our algorithm is based on that doing the sum in (18) from i =0 to N, we obtain an explicit expression for the total current density: J (t) =

N d(t) 1  + Ji→i+1 (Ei (t)). dt N +1

(24)

i=0

The initial system can be solved equivalently by solving, for i = 0, . . . , N, 1 dEi (t) d(t) = + dt dt N +1

N 

Jj →j +1 (Ej (t)) −

j =0,j =i

N Ji→i+1 (Ei (t)), N +1

(25)

together with the initial condition Ei (0) = (0), i = 0, . . . , N. Note that the bias condition is preserved: summing again in (25) yields ⎛ ⎞ N  N N N N  1 ⎝  d  ⎠ Ei − (N + 1) = Jj →j +1 − Ei , (26) dt N +1 N +1 i=0

which is equal to zero, so

i=0

N

i=0 Ei (t) = (N

j =0,j =i

i=0

+ 1)(t) plus a constant, which must be zero to fulfill the initial condition.

22

G. Dell’Acqua et al. / Journal of Computational and Applied Mathematics 204 (2007) 18 – 24

J (A/cm2)

3

1

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

2.7 2.4 2.1

0

3 J (A/cm2)

2

0.2 20

21

0.4 22

23

0.6 24

25

0.8 26

27

1 28

29

1.2 30

31

1.4 32

33

1.6 34

35

1.8 36

37

2 38

39 40

2.7 2.4 2.1

2

2.2

2.4

2.6

2.8

3

3.2

3.4

3.6

3.8

4

V (V) Fig. 2. Current–voltage diagram for T = 5 K and  = 25.2  m showing stationary branches Bi , i = 1, . . . , N , and regions of multistability.

We have used explicit and implicit methods to solve this system: an order 1 Euler method, an embedded Runge–Kutta method of order 7(8) with step-size control and error estimate, and BDF methods of order 1–4. Implicit methods are solved by means of Newton–Raphson iterations. Results. Fig. 2 shows the current–voltage characteristic curve for the SL values of Tables 1 and 2. For constant , the stable field profiles {Ei } are time-independent, step-like and increasing with i: typically they consist of two flat regions called electric field domains separated by an abrupt transition region called a domain wall or charge monopole [2]. The electric field profiles of each branch of solutions in Fig. 2 differ in the location of their domain wall: counting branches in the direction of increasing voltage, the profiles of the jth branch have their domain wall located in the (N − j + 1)th SL period. Notice that for certain values of the voltage, several branches with different current are possible (multistability). The central part of the first branches (from B1 to B21 in Fig. 2) are regions of monostability. The upper and lower limits of these branches overlap with the previous and the following ones, defining intervals of bistability. The 22th branch is the last having a region of monostability and the first whose upper limit overlaps with the two following branches, describing a region of tristability. Branches B22 –B39 are all bistable/tristable, and the last branch B40 has the three types of behavior: it is tristable in the lower limit, bistable in the central part and monostable in the upper limit. If we switch the voltage from a value Vini corresponding to one branch to a final value Vend =Vini +V corresponding to different branches, the domain wall has to relocate in a different SL period. During switching, V (t) = Vini + V˙ t, with V˙ = V /t, and t is the ramping time. Fig. 3 shows the current density and the electric field evolution during a relocation with t = 0. In this case, the current density exhibits an interval of double peaks, followed by an interval of single peaks, between both stationary states, accompanied in the electric field distribution by the nucleation of a dipole wave at the injecting contact, which crosses the sample and is finally absorbed by the high field domain. See Fig. 3(B). All these features are as observed in the experiments [1,3–6]. For larger values of V , this scenario is repeated a number of times equal to the number of branches crossed in V , provided t is greater than a critical value tc . Fig. 4 shows that depending on whether or not t is greater than tc the curve (V (t), J (t)) has time to describe the same number of pics than branches the voltage crosses. The system evolves with (t) near a stationary branch until the upper limit is reached; then it falls to the next branch. This change of branch is produced by means of the relocation of the electric field domain by injecting a dipole wave at the cathode. During this process, the current decreases to low values to allow the emission and the trip of the dipole

G. Dell’Acqua et al. / Journal of Computational and Applied Mathematics 204 (2007) 18 – 24

23

3.2 2.8 +0.2 V

Field

J (A/cm2)

2.4 2

40

V=1.0 V

1.6

80 60 40 20 35

1000

30 25 20 ll 15

1.2

800 600

we

0.8 -200

0

200

(A)

400

600

400

10 5

800 1000 1200

200

e

tim

0

(B)

t (ns)

3

3

2.5

2.5 J (A/cm2)

J (A/cm2)

Fig. 3. Relocation of the electric field wall for T = 5 K and  = 25.2  m. For t < 0, V = 1.0 V in B10 . At t = 0, we add V = 0.2 V and reach B12 . (A) Current density during relocation: J (t) starts at 2.44 A/cm2 , then falls to low values at t ≈ 20 ns, it describes the double/single peaks pattern observed in the experiments, and finally grows again to a stationary value of 2.47 A/cm2 . (B) Electric field, showing the domain relocation from well 31 to well 29.

2

1.5

1.5

1 (A)

2

0.85

1 0.95

1.05

1.15 V(V)

1.25

1.35

0.85 (B)

0.95

1.05

1.15

1.25

1.35

V(V)

Fig. 4. Current density curve (V (t), J (t)) (solid line) during the voltage switching from Vini = 0.83 V to Vend = 1.37 V for a large ramping time and the half: (A) tramp = 30 s and (B) tramp = 15 s. Also depicted (small circles) are the stationary branches of the I–V characteristic of Fig. 2.

wave. Once the wave is absorbed by the high field domain, the current returns to a stationary value located in the next branch. See Fig. 4(A) for t = 30 s and (B) for t = 15 s. Voltage values are Vini = 0.83 V and Vend = 1.37 V. Two different behaviors are observed: • When t > tc , the return of the current density to the I–V branch takes place in a region of monostability, allowing the curve to follow the branch again, until the upper limit is reached. This is repeated 4 times, and it is accompanied with the corresponding five emissions of a dipole wave. Fig. 5(A) shows one of these five relocation processes by injection of a wave corresponding to Fig. 4(A). • When t < tc , the curve (V (t), J (t)) returns to a range of bistability and falls to the lower branch. Fig. 4(B) shows that after having been following the third branch, the curve falls to low values of J to allow the injection of a dipole wave, as before. However, when it returns to stationary values, the curve has not time to reach the 4th branch before entering into the bistability region of branches 4 and 5. Then the curve goes to the 5th branch, and

24

G. Dell’Acqua et al. / Journal of Computational and Applied Mathematics 204 (2007) 18 – 24

80

80

60

60

40

40

20

20

3000 3500 4000 4500 5000 5500 (A)

5

10

15

20

25

30

35

40

2000 4000 6000 8000 10000 12000 14000 (B)

5

10

15

20

25

30

35

40

Fig. 5. Electric field distributions during voltage switching for both ramping times shown in Fig. 4. (A) Example of domain relocation by means of the absorption of the electric dipole wave. (B) Switching process for tramp = 15 s, showing the absence of the 4th and 6th waves.

the 4th branch has been dropped. The result is that the current density curve exhibits one maximum less, and the electric field has injected only four dipole waves during this shorter switching time. The same thing occurs with the 6th branch, which is also dropped; see again Fig. 4(B). In conclusion, we have presented a new formulation of a domain relocation problem and the numerical simulations which gives a more satisfactory explanations of a previously observed behavior of the device, in terms of the stability intervals of the applied voltage. These results have been obtained with a fast and effective numerical algorithm based on this new formulation. Acknowledgments We thank S.W. Teitsworth for fruitful discussions. This work has been supported by the MCyT Grant BFM200204127-C02-01, by the Spanish MECD grant MAT2005-05730-C02-01 and by the European Union under Grant HPRNCT-2002-00282. R. Escobedo has been supported by a postdoctoral grant awarded by the Consejería de Educación of the Autonomous Region of Madrid. References [1] A. Amann, A. Wacker, L.L. Bonilla, E. Schöll, Dynamic scenarios of multistable switching in semiconductor superlattices, Phys. Rev. E 63 (2001) 066207. [2] L.L. Bonilla, Theory of nonlinear charge transport, wave propagation and self-oscillations in semiconductor superlattices, J. Phys.: Condens. Matter 14 (2002) R341–R381. [3] K.J. Luo, H.T. Grahn, K.H. Ploog, Relocation time of the domain boundary in weakly coupled GaAs/AlAs superlattices, Phys. Rev. B 57 (1998) R6838–R6841. [4] M. Rogozia, S.W. Teitsworth, H.T. Grahn, K.H. Ploog, Statistics of the domain-boundary relocation time in semiconductor superlattices, Phys. Rev. B 64 (2001) 041308 (R). [5] M. Rogozia, S.W. Teitsworth, H.T. Grahn, K.H. Ploog, Relocation dynamics of domain-boundaries in semiconductor superlattices, Phys. Rev. B 65 (2002) 205303. [6] M. Rogozia, S.W. Teitsworth, H.T. Grahn, K.H. Ploog, Time distribution of the domain-boundary relocation in superlattices, Physica B 314 (2002) 427–430.