Bifurcation behavior of a superlattice model

BIFURCATION BEHAVIOR OF A SUPERLATTICE MODEL∗ M. MOSCOSO†,J.GAL ´ AN‡,AND L. L. BONILLA† SIAM J. APPL. MATH.c 2000 Society for Industrial and Applied Mathematics Vol. 60, No. 6, pp. 2029–2057 Abstract. We present a complete description of the stationary and dynamical behavior of semiconductor superlattices in the framework of a discrete drift model by means of numerical continuation, singular perturbation analysis, and bifurcation techniques. The control parameters are the applied DC voltage (φ) and the doping (ν) in nondimensional units. We show that the organizing centers for the long time dynamics are Takens–Bogdanov bifurcation points in a broad range of parameters and we cast our results in a φ-νphase diagram. For small values of the doping, the system has only one uniform solution where all the variables are almost equal. For high doping we find multistability corresponding to domain solutions and the stationary solutions may exhibit chaotic spatial behavior. In the intermediate regime of νthe solution can be time-periodic depending on the bias. The oscillatory regions are related to the appearance and disappearance of Hopf bifurcation tongues which can be sub- or supercritical. These results are in good agreement with most of the experimental observations and also predict new interesting dynamical behavior. Key words. Takens–Bogdanov bifurcation, current instabilities, numerical continuation, amplitude equation AMS subject classifications. 58F40, 34A12, 34E05 PII. S0036139998336754 1. Introduction. The term “superlattice” (SL) was proposed in 1969 by Esaki and Tsu to refer to a man-made semiconductor structure, created by the periodic deposition of different semiconductors [12, 13]. This results in a periodic array of quantum wells and barriers formed by the two distinct elements where the main charge transport mechanism is resonant tunneling (sequential or coherent). This gives rise to a strong nonlinearity in the transport process, with several regions of negative differential conductivity (NDC). The electronic properties can be further adjusted by means of an externally applied voltage. The optical and transport properties of an SL can be designed and controlled at will, by changing the composition, dimensions, and doping concentration of the regions during the growing process. Depending on all these conditions (number of quantum wells, doping density, applied voltage, etc.), interesting spatiotemporal structures have been observed in experiments. They include spatially uniform electric field, nonuniform electric field (domains), and self-sustained oscillations due to domain motion [5, 11, 16]. The current-voltage characteristic curve (I−φ) presents different structures in each case. Recently, it has been shown that devices based on semiconductor SLs can be used as tunable GHz oscillators and for higher temperatures [21] even spontaneous or driven chaotic current oscillations have been reported [47, 26]. Thus doped semiconductor SLs are a very interesting example of a nonlinear dynamical system with a large number of degrees of freedom. In principle, transport in semiconductor SLs should be investigated at the quantum mechanical level. A general approach is to derive and solve quantum kinetic ∗Received by the editors April 3, 1998; accepted for publication (in revised form) December 8, 1999; published electronically June 15, 2000. This work was partially supported by Junta de Andaluc´ıa and by DGES grant PB97-0088. http://www.siam.org/journals/siap/60-6/33675.html †Escuela Polit´ecnica Superior, Universidad Carlos III de Madrid, Butarque 15, Legan´es 28911, Spain ([email protected], b[email protected]). ‡Departamento de Matem´atica Aplicada II, Escuela Superior de Ingenieros de Sevilla, Camino de los Descubrimientos s/n, Sevilla 41092, Spain ([email protected]). 2029 Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 2030 M. MOSCOSO, J. GAL´ AN, AND L. L. BONILLA equations for nonequilibrium Green functions (NGF) [18]. A starting point of the NGF approach is to use a tight-binding Hamiltonian and a basis of Wannier functions localized at the different SL wells [44]. This setting implies that the spatial coordinates are discrete numbers representing well indices. In this representation, possible spatially nonlocal terms of the kinetic equations are replaced by sums, which are easier to handle. However promising, the NGF approach has been worked out only for stationary states and is subject to severe restrictions such as considering a constant electric field throughout the SL [44]. These restrictions preclude explaining the spatiotemporal structures mentioned above. Despite this, the NGF approach has already scored important successes such as a better quantitative calculation of static current-voltage characteristics of SL. Furthermore, this approach has been used to ascertain in which cases simpler theories hold. The main simpler theories are these (see Figure 1 of [44]): (i) Semiclassical calculations of miniband transport using the Boltzmann transport equation [30] or simplifications thereof, such as hydrodynamic [7] or drift-diffusion [38] models. These calculations hold for strongly coupled SL at low fields. In the miniband transport regime, electrons traverse the whole SL miniband thereby performing Bloch oscillations and giving rise to NDC for large enough electric fields [13]. The latter may cause self-sustained oscillations of the current due to recycling of charge dipole domains as in the Gunn effect of bulk n-GaAs [7]. (ii) Wannier–Stark (WS) hopping transport in which electrons move parallel to the electric field through scattering processes including hopping transitions between WS levels [35]. Calculations in this regime hold for intermediate fields, larger than those corresponding to collisional broadening of WS levels, but lower than those corresponding to resonant tunneling. (iii) Sequential tunneling calculations valid for weakly coupled SL (coherence length smaller than one SL period) at basically any value of the electric field [22, 1, 43]. On the other hand, the description of electric field domains and self-sustained oscillations in SL has been made by means of discrete drift models. These models use simplified forms of the sequential tunneling current through SL barriers and discrete forms of the charge continuity and Poisson equations [34, 3, 2]. Discrete drift models yield good descriptions of nonlinear phenomena in weakly coupled SL. Recently, these models have been derived from sequential tunneling microscopic descriptions [6]. We shall now briefly describe the discrete drift model of [2], indicating connections to the more detailed sequential tunneling model of [1] whenever appropriate. For a full description of this model and the derivation of discrete drift-diffusion models, see [6]. We consider that each quantum well jof the SL may be characterized by the local value of the electric field, Ej, and by the two-dimensional charge density, nj,at the well. In more detailed models, it is necessary to distinguish between drops of the electric potential at the barriers from those at the wells [1]. In the usual discrete drift model, these potential drops (proportional to the electric field) are not distinguished. However, we need to distinguish momentarily between potential drops at barriers and wells in order to understand the Poisson equation written below. Let us label SL barriers starting from the zeroth barrier separating the injecting contact from the SL. Then the jth barrier separates wells jand j+1. Ejis the average field at the region comprising the jth well and the jth barrier. For weakly coupled SL, electrons take a long time before they can escape well jby tunneling through the jth barrier to the adjacent well, j+ 1 [3, 2, 1]. Therefore at a high applied electric field, they tend to accumulate just at the end of the well, where the barrier begins. Thus we consider Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php BIFURCATION IN SUPERLATTICES 2031 that the electron density is singularly concentrated on a plane located at the end of the jth well [1, 6]. Then the averaged Poisson equation relating the two-dimensional charge density at well jto the electric fields at the barriers j−1 and j(which bound the well) is Ej−Ej−1=ν(nj−1) .(1.1) We are using the nondimensional variables defined in [21], where characteristic values of SLs used in experiments are given. The parameter νis the dimensionless doping at the wells. Notice that placing the electron plane at different locations inside the well will result in different forms of (1.1), corresponding to different discretizations of a spatial derivative. It is well known that differencing schemes may strongly affect the dynamics of a system described by ordinary or partial differential equations; see [46, 39, 40] and references therein. Our discrete drift model is a model on its own, not a numerical scheme for a continuous model of transport in weakly doped SL. To the best of our knowledge, a time-dependent, spatially continuous, quantum model of vertical transport in weakly coupled SL does not yet exist. Another equation of the discrete drift model is charge continuity, which establishes that the rate at which njchanges is proportional to the tunneling current arriving to the well jfrom j−1 minus the tunneling current from well jto well j+ 1. This equation is equivalent to Ampere’s [3, 2]: the total current density through the SL, I, is the sum of Maxwell’s displacement current, dEj/dt, and the tunneling current from well jto well j+1,Jj→j+1: dEj dt +Jj→j+1 =I.(1.2) The charge continuity equation may be obtained by differentiating (1.1) with respect to time and then substituting (1.2) in the result. The constitutive relation linking Jj→j+1 to Ejand njmay be obtained from quantum mechanical or phenomenological considerations; see [43] for discussions. Except at low biases where a reverse tunneling current should be included, the following constitutive relation [3] includes all the relevant physics: Jj→j+1 =njv(Ej).(1.3) Here v(Ej) has maxima at the (resonant) fields where energy levels of neighboring quantum wells are aligned. In this paper, we only consider voltages corresponding to fields between two successive resonances (not including the one at zero field). Then we can use a velocity curve with a single minimum as depicted in Figure 1, where the peak at low field values has been omitted. With the dimensionless units of [2], v(E) has a maximum vM=v(1) = 1 and a minimum vm=v(Em)≈0.32, with Em>1 (see Figure 2(a)). The negative slope between maximum and minimum is the NDC mentioned above. This is a purely quantum mechanical effect (due to the alignment or misalignment of energy levels) and it constitutes the essential quantum effect which necessarily has to be included in our model. In our calculation, v(E)isa phenomenological curve that reproduces the experimental current voltage curves and is consistent with the usual expressions in the literature [22, 1]. It has been computed with the following equation: v(E)= 1 voA0+A1 (E−E1)2+Γ 2 1 +A2 (E−E2)2+Γ 2 2,(1.4) Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 2032 M. MOSCOSO, J. GAL´ AN, AND L. L. BONILLA 0.0 1.0 2.0 E 0.0 0.5 1.0 1.5 2.0 V(E) v M v m E m Fig. 1.Velocity curve (v(E)) as a function of the electric field. It has a maximum vM= v(1)=1, and a minimum vm=v(Em)≈0.32, with Em>1. The negative slope between maximum and minimum is the negative differential conductivity mentioned in the text. where E1= 1 and E2=2.7 are the values of the energy levels in the well, A1=0.011 and A2=0.05 are proportional to the transmission coefficients across the barrier from the first and second energy levels, Γ1=0.3 and Γ2=0.6 are the broadening terms related to scattering, and A0and v0are normalizing factors that have been taken so that v(0) = 0 and v(1) = 1. In [6], it has been shown that the general form of the tunneling current density at low temperatures is Jj→j+1 =njv(Ej)−D(Ej)(nj+1 −nj). Here the discrete diffusion term is due to reverse tunneling current acting at low fields: the diffusion coefficient D(Ej) decreases rapidly to zero as Ejincreases. At the voltages we are interested in, D= 0, and we obtain (1.3). To complete the description of the model, we add the condition that the voltage drop across the SL is kept constant, 1 N N  j=1 Ej=φ,(1.5) where φis the external applied voltage and Nis the total number of quantum wells in the SL. Furthermore, we need a boundary condition for E0, the average field before the SL, (E1−E0)=cν(1.6) (equivalent to n1=1+c), where c>0 is a third parameter in our model. Equation (1.6) is a simplified model for the contact region before the SL: the most general boundary condition on physical grounds is a relation between n1,E1, and I[6]. The physical origin of cis clear for an electronically doped structure (n-i-n): due to the different electron concentrations at each side of the first barrier, some charge will be transferred from the contact to the first quantum well. This creates a small dipole field which cancels the electron flow. More complicated boundary conditions may be Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php BIFURCATION IN SUPERLATTICES 2033 proposed based upon models of the injecting and receiving contact regions outside the SL [1, 6]. There are wavefront solutions of (1.1)–(1.6) traveling from left to right across the SL for appropriate doping. A wavefront has a step-like profile of the electric field, which increases from a constant value to another one as the well index increases. The domain wall joining regions (electric field domains) where the field is constant contains an excess of electrons, and the wavefront is called a (charge) monopole [5, 21]. The boundary condition at the injecting contact selects whether self-sustained oscillations of the current are due to recycling of monopoles or dipoles (pulses of the electric field) [36]. The condition (1.6) selects monopoles as responsible of self-sustained current oscillations, which is consistent with the highly doped contact regions used in the experiments reported in [21]; see [36]. A boundary condition at the collector region (receiving contact) is not needed for the model (1.1)–(1.3). Only at low enough electric fields does the presence of the diffusion term in the constitutive relation require using another boundary condition in addition to (1.6) [6]. Equations (1.1)–(1.6) constitute the simplest model that retains the essential physics of vertical sequential transport in an SL, allows a thorough analytical study [42, 5], and successfully compares with the experiments reported in [16, 21]. Recently, the bifurcation diagrams obtained in this paper were used to interpret experimental observations of current self-oscillations in undoped, photoexcited SL [27]. Furthermore, B. Sun et al. observed how a transverse magnetic field induced self-sustained oscillations of the current in a doped SL which had static current-voltage characteristics at zero magnetic field [41]. These authors modified the velocity curve v(Ej)of the constitutive relation (1.3) to account for the influence of the magnetic field. Then they showed that simulations of the model (1.1)–(1.6) yield excellent agreement with observed results [41]. Other models proposed recently include Laikhtman’s equivalent circuit model [25], which does not use Poisson’s equation. For this model, the only coupling between wells comes from the bias condition. This produces spurious nonmonotonic stable stationary field profiles. Prengel, Wacker, and Sch¨oll [34] tried a rough elementary derivation of the constitutive relation for a more complicated two-level model of transport between wells. This model (or successive elaborations thereof, including randomness in the well doping) also explains the shape of stationary electric field profiles in the SL and the self-sustained oscillations (see [43]), but it is too complicated to perform a detailed analysis. Numerical studies of the phase diagram (bifurcation set) for this model are presented in [31, 32, 37]. In this work we undertake a complete description of the phase diagram of the model (1.1)–(1.6) partially covered in previous works [3, 15, 21, 42]. We will show that the dynamical behavior is very rich and the long time dynamics is governed by Takens–Bogdanov bifurcation points in a broad range of parameters. This fact was already noted by Patra et al. in their models [31, 33], Patra, Schwarz, and Sch¨oll [32], and Sch¨oll et al. [37], although they had to include randomness in the doping to obtain similar dynamic behavior. We have found new interesting phenomena that could be relevant to the analysis and design of future devices: the existence of windows of oscillatory behavior as the applied voltage is swept, subcritical as well as supercritical Hopf bifurcations, and the possibility of spatially chaotic stationary solutions. The rest of the paper is structured as follows. In section 2, the stationary and periodic solutions of the model are investigated. In section 3 we present numerical continuation results in the form of a phase diagram. In section 4 we present an analytical study of the onset and end of the oscillatory solutions by means of an Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 2034 M. MOSCOSO, J. GAL´ AN, AND L. L. BONILLA amplitude equation, and we summarize our results in section 5. The details of the calculations of section 4 are given in the appendices. In Appendix A we study the linearization of the model and introduce some notation. The linear system with a constraint has been solved in detail because this problem appears again in the hierarchy of linear equations. In Appendix B we derive the solvability condition that ensures the absence of secular terms. In Appendix C we derive the amplitude equation with a multiple scale expansion. 2. Bifurcation diagram. In this section we study the equilibria and the time periodic solutions of our model. It is convenient to eliminate the charge densities nj from the Ampere equation (1.2) by means of Poisson’s equation (1.1). The result is dEj dt +1+Ej−Ej−1 νv(Ej)=I.(2.1) The N+ 2 equations, (2.1), (1.5), and (1.6), contain N+ 2 unknowns (E0,E1,E2, ..., EN,I). Thus these equations and a set of initial conditions, Ej(0), constitute a well-posed problem. The main parameters of our system are νand φ. They are proportional to the doping density in the quantum wells and to the external DC voltage applied between the two SL ends, respectively. Note that the doping density is fixed once the growth process has finished, while the external voltage can be easily controlled in the laboratory. 2.1. Construction of stationary solutions and their current-voltage characteristics. It is straightforward to show that the stationary solutions of the system can be found by analyzing the following discrete map: Ej−1=f(Ej,ν,I)(2.2) with f(E,ν,I)≡E+ν1−I v(E) (2.3) for j=2,3,...,N. The boundary condition implies a relation between E1and I: I=(1+c)v(E1)⇔f(E1,ν,I)=E1−cν .(2.4) In the experimental situation, for a fixed ν, the current is measured as the applied voltage is swept. To solve our stationary problem it is easier to fix I(or equivalently E1), find the stationary electric-field profile {Ei},i=1,2,...,N, by forward iterating the map (2.2), and calculate the corresponding value of φ. Note that the constant voltage condition (1.5) forces the current for each electric field profile to be I=1 N N  i=1 1+Ej−Ej−1 νv(Ej).(2.5) Thus, Iis a good scalar measure of the solution of the system, and therefore the I−φ curve corresponds to a stationary bifurcation diagram. In Figure 2 we show, in the solid lines, three different bifurcation diagrams for ν=0.1, 0.175, and 0.5(N= 20 and c=10 −4fixed). The different structures of these diagrams are related to the properties of f(E,ν,I). This nonlinear map has been investigated in detail in references [3, 42] and we include here a brief discussion for the Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php BIFURCATION IN SUPERLATTICES 2035 Fig. 2.Stationary bifurcation diagram for three different values of ν. Figure (a) corresponds to a value of νfor which the unique stationary solution for each φis uniform (the stationary bifurcation diagram coincides with the drift velocity v(E)). Figure (b) displays a Z-shape form and presents oscillatory solutions for intermediate values of φ. The dotted lines are the maximum and minimum of the oscillatory current. Note in the inset that there are two windows of oscillations. Figure (c) exhibits several branches of domain solutions. In Figure (b) we have included the v(E)curve in thick dashed lines. sake of completeness. In the upper panel of Figure 3 we plot the map ffor I=0.6 and the same three values of the doping as in the previous figure. The stationary solution is constructed and depicted in the lower panel. Figure 3(a) corresponds to a uniform solution and Figures 3(b) and 3(c) to a nonuniform solution which presents two regions where the electric field is almost constant (the domains). We call it a domain solution. The map will not be invertible if its derivative, ∂f(E,ν,I) ∂E =1+νI v(E)2 dv(E) dE ,(2.6) changes sign on an interval of E. Due to the existence of the NDC (v(E) has a negative slope on a subinterval of E>1), this occurs for large enough ν. In the case of a noninvertible f(E,ν,I) (see Figure 3(c)), multiple stationary solutions are possible for the same value of I. Among these, there are domain solutions, and the stationary I−φcurve oscillates with φ(Figure 2(c)). For the lowest νvalue in Figure 2, the I−φcharacteristic curve coincides with the v(E) curve (dashed line), which means a uniform stationary solution for all the voltages φ. When the value of νis increased (ν=0.175 in Figure 2), the characteristic curve develops a Z-shaped form, indicating nonuniform solutions for a wide range of φ values. Furthermore there are time-periodic solutions for which we have depicted their maximum and the minimum value of the current (obtained by direct simulation of the system) as dotted lines in Figure 2. The oscillatory solution branches begin and end via supercritical Hopf bifurcations. In the inset we show that there are two voltage intervals for which time-dependent solutions exist. Finally, for still larger values of the doping, ν=0.5, the map f(E,ν,I) is noninvertible and the I−φcurve has 2N+1 branches corresponding to domain solutions. For a large enough I, their electric field profiles consist of two regions with Ejroughly constant and approximately equal to one of the three roots of the equation v(E)=I(the fixed points of f(E,ν,I)). Both domains are joined by a transition layer in which the field takes on intermediate values between the low field and the high field domains. We may define the domain boundary as the first well after the low field domain has ended, for which the field takes values on another branch of f(E,ν,I). If the electric field takes values only on branches of v(E) with positive slope, the stationary state is stable [42]. The voltage (1.5) depends on the position of the domain boundary, and this position characterizes each of the Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 2036 M. MOSCOSO, J. GAL´ AN, AND L. L. BONILLA 010 20 0.6 0.8 1 1.2 1.4 1.6 010 20 1 1.5 2 2.5 010 20 0.5 1 1.5 2 2.5 0.5 11.5 2 0.5 1 1.5 2 2.5 0.5 11.5 2 0.5 1 1.5 2 2.5 0.5 11.5 2 0.5 1 1.5 2 2.5 (a) Ej1 Ej j Ej (b) (c) Fig. 3.Construction of the stationary state iterating map fequation (2.2) for I=0.6,c= 10−4, and N=20. The upper panel is the map and the lower panel is the electric field profile as a function of the site for different values of ν. Column (a) corresponds to an uniform solution, column (b) to an unstable stationary solution in the oscillatory regime, and column (c) to a domain solution. solution branches in Figure 2(c). The field value jumps from the first to the third branches of v(E) for the N+1 stable solution branches of this diagram. Between two stable branches there is a branch of unstable stationary solutions which have at least one field value in the NDC region [42]. We have followed the solution branches bifurcating from the diverse stationary states by numerical continuation. The resulting diagrams for large νare not simple. In Figure 4 we show a blow-up of a peak in the I−φcharacteristic curve 2 for 1.1<φ<1.23 and ν=0.5, and we have marked the bifurcation points by arrows. The stable stationary solution (solid line) undergoes a Hopf bifurcation (HB) just before a limit point (LP) is reached, and the stationary solution becomes unstable (dashed line). After an additional limit point (not shown in the figure) the stationary solution becomes stable again and heads to the next peak. Direct simulation of the equations misses, as in the experimental situation, the unstable branch. We will explain this behavior in section 3 in a more complete way when analyzing the structure of the phase diagram. 2.2. Spatially chaotic stationary solutions. Note that the previous discussion holds for a given value of the parameter c. We have chosen cpositive and small as an approximation for the real contact. However, these results do not depend qualitatively on the value of c. A more detailed characterization of the contacts would involve a quantum mechanical microscopical calculation [1]. Let us relax this boundary condition and accept any value of E1−E0resulting from solving the discrete map (2.2) backward. Thus we will fix Iand ENand obtain directly EN−1,E N−2,...,E 1via (2.2) and the corresponding value of φ. The result of this construction is that the stationary problem has spatially chaotic solutions if Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php BIFURCATION IN SUPERLATTICES 2037 1.10 1.15 1.20 φ 0.98 0.99 1.00 1.01 I HB LP Fig. 4.Blow-up of the I−φcharacteristic curve for 1.1<φ<1.23 and ν=0.5. The stable stationary solution (solid line) undergoes a Hopf bifurcation (HB) just before a limit point (LP) is reached and becomes unstable (dashed line). The squares are from numerical simulation of the system and miss the unstable solutions. the number of quantum wells Nis large enough. This could have some relevance for the dynamics. Actually, if we increase the value of νwe reach a point where the derivative of the map (2.6) is equal to −1 and the map undergoes a period doubling bifurcation, and a period doubling cascade. Figure 5 shows the possible values of Ejfor each value of ν(at I=0.6 fixed). For a given ν, we eliminate the first 50 transient values and we depict the next 450 values of Ej. There are values of νfor which a period doubling bifurcation occurs (note where the single line splits into two branches). The period doubling occurs for ν=0.716, the period 4 for ν=0.906, and for ν>0.956 there are spatially chaotic stationary solutions. For the latter case, all values of Ejin a certain range are possible. In Figure 6 we have depicted the electric field profiles corresponding to the abovementioned cases for different values of ν. Before leaving this subject let us note that these solutions could be observed for very high doping (difficult to achieve in experiments), but a spatially chaotic electric field domain configuration would be, presently, of doubtful interest for a device. 3. Phase diagram. In the previous section we analyzed the construction of stationary solutions and found that there are three different qualitative regimes: •for small values of νwe have only stable uniform solutions; •for large values of νwe have multistability and domain solutions; •for intermediate values of νand appropriate values of φwe find self-sustained oscillations. In this section we will locate these regions more precisely and cast the results in aφ−νphase diagram for fixed values of cand N(c=10 −4,N= 20). We will use Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 2044 M. MOSCOSO, J. GAL´ AN, AND L. L. BONILLA 1.0 1.2 1.4 1.6 1.8 2.0 φ 0.10 0.15 0.20 0.25 ν TB DH Fig. 11.Total phase diagram of the model for N=20and c=10 −4. The dotted lines are curves of stationary saddle-nodes. For the sake of clarity, we have plotted only the main line of homoclinic orbits which sprout from the TB point (thin solid line). We have not shown other homoclinic orbits: There is one curve of homoclinic orbits for each Hopf curve as shown in Figure 10. are very close to each other and all tend to infinity as the doping νincreases. This suggests that there is a TB point at ν=∞, at which all these curves intersect tangentially. In a very narrow bias interval the stable equilibrium undergoes an HB followed by LP. The evolution of the leading eigenvalues determining the linear stability of the stationary state and the corresponding phase portrait is shown schematically in Figure 12. First, two complex conjugate eigenvalues cross the imaginary axis and the stationary solution becomes an unstable focus (HB). Very soon after the Hopf bifurcation the imaginary parts of the eigenvalues vanish and the stationary state changes to an improper node (IN). That means that only one unstable eigendirection exists; the linearization matrix corresponding to the leading eigenvalues is nondiagonalizable. At this point, one of the eigenvalues heads on toward the real axis toward the stable region and the other moves in the opposite direction. When the first eigenvalue crosses the imaginary axis we have a limit point (LP), but the stationary state remains unstable. The remaining unstable eigenvalue returns later to the stable region with another limit point bifurcation (not shown in Figure 4 but visible in Figure 2(c)) and continues to the next peak in the I−φbifurcation diagram. The periodic orbit that was born at the Hopf bifurcation collides with a curve of homoclinic orbits and disappears at a bias very close to HB and LP. 3.4. Comparison with experiments. As mentioned in section 2, our model explains in a simple way the shape of the stationary characteristic curve for different values of the doping. In particular, for high doping values there is a set of equidistant stable branches in the experimental characteristic curve, which are related to the existence of domain solutions. We have shown that the basic element is the presence of an NDC in the charge transport mechanism. All the proposed models in Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php BIFURCATION IN SUPERLATTICES 2045 IN Im Im Im Im Im Re Re Re Re Re λ λ λ λ λ λ λ λ λ λ LPHB Fig. 12.Scheme of the evolution of the leading eigenvalues of the equilibrium around the HB and LP in Figure 4. IN stands for improper node. the literature include this effect either from quantum mechanical calculations or from phenomenological considerations. On the other hand, our model has the potential of displaying the different dynamical behavior experimentally observed. Recent measurements show that the oscillatory behavior can start via a supercritical [21] or a subcritical [48] Hopf bifurcation. In the first case, the oscillation of the current has an amplitude increasing as (φ−φα)1/2, while the oscillation frequency is finite and different from zero at the instability onset. In the second case, hysteresis is observed in the time-averaged I-V characteristics. By sweeping up the voltage interval, the current exhibits a jump at a certain voltage value, while another jump occurs at a smaller voltage value when the direction of the sweeping is reversed. Note that the doping cannot be changed once the sample has been grown, and therefore all the experiments are carried out for a fixed ν. However, in our study we can sweep in νfor fixed φand even vary both parameters simultaneously. This provides us a complete understanding of the whole parameter space and the knowledge of the corresponding behavior for any values of νand φ. Moreover, from the computed phase diagram we predict the existence of voltage windows of self-sustained current oscillations for values of νaround ν=0.20. In fact, these windows have already been observed in photoexcited1SLs [29]. 4. The onset of oscillations: Amplitude equations. In the previous section we saw that the bifurcation on the main Hopf curve is supercritical below a certain value of νand subcritical above it. We computed this critical value numerically and got νDH =0.193 and φDH =1.168. Our aim in this section is to describe, for νfixed, the solutions of (1.5)–(2.1) perturbatively when the control parameter φexceeds a certain critical value φcwhere the stationary state, calculated in section 2, loses its stability. Close to this threshold, the time evolution is described as a product of a slowly varying amplitude and a pattern with faster time-dependence which comes from 1It is possible to change the carrier density by illuminating the sample with a laser with frequency larger than the semiconductor gap and creating electron-hole pairs. The laser power plays a similar role to the doping density, but a detailed study is needed to predict theoretically the bias windows in the photoexcited case. Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 2046 M. MOSCOSO, J. GAL´ AN, AND L. L. BONILLA 1.160 1.165 1.170 1.175 φ 20 10 0 10 20 Reγ Reλ 1 DH Fig. 13.Coefficients of the amplitude equation as a function of φalong the main Hopf curve. For Re(γ)<0the Hopf bifurcation is subcritical, whereas for Re(γ)>0it is supercritical. DH stands for degenerate Hopf bifurcation point. Note that Re(λ1)is always positive at the onset of oscillatory solution and represents the velocity at which the eigenvalues cross the imaginary axis. the linear analysis of the equations. This leads to a multiple time scales expansion, and to the determination of an equation of motion for the amplitude that describes slow modulations in time of the variables {Ei},i=1,2,..., and I[9, 4]. We will derive this amplitude equation for the oscillatory solution that bifurcates from the equilibrium and will precisely locate the position of the degenerate Hopf bifurcation point. This approach is equivalent to the reduction to the center manifold and the computation of the Hopf bifurcation equation [45, 28]. In our case, the latter would be a formidable task due to the dimension and structure of the equations. An additional advantage of our approach is that it provides an approximation for the periodic orbit that can be used to predict new bifurcations (local or global) and can be extended to higher orders in a systematic way. The details of the calculation are given in the appendices. We just rewrite the final result. The time evolution for the Ejand the current Ifor a voltage close to the threshold (φ=φc+ϕ2) shows two times scales: a slow modulation of the amplitude (2t) and a fast oscillation (t) associated with the eigenvalue involved in the Hopf bifurcation. The parameter ϕis +1 at the onset (φc=φα) of the oscillatory solution and −1at the end (φc=φβ). The ansatz in the calculation is, E∗ i(t;)=E∗ i,c +E(1) i(t, T)+2[E(2) i(t, T)+E∗ i,2]+3E(3) i(t, T)+O(4), I∗(t;)=I∗ c+I(1)(t, T)+2[I(2)(t, T )+I∗ 2]+3I(3)(t, T)+O(4), t=t, T =2t,(4.1) Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php BIFURCATION IN SUPERLATTICES 2047 1.176 1.180 1.184 φ 0.126 0.128 0.130 0.132 0.134 Frequency ω+ε 2 θ(φ) Fig. 14.Comparison of the oscillation frequency from the simulation (thick line) with the prediction of the amplitude equation (thin line) for ν=0.1819. The distance in voltage to the Hopf bifurcation is measured by 2. 0 100 200 300 t 1.06 1.08 1.10 1.12 E 1 |A(T)| θT T T Fig. 15.Time evolution of the electric field at site 1(E1(t)) close after the onset of the Hopf bifurcation (ν=0.1819,φ =1.18). The two time scales are clearly visible. The upper inset shows the evolution of the module of the amplitude and the lower inset is the phase θT. Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 2048 M. MOSCOSO, J. GAL´ AN, AND L. L. BONILLA with E(1) i(t, T)=A(T)eiwtξi+c.c. , i =1,...,N, I(1)(t, T)=A(T)eiwt +c.c.(4.2) In these expressions, c.c. means complex conjugate of the previous term, 2is the distance in voltage to the bifurcation, A(T) is the amplitude of the oscillation, wis the imaginary part of the eigenvalue that crosses the imaginary axis at φ=φc,ξiis the corresponding eigenvector, and the E∗ i,j are stationary solutions of the perturbed system (Appendix C). The amplitude equation obtained imposing a solvability condition is dA dT =ϕλ1A−γA|A|2,(4.3) where λ1and γare complex numbers given by (C.17) in Appendix C. Note that it is nothing but the complex form of the equation we get for the Hopf bifurcation with the normal form technique. Equation (4.3) describes the time evolution of the amplitude of the solution close to the bifurcation in the slow time scale. It has a periodic solution of the form A(T)=Re iϕθ(T−T0),T 0constant, R=ϕRe(λ1) Re(γ),θ=Im(λ1)−Im(γ)Re(λ1) Re(γ).(4.4) In Appendix C it is shown that ϕRe(λ1>0) at the onset and at the end of the oscillations. Therefore, the stability of the oscillatory solution depends only on the sign of Re(γ). Let φc=φα;ifRe(γ)>0atφ=φc, then the system bifurcates toward a stable periodic solution with amplitude R(the stationary solution is stable for φ<φ cand unstable if φ>φ c). This is a supercritical Hopf bifurcation. On the other hand, if Re(γ)<0atφ=φc, then the branch of oscillatory solutions bifurcates to φ<φ cand the Hopf bifurcation is subcritical. The case φc=φβis similar. The point where Re(γ) = 0 is the degenerate or generalized Hopf bifurcation point. The divergence in the expression of Rindicates that we should calculate higher order in the amplitude equation to analyze the oscillatory behavior, in a similar way as we would do in normal form theory. Following the previous scheme we have computed λ1and γfor the main Hopf curve of the phase diagram. In Figure 13 we plot Re(γ) and Re(λ1) as a function of φin the neighborhood of the degenerate Hopf point (which is highlighted by a small circle). The numerical prediction is in very good agreement with the crossing point [(φ, ν)=(1.1682,0.1928)]. In summary, the previous analysis tells us that along the main Hopf curve the Hopf bifurcation is supercritical for ν<0.1928, whereas it is subcritical for ν>0.1928. Also, it provides the first correction to the frequency of the oscillation as we depart from the bifurcation point and the curvature of the square root behavior of the amplitude as a function of the bias. In Figure 14 we show the comparison of the numerically computed frequency and the one given by w+2θas a function of the bias for a fixed value of ν. The interpretation of the two time scales is clear. If the bifurcation is supercritical, above the critical value, any perturbation will take the system from the unstable Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php BIFURCATION IN SUPERLATTICES 2049 equilibrium point to the limit cycle. The amplitude of the oscillation will slowly increase as the variable describes a spiral around the stationary solution. The final stable periodic solution is given by (4.4). This is depicted in Figure 15 for our system, where we plot the electric field at the first site [E1(t)] as a function of time for ν=0.1819 and φ=1.18. The insets are the time evolution of the modulus of the amplitude |A(T)|and its phase θT, obtained by solving (4.3). 5. Conclusions. In this work we have studied in detail the bifurcation behavior of a drift model for a semiconductor SL. Our results can be used for a better understanding of the experimental data. In fact, the information contained in the bifurcation set (Figure 11) can be used to design new samples and predict bifurcation scenarios [27]. We have shown that qualitatively different behavior can be reached depending on the values of φand ν, which are proportional to the external DC voltage applied and the doping density in the quantum wells. The former can be easily controlled in the laboratory, but the latter is fixed once the SL has been grown. Therefore, it is necessary to know a priori an appropriate value of the doping density depending on the purpose of the device. The possibility of building a tunable room temperature GHz oscillator is of particular technological interest. The stationary solutions are controlled by a nonlinear map and explain the qualitative behavior found in the experiments. The dynamical behavior presents oscillatory solutions for certain parameter ranges and the bifurcation to the oscillatory solution can be sub- or supercritical. The interplay between Hopf, homoclinics, and saddlenode bifurcations gives rise to a rich phase diagram whose details have been thoroughly investigated. The presence of domain solutions is crucial for the disappearance of the oscillatory behavior. The organizing centers for the long time dynamics in a broad range of parameters are multiple Takens–Bogdanov bifurcation points and degenerate Hopf bifurcation points. The system shows a robust planar behavior typical of a two-dimensional dynamical system; although the system is high dimensional and strongly nonlinear the bifurcations that take place are those that can be found in a planar system. For instance, we have not found any signature of ˇ Sil’nikov behavior for the curves of homoclinic orbits that emanate from the Takens–Bogdanov points [8]. The calculations have been restricted for computational reasons to a number of sites N= 20, but we have checked that our results do not change qualitatively as the number of sites is increased. We have also shown that the boundary condition imposed in the model has nontrivial consequences. In fact, the system may exhibit stationary spatially chaotic solutions for high values of ν, but, on the other hand, we have not found any signatures of temporal chaotic solutions. In conclusion, we have obtained a fairly complete understanding of the dynamics of the model and predicted new interesting dynamical phenomena that could be relevant for the design of new semiconductor devices. Possible homoclinic bifurcations and the role in the dynamics of the spatially chaotic stationary solutions found in section 2.2 are still under study. Appendix A. Linearization. We linearize (2.1)–(1.6) around a stationary solution, {E∗ i},i=1,...,N,atφ=φc. Substituting Ei(t)=E∗ i+ε˜ei(t),I(t)=I∗+ε˜ j(t),0<ε1,(A.1) in (2.1), we obtain d˜ei dt +I∗v(E∗ i) v(E∗ i)˜ei+˜ei−˜ei−1 νv(E∗ i)−˜ j=0, Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 2050 M. MOSCOSO, J. GAL´ AN, AND L. L. BONILLA ˜e1−˜e0=0; N  i=1 ˜ei=0.(A.2) This system can be solved by separation of variables ˜ j=ˆ je λt ,˜ei=ˆeieλt .(A.3) Insertion of (A.3) into (A.2) yields [λ+b(0) i+a(0) i]ˆei−a(0) iˆei−1−ˆ j=0, ˆe1−ˆe0=0; N  i=1 ˆei=0,(A.4) where we have defined the following quantities that depend only on the stationary electric field profile: b(k) i=I∗v(k)(E∗ i) k!v(E∗ i),a (k) i=v(k)(E∗ i) ν,k=0,1,2,... .(A.5) In these expressions v(k)is the order kderivative of the nonlinear function v(E). In order to compute λas an eigenvalue of a matrix, we have to eliminate ˆ jand one of the ˆei. Taking ˆeN=N−1 i=1 ˆeiand summing (A.4) over iwe get ˆ j=1 N(b(1) 1−a(0) 2−b(1) N−a(0) N)ˆe1 (A.6) + N−1  k=2 (b(1) k+a(0) k−a(0) k+1 −b(1) N−a(0) N)ˆek and the eigenvalue equation Aˆe=−λˆe,(A.7) where Ais an (N−1)×(N−1) real matrix and ˆeisa(N−1)-dimensional eigenvector associated to the eigenvalue +λof (A.4). Defining c1=b(1) 1−a(0) 2−b(1) N−a(0) Nand ck=b(1) k+a(0) k−a(0) k+1 −b(1) N(2 ≤k≤N−1) for the coefficients in expression (A.7), Ahas the following form: A=−1 N     c1c2... c N−1 c1c2... c N−1 . . .. . .. . . c1c2... c N−1      (A.8) +         b(1) 100... 0 −a(0) 2b(1) 2+a(0) 20... 0 0−a(0) 3b(1) 3+a(0) 30... 0 . . .. . ........ . . 00 ... −a(0) N−1b(1) N−1+a(0) N−1.         . Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php BIFURCATION IN SUPERLATTICES 2051 Appendix B. Solvability condition. As we are interested in studying the onset and end of oscillatory solutions let us consider the case in which a pair of eigenvalues λ=±iw crosses the imaginary axis for a critical value of the voltage φc. In the perturbative expansion we will get a hierarchy of linear equations and we must impose a solvability condition to avoid secular terms. These would give rise to terms proportional to teiwt which become unbounded for t→∞and we are looking for a solution of the nonhomogeneous linear problem bounded and periodic in time. The nonlinear homogeneous problem is d˜ei dt +[b(1) i+a(0) i]˜ei−a(0) i˜ei−1=˜ h+fieiwt , ˜e1−˜e0=ge iwt ; N  i=1 ˜ei=0(B.1) for i=1,...,N. Using separation of variables as in (A.3) we get [iw +b(1) i+a(0) i]ˆei−a(0) iˆei−1=ˆ h+fi, ˆe0=ˆe1−g; N  i=1 ˆei=0.(B.2) This can be expressed in matrix form as (M+iwI))e =) h+) f+)g,(B.3) N  i=1 ˆei=0,(B.4) where ) h=( ˆ h, ˆ h,...,ˆ h)t,) f=(f1,f 2,...,f N)t,)g =(−a(0) 1g,0,...,0)t,Iis the N×N identity matrix and equals           b(1) 100... 0 −a(0) 2b(1) 2+a(0) 20... 0 0−a(0) 3b(1) 3+a(0) 3... 0 00 −a(0) 4b(1) 4+a(0) 40 . . .. . ........ . . 00 ... −a(0) Nb(1) N+a(0) N           .(B.5) The determinant of the matrix (M+iwI) is the product of the diagonal elements, and therefore it is always invertible because Mis real and λis complex by hypothesis in a Hopf bifurcation point. Then, the solution of (B.3) is )e =(M+iwI)−1() h+) f+)g), which has to fulfill (B.4). Let us first consider the term proportional to ) h. If we solve (A.4) for λ=iω in matrix form we get )e =[M+iwI]−1) j(B.6) with ) j=ˆ j(1,1,...,1)t. Taking into account that N i=1 ˆei= 0 we conclude that the sum of all the elements of the matrix [M+iwI]−1must vanish, i.e., N  j N  k [M+iwI]−1jk =0.(B.7) Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php 2052 M. MOSCOSO, J. GAL´ AN, AND L. L. BONILLA We can use expression (B.7) to show that the term proportional to ) h=ˆ h(1,1,...,1)t already fulfills the bias condition. Therefore, to avoid secular terms in the solution of (B.1), fiand gshould satisfy the following solvability condition: N  i=1 N  j=1 (M+iwI)−1ij(fj+gj)=0.(B.8) We will use (B.8) to get the amplitude equation in the next section. Appendix C. Derivation of the amplitude equation. We want to construct a periodic solution that bifurcates from φ=φcby means of a multiple scale expansion. Let 1 be a small parameter that measures the separation from the critical voltage φc: φ=φc+2ϕ,(C.1) where ϕ=+1orϕ=−1ifφccorresponds to the value of the voltage where the oscillatory branch begins or ends, respectively. The stationary electric field profile and the current for this value of φis E∗ i()=E∗ i,c +2ϕE∗ i,2+O(4),i=1,...,N, I∗()=I∗ c+2ϕI∗ 2+O(4).(C.2) The terms of order O(1) are known, and those of order O(2) satisfy the following equations: [b(1) i+a(0) i]E∗ i,2−a(0) iE∗ i−1,2−I∗ 2=0, E∗ 1,2−E∗ 0,2=0; 1 N N  i=1 E∗ i,2=1,(C.3) which can be solved numerically. To compute the solution that bifurcates at φcwe assume the usual Hopf multiscale ansatz [23, 4]. E∗ i(t;)=E∗ i,c +E(1) i(t, T)+2[E(2) i(t, T)+ϕE∗ i,2]+3E(3) i(t, T)+O(4), I∗(t;)=I∗ c+I(1)(t, T)+2[I(2)(t, T)+ϕI∗ 2]+3I(3)(t, T)+O(4), t=t, T =2t,(C.4) where i=1,...,N. The slow scale has been chosen so that the secular terms that first appear at order O(3) can be removed. In reference [19] it is shown that for the Hopf bifurcation the power series of the bifurcation parameter (φ) and the slow time scale (T) has to be even in . Substituting (C.4) in (2.1)–(1.6) and comparing powers of we get the following hierarchy of equations: O(): ∂E(1) i ∂t +[b(1) i+a(0) i]E(1) i−a(0) iE(1) i−1−I(1) =0, E(1) 1−E(1) 0=0; N  i=1 E(1) i=0.(C.5) Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php BIFURCATION IN SUPERLATTICES 2053 O(2): ∂E(2) i ∂t +[b(1) i+a(0) i]E(2) i−a(0) iE(2) i−1−I(2) =−[b(2) i+a(1) i](E(1) i)2+a(1) iE(1) iE(1) i−1, E(2) 1−E(2) 0=0; N  i=1 E(2) i=0.(C.6) O(3): ∂E(3) i ∂t +b(1) i+a(0) iE(3) i−a(0) iE(3) i−1−I(3) =−∂E(1) i ∂T −ϕ(2 b(2) i+a(1) iE∗ i,2−a(1) iE∗ i−1,2)E(1) i+ϕa(1) iE∗ i,2E(1) i−1−2b(1) i+a(1) iE(1) iE(2) i +a(1) iE(1) i−1E(2) i+a(1) iE(1) iE(2) i−1−b(3) i+a(2) i 2(E(1) i)3+a(2) i 2(E(1) i)2E(1) i−1, E(3) 1−E(3) 0=0; N  i=1 E(3) i=0.(C.7) We will solve these systems and apply the solvability condition (B.8) when we find a term proportional to e±iwt. We make the following ansatz for the form of the order O() time dependent electric fields and current, E(1) i(t, T)=A(T)eiwtξi+c.c., i=1,...,N, I(1)(t, T)=A(T)eiwt +c.c.,(C.8) where c.c. means complex conjugate of the previous term, and ξiis the eigenvector of (A.4) with λ=iw, i.e., it is the eigenvector corresponding to the eigenvalue which crosses the imaginary axis at φ=φc. This form is equivalent to a projection into the tangent space along the bifurcating eigendirection. Note that, at this level, A(T)is still undetermined. The eigenvector ξihas to be properly normalized so that the first equation in expression (A.4) can be satisfied for i= 1. This gives the condition ξ1(iw +b(1) 1)=1.(C.9) The equations to order O(2) are (C.10) ∂E(2) i ∂t +[b(1) i+a(0) i]E(2) i−a(0) iE(2) i−1−I(2) ={−2[b(2) i+a(1) i]|ξi|2+2a(1) iRe(ξi¯ ξi−1)}|A|2 +{−[b(2) i+a(1) i]ξ2 i+a(1) iξiξi−1}A2ei2wt +c.c., E(2) 1−E(2) 0=0; N  i=1 E(2) i=0. There are no resonant terms at this order either; we seek solutions of the form E(1) i(t, T)=α0,i|A(T)|2+α2,iA(T)2ei2wt +c.c., i=1,...,N, I(2)(t, T)=β0|A(T)|2+β2A(T)2ei2wt +c.c..(C.11) Downloaded 04/17/17 to 150.214.182.208. Redistribution subject to SIAM license or copyright; see http://www.siam.org/journals/ojsa.php