University of Southern California Dissertations and Theses

Search, synthesis, and analysis for semiconductor device design

(USC Thesis Other)

Search, synthesis, and analysis for semiconductor device design

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content SEARCH, SYNTHESIS, AND ANALYSIS FOR SEMICONDUCTOR DEVICE DESIGN by Jason Thalken _____________________________________________________________ A Dissertation Presented to the FACULTY OF THE GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulfillment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (PHYSICS) December 2006 Copyright 2006 Jason Thalken ii Acknowledgements Thanks must be given first and foremost to my PhD advisor, Stephan Haas, who has not only provided me with sound academic guidance, but a great friendship as well. He made sure I never forgot to enjoy my work and my life, which can sometimes be tough for a grad student wrapped up in his research. I am also thankful for the opportunities and steady guidance given to me by Professor Anthony Levi, who taught me the importance of presentation in science. Professors Gene Bickers, Hans Bozler, and Daniel Lidar all deserve thanks as well for serving on my dissertation committee and for being so accommodating and supportive despite their own busy schedules. I am always grateful for the constant support I have received from my family. My parents, my brother and my sister have all been unconditionally encouraging, which has enabled me to spend my education in pursuit of knowledge rather than money. Finally, I cannot ever forget the care and thoughtfulness of my lovely wife, Cathy, whose companionship makes even the toughest of times enjoyable, and who could very well be right when she claims that the support she has given me has earned her at least 50% of my degree. iii Contents Acknowledgments ii List of Tables v List of Figures vi Abstract x 1 Introduction 1 2 Optimization and the Inverse Problem 4 2.1 The inverse problem and adaptive quantum design 4 2.2 Example 1: Adaptive quantum design of atomic clusters 6 2.2.1 Physical model: Long-range tight-binding model 7 2.2.2 Example target functions and convergence criterion 9 2.2.3 Design of tight-binding clusters in continuous configuration space 12 2.2.4 Adaptive design in discrete configuration space 18 2.2.5 Summary and conclusions for example 1 24 2.2.6 Appendix to example 1: Optimization and search algorithms 26 2.3 Example 2: Adaptive design of excitonic absorption in broken-symmetry quantum wells 30 2.3.1 Physical model: Excitonic absorption in Al x Ga 1-x As quantum wells 31 2.3.2 Frequency switch: delocalized electron solution 36 2.3.3 Frequency switch: ionized hole solution 38 2.4 Search: Motivations, terminology, and terrain 40 2.4.1 Adaptive quantum design search conditions 41 2.4.2 Parameters, spaces, and mapping 41 2.4.3 Assessing fitness 43 2.4.4 Dimensionality reduction 46 2.4.5 Set-valued target functionals 47 2.5 Search: Methodology and implementation 49 2.5.1 Downhill methods 49 2.5.2 Simplex methods 50 2.5.3 Gradient methods 50 2.5.4 Annealing methods 52 2.5.5 Genetic algorithm methods 53 2.5.5.1 Encoding: Chromosomes 54 2.5.5.2 Selection schemes 55 2.5.5.3 Crossover 60 2.5.5.4 Mutation 62 2.5.6 Hybrid techniques 63 2.5.7 Reset techniques 64 iv 3 Optimization without a target function 66 3.1 Example 3: Quantum well electroabsorption modulator 66 3.2 Motivations and methodology 68 3.2.1 Defining utility 70 3.2.2 Classifying solutions 72 3.3 Example 4: Double well electroabsorption modulator 77 3.4 Example 5: Initially separated double well solution 79 4 Experimental efforts and analysis 83 4.1 Incorporating fabrication and pragmatics into the model 83 4.1.1 Poisson solver 84 4.1.2 Sensitivity analysis for inaccuracies in well thickness 86 4.2 Experimental results 89 4.2.1 1550 nm MOCVD Growth InGaAsP integrated devices 89 4.2.2 850 nm MOCVD Growth GaAs device 96 4.2.3 840 nm MBE Growth GaAs devices 99 5 Summary and conclusions 102 References 107 Bibliography 111 Appendix: Quantum well device synthesis program guide 114 A-1 Modeling the materials 115 A-2 Setting up the quantum well structure 121 A-3 Calculating the absorption 126 A-4 Data presentation and analysis 130 A-5 The Poisson solver 138 v List of Tables Table A1: 122 Global variables defined before absorption is calculated Table A2: 128 Global variables generated during absorption calculation Table A3: 129 Variables, default values, and explanations for absorption.m Table A4: 140 Global variables created by read_poisson.m Table A5: 142 Global variables created by poisson_solver.m Table A6: 143 Variables, default values, and explanations for poisson_solver.m vi List of Figures Figure 1.1: 2 Flow chart outlining the steps taken in the iterative design processes addressed in this work. Figure 2.1: 8 Density of states in spatially invariant tight-binding systems with periodic boundary conditions. Figure 2.2: 10 Target densities of states. Figure 2.3: 13 Adaptive quantum design applied to 16 atoms in a two-dimensional box with periodic boundary conditions. Figure 2.4: 14 Same as Fig. 2.3, but for a symmetric two-peak target function. Figure 2.5: 15 Spectral response of tight-binding dimers, trimers, and quadrumers. Figure 2.6: 17 Asymmetric two-peak quasiparticle density of states. Figure 2.7: 19 Adaptive design of clusters with 48 atoms on a discrete lattice. Figure 2.8: 21 Effect of the coarseness of the underlying lattice on the formation of building blocks. Figure 2.9: 24 Effect of the density of atoms and the power-law dependence of the tight- binding overlap integral the convergence of adaptive quantum design. Figure 2.10: 27 Comparison of the convergence of several optimization algorithms. Figure 2.11: 32 Changes in excitonic absorption spectrum for a symmetric quantum well with increasing electric field. vii Figure 2.12: 35 Exciton binding energy calculated using the tight binding Hamiltonian as a function of well width for Al 0.3 Ga 0.7 As/GaAs quantum well structures. Figure 2.13: 37 Broken-symmetry double quantum well obtained from numerical optimization of well width and depth parameters. Figure 2.14: 39 Broken-symmetry quantum well structure obtained from numerical optimization. Figure 2.15: 48 An illustration of a set-valued target functional applied to a specific solution trajectory. Figure 3.1: 67 Broken-symmetry quantum well structure obtained by modifying the structure from Fig. 2.14 to fit a new target functionality. Figure 3.2: 71 Band-edge potential profiles and the corresponding absorption peak paths selected by a machine-based search of configuration space without a human- specified target function. Figure 3.3: 73 Various categories for exciton absorption peak paths exhaustively sampled for the case of quantum well structures generated from two equipotential regions (N=2). Figure 3.4: 75 The paths of the electron and hole wave function peaks as a function of applied field, F, for a square well 26 monolayers wide. Figure 3.5: 77 The paths of the electron and hole wave function peaks as a function of applied field, F, for a step well 52 monolayers wide, with a 40 monolayer step of 15% Al. Figure 3.6: 78 Asymmetric double quantum well intensity modulator discovered by machine-based exploration without human specification of a target function. Figure 3.7: 81 Electron and hole probability densities for the ISDW solution. viii Figure 3.8: 82 Absorption spectrum for the ISDW solution. Figure 4.1: 86 The conduction and valence band edge profiles for a 6 quantum well AlGaAs/GaAs electroabsorption modulator grown in a PIN layer structure. Figure 4.2: 87 Absorption spectra plotted at F = 0 kV/cm for conventional Al 0.3 Ga 0.7 As/GaAs structures of varying well widths. Figure 4.3: 88 Insensitivity of ISDW structure to monolayer fluctuations, despite narrow well width. Figure 4.4: 90 Measured and simulated absorption coefficients at 1550 nm for a conventional square well device with integrated laser. Figure 4.5: 91 Conduction and valence band profiles for InGaAsP Step-Well 1. Figure 4.6: 92 Predicted absorption as a function of applied field, F, at 1530 nm for InGaAsP Step-Well 1. Figure 4.7: 93 Measured absorption as a function of V bias at 1530 nm for InGaAsP Step- Well 1. Figure 4.8: 95 Conduction and valence band profiles for InGaAsP ISDW structure. Figure 4.9: 97 Predicted absorption spectrum for AlGaAs/GaAs step well structure. Figure 4.10: 98 Measured absorption spectrum for AlGaAs/GaAs step well structure. Figure 4.11: 99 Predicted absorption spectrum for conventional Al 0.2 Ga 0.8 As/GaAs square well structure. ix Figure 4.12: 100 Measured absorption spectrum for conventional Al 0.2 Ga 0.8 As/GaAs square well structure. Figure A1: 133 Output from the function graph_structure.m after data text has been moved. Figure A2: 134 Sample figure generated by graph_wavefunctions.m using a conventional GaAs quantum well structure. Figure A3: 135 Absorption spectrum for a conventional well structure. Figure A4: 136 Absorption at 840 nm as a function of applied electric field for a conventional quantum well structure. Figure A5: 144 The potential profile for the conduction band at each applied voltage step. Figure A6: 145 The potential profile for the conduction and valence bands with zero applied voltage. Figure A7: 146 The electric field profile for the structure at each applied voltage step. Figure A8: 147 The electric at zero applied volts with the well structure superimposed. Figure A9: 148 The charge density as a function of position for the structure at each applied voltage step. x Abstract Synthesis of semiconductor device design requires access to realistic physical models and adaptive algorithms. To demonstrate that such synthesis is feasible, atomic clusters with specified quasiparticle densities of states are designed using advanced search algorithms, and broken-symmetry quantum-well potential profiles are developed with optical response properties superior to previous ad hoc solutions. Advanced stochastic search algorithms are introduced in the context of a guide to efficient implementation of inverse problem solving methods, and then, it is shown that automated searches of configuration space may be performed without the need to input a specific target function. Improved electroabsorption modulator designs resulting from these machine-based searches are given as a demonstration of the utility of targetless exploration. Finally, the results of experimental testing of the asymmetric quantum well structures discovered by machine are presented, and it is established that in order to fabricate functional devices from these designs, significant improvements in the growth and measurement processes are required. 1 Chapter 1 Introduction The design of both photonic and electronic devices has evolved over many years in an essentially ad-hoc way. There is, for example, a long history of laser diode and photodiode design that incrementally builds on previous designs and concepts. Now an opportunity exists to harness a combination of modern compute power, adaptive algorithms, and realistic physical models, to seek robust, manufacturable designs for components that meet previously unobtainable system specifications. In such a scheme, device design is synthesized. In its simplest form, semiconductor device synthesis software solves the inverse problem by identifying the best configuration to achieve a user-defined specification or target function. Typically, this involves searching for a (broken-symmetry) spatial configuration of a semiconductor structure that produces a desired target function response. Using such an approach it is possible to break the cycle of piecemeal design methods and discover completely new functionalities for both photonic and electronic semiconductor devices. The purpose of this dissertation is not only to demonstrate 2 that such device synthesis is, in fact, feasible, but also to serve as a guide for researchers attempting to implement these methodologies in new physical systems. The design processes developed and implemented in this dissertation involve an iterative scheme where a desired objective is defined, and input parameters to a physical model are adjusted to find the best fit to the objective. This scheme is outlined in Fig. 1.1, and will be elaborated upon throughout this dissertation. Figure 1.1: Flow chart outlining the steps taken in the iterative design processes addressed in this work. First an objective is defined, and an initial guess for physical input parameters is fed into the physical model. The results of the physical model are compared to the target using a fitness value, and if the objective is not yet met, the input parameters to the physical model are changed and the evaluation process begins again. Chapter 2 opens with an introduction to the inverse problem and the process of adaptive quantum design through the example of arranging tight-binding atoms on a Input objective function Input initial parameter values No Output optimized result Forward solve: Physical model based Yes Evaluation of fitness Change parameter values Objective met? 3 1 or 2 dimensional surface. The more complex example of excitonic absorption in broken-symmetry quantum wells is then presented to serve as a further example of how the design process is implemented. The process of searching for an acceptable solution is then introduced, and various search algorithms are presented in the form of a guide to choosing the best method and parameters for a new design investigation. Chapter 3 focuses on the idea of performing a search for any possible interesting or useful device, without a human-specified target, and Chapter 4 presents the accounts and results of the efforts made to experimentally verify the asymmetric quantum well solutions discovered during the search and classification in chapters 2 and 3. The appendix is a reference guide to the suite of matlab functions developed for the design and evaluation of solutions to be grown and experimentally tested. It includes functions specifying material and experimental parameters, the forward problem solver, a Poisson solver incorporating the entire PIN structure encasing the quantum well of the forward problem in the I-region, and several other functions developed for the purpose of visualizing an understanding a particular solution. It is not meant to be read in a single sitting, but rather to be referred to while using the suite to develop devices for fabrication. 4 Chapter 2 Optimization and the Inverse Problem 2.1) The inverse problem and adaptive quantum design In the near future, it will likely be possible to control the precise spatial positions of atoms and molecules using the experimental techniques now being developed by nanoscience [1, 2, 3, 4]. To complement these emerging capabilities it is clear that a new set of theoretical tools has to be developed to assist in the exploration of a potentially vast number of atom configurations and a corresponding enormous range of physical properties. In this section, I will outline an approach, called “adaptive quantum design” [5], which sets out to address this challenging task. In contrast to classical systems, atomic scale devices exhibit quantum fluctuations and collective quantum phenomena caused by particle interactions. Besides offering an excellent testing ground for models of correlated electrons, they also force us to reconsider conventional paradigms of condensed matter physics, such as crystal symmetries that are imposed by nature. In some instances, such symmetries need to be explicitly broken in order to enable or optimize a desired system response. Consider for 5 instance the quasiparticle density of states in tight-binding systems, which is to be the subject of the first example. In translationally invariant structures, i.e. crystals, it is well known that the spectral response function exhibits van-Hove singularities at positions of low dispersion, such as the band edges in a one-dimensional chain or the band center in a two-dimensional square lattice. These enhancements of the density of states can be very useful in amplifying system responses such as optical conductivity at specific incident energies. It is therefore important to be able to control the positions and shapes of such features by adaptive design techniques applied to models which capture the essential degrees of freedom of interacting atomic clusters. Traditional ad-hoc methods for the design of nanoscale devices will likely miss many possible configurations. At the same time, it is unrealistic to expect individuals to manually explore the vast phase space of possibilities for a particular device function. The proposed solution to this difficult design problem is to employ machine-based searches of configuration space that enable user-defined target functions. Adaptive quantum design solves the inverse problem by numerically identifying the best broken-symmetry spatial configuration of atoms and molecules that produce a desired target function response. The two major ingredients of adaptive quantum design are the physical model, which in this example evaluates the electronic density of states for a particular spatial arrangement of the atoms, and the search algorithm that finds the global minimum in 6 the parameter space of all possible configurations. This problem is typically highly underdetermined in the sense that there can be several atomic configurations that yield a system response very close to the desired target function. Often, the associated landscape of solutions is shallow and has many nearly degenerate local minima. Adaptive quantum design therefore relies heavily on efficient algorithms that accurately model the interacting system and find its optimal configuration, i.e. the global minimum in the available parameter space that yields the best match to the desired target function. This target response may be a specific angular transmission of a photonic crystal, a fine-tuned Josephson current along a superconducting junction array, or an energy-dependent density of states profile of an atomic cluster - the case on which we focus here. In this example, various numerical search tools are applied, including guided random walk, simulated and triggered annealing, and genetic algorithms. In particular the last set of techniques is most efficiently implemented on parallel computers. It is observed that in practice hybrids of these methods yield the best results, i.e. fast convergence towards a global minimum. 2.2) Example 1: Adaptive quantum design of atomic clusters The particular example of a long-range tight-binding model is chosen for this study because it captures essential features of correlated quantum mechanical systems, and yet permits fast numerical diagonalization of relatively large clusters with broken translational symmetry. In a typical optimization run, these “function calls” occur 7 100-10,000 times. This model is therefore suitable for developing and testing adaptive design algorithms, and should be viewed as an initial step towards the design of atomic clusters. The purpose of this example is to show the applicability and limits of adaptive design techniques on a simple but non-trivial quantum system. 2.2.1) Physical model: Long-range tight-binding model The tight-binding approach is an effective tool to describe the band structure of electronic systems. It is commonly used to model the relevant bands close to the Fermi level, obtained from complex density functional theory calculations. Since in this work symmetry-breaking, non-periodic configurations are considered, a long- range variant of the tight-binding model with overlap integrals depending on the variable inter-atomic distance has to be used. Its Hamiltonian is given by ( ) ∑ + − = j i j i j i j i c c c c t H , † † , , (2.1) where c i † and c i are creation and annihilation operators at a site r i , and the sum over pairs of atoms is restricted to avoid double-counting. Here the spatial decay of the overlap integral t i,j is parameterized by a power-law, α j i j i t t r r − = , , (2.2) where the exponent α is taken to be 3.0 throughout this example unless mentioned otherwise [6]. This parameterization reflects an algebraic variation of the overlap integral with inter-atomic separation. The choice of sign in the Hamiltonian follows the convention for s-orbitals. However, this simple implementation of the tight- 8 binding model does not account for the orbital directionality of realistic Au, Ag, or Pd atomic wave functions. More sophisticated and numerically expensive techniques, such as the local density approximation, would be needed to make quantitative predictions for these systems. Figure 2.1: Density of states in spatially invariant tight-binding systems with periodic boundary conditions. (a)Nearest-neighbor chain (no long-range overlap integral) with 30 atoms. (b)Same as (a), but with long-range overlaps according to Eq. 2.2. (c) Nearest- neighbor square lattice with 400 atoms. (d)Same as (c), but with long-range overlaps. The Hamiltonian matrix of the long-range tight-binding model in the basis of single- particle states is non-sparse, and only its diagonal matrix elements vanish. In order to obtain the spectrum, the matrix is diagonalized numerically for finite clusters. In Fig. 9 2.1, the resulting densities of states of translationally invariant chains and square lattices are shown for the nearest-neighbor tight-binding model with t i,j = tδ i,j in Figs. 2.1(a) and (c), and for the case of long-range overlap integrals in Figs. 2.1(b) and (d). Characteristic van-Hove singularities are observed, in one dimension at the band edges, and in two dimensions at the band center. For the long-range model, the particle-hole symmetry is broken because of frustration introduced by the competing overlaps, leading to asymmetries in the density of states. While the system sizes in Fig. 2.1 are chosen to be rather large in order to make contact with the familiar thermodynamic limit, there are still some visible finite-size remnants, i.e. a faint pole structure due to the discreteness of the system. These features become much more pronounced for the few-atom clusters that are studied in this example. 2.2.2) Example target functions and convergence criterion While “nature” gives us densities of states that are constrained by the dimensionality and symmetry of the underlying lattice, our objective is to engineer spectral responses with specific shapes that are useful in designing nanoscale devices. For example, we may wish to produce a quasi-two-dimensional spectrum in a one- dimensional system or to concentrate spectral weight in particular energy windows. These goals are achieved by placing the atomic constituents into optimized symmetry-breaking configurations which are determined by numerical searches. 10 Figure 2.2: Target densities of states used in this work: (a) top hat function, centered at E = 0, (b) particle-hole symmetric two-peak function, and (c) asymmetric two-peak function. For systems with finite numbers of particles these shapes are approximated by quasiparticle peaks. The specific target functions that are studied in this work are shown in Fig. 2.2. They are (a) a flat top hat density of states centered at E = 0, N(E) target = θ(E - E c )θ(E + E c ), (2.3) where θ(x) is the Heavyside function, and E c is an energy cut-off, (b) a symmetric two-peak function, centered at E = 0, N(E) target = θ(E - E c2 )θ(E + E c1 ) + θ(E - E c1 )θ(E + E c2 ), (2.4) 11 and (c) a particle-hole-symmetry-breaking function with two unequal peaks, i.e. more spectral weight on the quasi-hole side (E < 0) than on the quasi-electron side (E > 0) of the spectrum, N(E) target = θ(E - E c2 )θ(E + E c1 ) + θ(E - E c4 )θ(E + E c3 ). (2.5) In systems with finite numbers of tight-binding atoms, these continuous shapes are approximated by equally spaced poles within the energy windows where N(E) target is a non-vanishing constant. Here, the delta-functions are given a finite width of 0.02t. As more atoms are added to the system, these peaks merge together, approaching the bulk result. For other non-flat target functions the quasiparticle peak spacing can be varied, e.g. following a Gaussian or Lorentzian shape. Naturally, not all targets can be achieved equally well. Factors that influence the achievable match to a target include the number of available atoms and, as will be shown, a continuous or discrete number of accessible spatial positions. The dimensionality of the system poses additional constraints. In particular, there are less available configurations in lower dimensions. This study focuses on three prototype spectral responses which are targeted by numerically optimizing configurations of clusters with up to 48 atoms in two spatial dimensions. The optimization algorithms seek to minimize the deviation from a given target density of states, defined by the error function 12 ( ) ( ) [ ] ∫ ∞ ∞ − − = Δ 2 target E N E N dE (2.6) which is the least-square difference between the system response for a given configuration and the target response function. We have explored a number of numerical techniques, including the Newton-Raphson steepest descent method, simple downhill random walk, simulated and triggered annealing, and genetic algorithms. The advantages and disadvantages of these techniques are discussed briefly in section 2.2.6, and will be elaborated upon throughout this work. In general, it is found that hybrids of these methods tend to work best. In the next section, we focus on adaptive design of atomic clusters in continuous configuration space without any restrictions to underlying discrete lattices. In this case, the search space is infinite which generally allows better convergence to a given target response than in finite configuration space. However, some experimental realizations of such structures require deposition of atoms on substrates with discrete lattice structures [7, 8, 9, 10, 11] Therefore, this case is addressed separately in the subsequent section. 2.2.3) Design of tight-binding clusters in continuous configuration space The first target density of states we would like to study in detail is the particle-hole symmetric top hat function, centered at energy E = 0, shown in Fig. 2.2(a). This function represents a constant density of states for a bulk solid between the energy cut-offs ±E c , giving a bandwidth W = 2E c . Here, we choose E c = 3t. However, it should be noted that with the adaptive quantum design approach we are not restricted 13 to this choice, and target densities of states with quite different bandwidths can be matched, although often to a lesser degree of accuracy. Figure 2.3: Adaptive quantum design applied to 16 atoms in a two-dimensional box with periodic boundary conditions. The target density of states is a symmetric top hat function with bandwidth 6t. The atomic configurations are optimized by applying a guided random walk algorithm to the long-range tight-binding model. (a) The best matching solution. (b) Contour plot of the potential of the resulting spatial configuration. (c) Convergence to the target function, log 10 (Δ), with the number of updates. In Fig. 2.3, the solution for a system with 16 atoms confined to a box with periodic boundary conditions is shown. The guided random walk method is applied to optimize the configuration of atoms by iterative local updates of their positions in order to match the top hat density of states. As observed in Figs. 2.3(a) and (c), good convergence to the target function can be achieved for this case after less than 2000 updates. A contour plot of the potential Σ ij - t ij / |r i - r j | α for the resulting spatial configuration is shown in Fig. 2.3(b). Here, equipotential lines are used to denote the positions and overlapping wave functions of the atomic constituents in the system. For this target function one discovers the formation of dimers with a wide range of 14 inter- and intra-dimer spacings. These self-assembled building blocks have variable directional orientation, and are closely packed. Figure 2.4: Same as Fig. 2.3, but for a symmetric two-peak target function with peaks of bandwidth 2t centered at -3t and 3t. To explore the generic aspects of this result, let us now turn to the two-peak target function shown in Fig. 2.2(b). This density of states has a gap in the center of the spectrum, as may be desired for the construction of two-state systems such as q-bits. As shown in Figs. 2.4(a) and (c), the convergence to this target function is not quite as good as for the top hat function, but saturation occurs already at half the number of iterations compared to the previous example. Interestingly, the optimized spatial configuration that is found in Fig. 2.4(b) also displays a preference for dimer formation. 15 Figure 2.5: Spectral response of tight-binding dimers, trimers, and quadrumers. (a) Potential contour plot of dimer, trimer, and quandrumer of size L 0 . (b) In the dimer, the single-atom level (dashed line) is split symmetrically (solid line). Competing interactions (frustration) in the trimer and quadrumer lead to asymmetric densities of states. The strength of the frustration depends on the separation L. (c) Spectral peak positions as a function of normalized distance L/L 0 for dimer, trimer, and quadrumer. The spectra for L/L 0 = 1 are shown as solid lines and L/L 0 = 3 as dashed lines in (b). In the long-range interacting systems we are studying it is, at least at first sight, not obvious why these target functions prefer dimer building blocks. Let us address this question by examining the individual spectra of the molecular building blocks shown in Fig 2.5. Each dimer molecule contributes to the density of states a positive and negative pole with energies E = ±t 12 , where t 12 is the hopping integral between the two participating atoms. The zero-energy quasiparticle peaks of two isolated atoms are split into bonding and anti-bonding combinations once they form dimer molecules. Therefore, isolated dimers are ideal building blocks for particle-hole 16 symmetric densities of states, such as the top hat and the two-peak target function. The intra-dimer spacing determines the positions of the E = ±t 12 poles via Eq. 2.2. With an appropriate distribution of these distances the full target spectrum of the top hat function can be covered. The poles close to the band center (E = 0) are provided by the less tightly bound dimers. For the two-peak target function the dimer building blocks required to realize the target spectrum are more tightly bound, and the intra- dimer spacings need to vary less to achieve this target. The idea that dimers can be used as building blocks for the particle-hole symmetric target functions strictly applies only to isolated dimers, i.e. the dilute limit, or when potential gradients across dimers pairs from the presence of adjacent atoms do not break particle-hole symmetry. The absence of such gradients, even for relatively high atom densities in a long-range interacting system, accounts for the success of dimers in satisfying the target function. Lower symmetry building blocks, such as trimers and quadrumers shown in Fig. 2.5, can achieve more complex target functions, in particular those with broken particle- hole symmetry. Due to frustration, trimer molecules intrinsically have asymmetric densities of states with unequal spectral weights on the electron and the hole side of their spectra. Quadrumers have a symmetric spectral response in the absence of longer-range frustrating interactions. As an example of how these building blocks enable more complex target functions, let us consider the asymmetric two-peak density of states given by Eq. 2.5 and shown in Fig. 2.2(c), which has a narrow upper 17 peak and a wider lower peak, separated by a gap. In Fig. 2.6 it is observed that an approximate match can be achieved by the adaptive method. As expected, the building blocks for this particle-hole-asymmetric target function are combinations of dimers, trimers, and quadrumers, which partially recombine into larger clusters. Obviously, in order to achieve asymmetric target functions more complex building blocks are required. Especially for systems which contain only a small number of atoms it is important whether the required building blocks are available, and whether there are non-participating unbound atoms that may deteriorate convergence and match to the target function. Figure 2.6: Asymmetric two-peak quasiparticle density of states. The solution contains dimer, trimer, and quadrumer building blocks. Let us finish this section by addressing the convergence properties of the three cases that were discussed. So far, only the simplest guided random walk optimization method was considered in which every “downhill step” is accepted. These “downhill steps” are local updates of individual atomic position that lead to a better match of 18 the spectrum to the target density of states. Especially for the most symmetric target function, this algorithm converges efficiently, with a small remaining error. For the symmetric two-peak function, convergence occurs even faster because the required dimer building blocks are more uniform than for the first case. However, the remaining error is slightly larger, indicating that this may be a metastable solution which could be improved by global updates in which whole subclusters are simultaneously updated. Finally, the convergence plot for the asymmetric two-peak target function (Fig. 2.6(c)) shows several plateaus, indicating metastable configurations that exhibit a high resistance against local updates. For this case, a much larger number of local updates is required to achieve an acceptable match. Hence, more complex target functions clearly call for more sophisticated numerical search tools, including annealing steps, parallelization, and global updating schemes when available. 2.2.4) Adaptive design in discrete configuration space In the previous section local updates were considered which allow atomic positions to change continuously within a given radius. However, for the case of atoms deposited on a substrate with a given lattice structure, the set of available positions is usually discrete, although it may be very large. This has significant consequences for the adaptive design approach. The search space of solutions is finite in this case, which makes it feasible to study more atoms with similar computational effort compared to the continuous case. At the same time, the discreteness of the lattice can prohibit favorable configurations that are available in the continuous case, thus 19 deteriorating convergence properties of the optimization procedure. In practice, we find that the feasibility of computations for larger numbers of atoms in the discrete case helps to achieve better matches, as long as the lattice spacings remain sufficiently small. In this section, we explore the effects of an underlying grid on the adaptive design procedure. Also, more advanced techniques are implemented for the numerical optimization, including hybrids of the genetic algorithm, simulated annealing, and the guided random walk method. Figure 2.7: Adaptive design of clusters with 48 atoms on a discrete lattice. The (a) flat, (b) symmetric two-peak, and (c) asymmetric two-peak densities of states are the same as discussed in the previous section. The corresponding atomic configurations are shown in (d), (e), and (f). 20 In Fig. 2.7, adaptive design results are shown for systems with 48 atoms on a square lattice with spacing 0.01. The length scale is set by the linear size of the two- dimensional box (48 × 48) to which the particles are confined. These systems are in the dilute limit, and hence the convergence to the target functions, chosen to be the same as in the previous section, is good. Also, because of the larger number of particles, there are less finite size effects. For the top hat target function one obtains a small matching error of Δ = 0.000643, for the symmetric two-peak function one finds Δ = 0.000605, and for the asymmetric two-peak function the error is Δ = 0.150347. Thus, analogous to the continuous case, the more symmetric targets are easier to achieve. Again, clustering into dimers is observed for the particle-hole symmetric target functions, whereas trimers are the preferred building blocks for the asymmetric target. Since atomic densities in these examples were chosen to be in the dilute limit, boundary effects and interactions between the building blocks are small. The resulting configurations have the character of liquids, governed by weak interactions between the molecular building blocks, and relatively strong confining forces that lead to the formation of dimers and trimers. 21 Figure 2.8: Effect of the coarseness of the underlying lattice on the formation of building blocks. 32 atoms are confined to a 32 × 32 square box. The target function is the top hat density of states. The grid spacing is varied: (a) 0.01, (b) 0.3, (c) 0.7, and (d) 1.0. Next, let us explore the dependence of these solutions on the coarseness of the underlying lattice. In Fig. 2.8, the grid spacing is varied over two orders of magnitude from 0.01 up to 1.0. As expected, the convergence to the target top hat function deteriorates dramatically as the substrate is made coarser. For the smallest spacing of 0.01 one converges to a final error of Δ = 0.28457, for a spacing of 0.3 the 22 error is Δ = 2.45047, for a spacing of 0.7 the error becomes Δ = 3.81779, and for the coarsest case with spacing 1.0 the error is Δ = 36.3721, indicating failure of convergence. The corresponding configurations in Fig. 2.8 show a strong dependence of the clustering sizes as the system is trying to cope with less available phase space to match the target function. For the finest grid spacing (Fig. 2.8(a)) one observes almost entirely isolated dimers. As the spacing is increased (Fig. 2.8(b)), a few small strings and groups are formed. At grid spacing 0.7 (Fig. 2.8(c)), the solution is made up mostly of long strings and dimers in close proximity to each other. Ultimately, for the coarsest grid spacing of 1.0, (Fig. 2.8(d)), the final configuration consists of square and rectangular blocks of atoms. This result demonstrates that there is a hierarchy of building blocks. The most suitable building blocks for the given target function are dimers. When these become less available due to lattice constraints, the adaptive method selects higher order solutions, i.e. larger size clusters, in order to cope with the more restricted phase space of possible configurations. Hence, for each level of coarseness, adaptive design discovers solutions that enable -up to a given degree of accuracy -a targeted system response. By widening the grid spacings and keeping the box size constant, the density of atoms in the system is effectively increased, and simultaneously the available energy levels are spaced further apart. Let us examine these dependencies by varying the power-law governing the atomic overlaps (Eq. 2.2), and by increasing the number of atoms, i.e. the filling fraction, on a fixed lattice. Results for the achieved convergence are shown in Fig. 2.9. These demonstrate that excellent target matches 23 can be achieved in the dilute limit with filling densities of a few percent. For larger coverages, the numerical search becomes exponentially less effective, indicating increasing frustration effects that need to be addressed by global updating schemes. Higher power-laws imply effective shorter-range atomic overlaps, thus rendering the system more dilute. This is reflected in Fig. 2.9, where the departure from the regime of negligible matching errors is pushed toward higher filling percentages as α in Eq. 2.2 is increased. Some aspects of the long-range tight-binding model on a substrate have already been confirmed experimentally using scanning tunneling microscopy (STM) to precisely position gold atoms on the surface of a nickel-aluminum crystal. In a recent study performed at UCI by Dr. Ho’s group [9, 10], STM measurements show that the splitting in the value of eigenenergies for Au dimers on NiAl depends inversely on Au atom separation corresponding to α = 1. Here, the grid spacings are dictated by the 0.29 nm lattice periodicity of available add-atom sites on the NiAl substrate. Remarkably, the power-law dependence of the effective overlaps t ij between the deposited atoms takes into account interactions with the substrate which are typically difficult to model by first principle computations. 24 Figure 2.9: Effect of the density of atoms and the power-law dependence of the tight- binding overlap integral the convergence of adaptive quantum design. 2.2.5) Summary and conclusions for example 1 In this example, adaptive quantum design techniques were applied to tailor the quasiparticle density of states of atomic clusters, modeled by the long-range tight- binding Hamiltonian. Broken-symmetry spatial configurations of atoms were optimized to match target spectra. By applying adaptive search algorithms, it was shown that matches to target responses can be achieved by forming hierarchies of 25 molecular building blocks that depend on system constraints. For example, symmetric top hat and two-peak target densities of states can be achieved by forming lattices of weakly interacting dimers. While these are the elementary building blocks for particle-hole symmetric case target functions, more complex molecules, such as trimer and quadrumers, are found to dominate the solutions for asymmetric target functions. Implementation of this approach on discrete lattices, corresponding to finite configuration space, introduces frustration effects that destabilize the elementary building blocks in the limit of coarse grid spacings. The core task of adaptive quantum design is the numerical search for global minima in typically shallow landscapes of configurations with many local minima. Since this procedure typically requires many function calls, an efficient implementation on parallel computers is necessary. For the purposes of this example, the complexity of the physical model was minimized in order to limit the computational expense. While the long-range tight-binding model can be viewed as a semi-realistic testing ground for adaptive quantum design techniques, it is crucial to apply these algorithms to more sophisticated models that include, among other ingredients, orbital directionality, spin degrees of freedom, and electronic correlations. Furthermore, adaptive design is applicable to related areas in nanotechnology, including the design of electro-optic components [13] and RF systems [15]. The enthusiasm for adaptive quantum design is, in part, driven by the potential for scientific and technological discovery built into the methodology. It is appealing to 26 perform machine-based searches for new configurations, new components, and new sub-systems, as these will inevitably create new understanding, heuristics, and intuition. The fact that searches are performed in an ultra-large configuration space virtually guarantees the discovery of completely new designs and a much more thorough exploration of the possibilities and capabilities of nanotechnology. 2.2.6) Appendix to example 1: Optimization and search algorithms While this appendix to example one is meant to serve as an overview and benchmark for the search methods presented in this example, a more in-depth discussion of these and other search methods can be found in section 2.5. The optimization algorithms used in this example are based on local updates of atomic positions in order to minimize the error Δ defined in Eq 2.6. Each atom in the system is visited periodically, and a trial change of its position is attempted. Depending on the response in Δ and on the specific algorithm, this trial step is either accepted or rejected. For the case of continuous configuration spaces, these local updates are random shifts of positions within a given radius. A hard core constraint is implemented which forbids atoms to be placed on top of each other. For discrete configuration space, a stochastic distribution function is used to decide which sites in the neighborhood of the original position of an atom should be visited in a trial step. While this function is naturally peaked at nearest-neighbor sites, it has to include a finite probability of longer range updates in order to avoid getting stuck in local minima of search space. The particular results discussed in this appendix are obtained for the continuous case. 27 Figure 2.10: Comparison of the convergence of several optimization algorithms. The error Δ is shown as a function of computer run time. The inset uses a logarithmic scale for Δ. The Newton-Raphson and Broydn methods are multi-dimensional generalizations of the one-dimensional secant method [12]. Unfortunately, their global convergence is rather poor for more than 20 variable parameters. Therefore, they are only applicable to the smallest system sizes, and are not useful for the non-linear multi-parameter searches required for adaptive quantum design. 28 In the guided random walk or “downhill method”, each random step that results in a smaller target error Δ is accepted [13]. Random steps are trial spatial variations about a particle’s position. This guided random walk technique is a quickly implemented power horse. However, especially for shallow landscapes of solutions it gets easily stuck in local minima. Simulated annealing uses an effective temperature representing the likelihood of accepting a step that does not minimize the function. This temperature is lowered slowly - at a rate of 10% - with each iteration. The initial temperature is taken as T init =5t. This method is better at avoiding local minima than the previous techniques. However, it takes a relatively long time to converge. The triggered annealing method is a hybrid of the downhill and the simulated annealing method. It implements the downhill method until minimizing steps become hard to find, at which point the simulated annealing method is used to escape from local minima. Parameters are chosen the same as for the simulated annealing method. Triggered annealing tends to converge relatively quickly if there are only a few local minima. The particle replacement method uses the simple downhill method for guided random updates. In addition, it identifies particles that have not been updated for an extended period, because of being stuck in a local minimum, and assigns them to a 29 new random position within the lattice boundaries. In our implementation, a particle is replaced if there are10 idle iterations without a successful downhill update for that particular particle. Note that this particle replacement update is different from random step updates because it is independent of the previous position of the particle. In genetic algorithms a population of possible solutions is created. Those that best minimize the function are allowed to take part in creating a new generation of possible solutions [14]. These methods are generally good at avoiding local minima, and are also easily implemented on parallel computers. They typically require more function calls than other search algorithms. In order to illustrate the efficiency of these various approaches each method is used to match a flat top hat target function on a one-dimensional lattice with 24 tight- binding atoms and a box size of 96. Particles exiting the box on one end enter it from the other side via periodic boundary conditions. The target density of states is chosen to be the symmetric top-hat function, with poles evenly spaced between E = -4t and 4t. All of these benchmark runs are started with atoms placed randomly along a chain. The computer is a Pentium III 1 GHz with 2GB pc133 memory, and the additional seven nodes used by the genetic algorithm method are Pentium III 850MHz processors with 1GB memory. Each method is run for 600 seconds. As shown in Fig. 2.10, all methods, with the exception of the genetic algorithm, converge rapidly within the first minute of run time, and show only relatively small 30 corrections afterwards. The two techniques with the fastest convergence (inset of Fig. 2.10) are the downhill and the simulated annealing methods. However, their asymptotic error functions remain relatively large, indicating that they get easily stuck in local minima. In contrast, the triggered annealing and particle replacement methods yield much better matches to the target function, while still converging relatively fast. As shown in the inset, the annealing methods sometimes accept trial steps in the “wrong” direction in order to avoid local minima. Finally, the genetic algorithm takes a relatively long time to converge. However, it yields by far the best match to the target density of states after about 7 minutes run time. In order to ensure best matches to the target, this last method is therefore used whenever the computational effort allows it. 2.3) Example 2: Adaptive design of excitonic absorption in broken- symmetry quantum wells This example is presented as a further demonstration of the methodologies introduced in the first example, and uses adaptive quantum design to identify broken- symmetry quantum-well potential profiles with optical response properties superior to previous ad hoc solutions. This technique allows for the engineering of many- body excitonic wave functions and thus provides a new methodology to efficiently develop optimized quantum-confined Stark effect device structures. 31 2.3.1) Physical model: Excitonic absorption in Al x Ga 1-x As quantum wells The absorption of photons in a direct band gap semiconductor involves creation of electron-hole pairs. For uncorrelated electron-hole pairs, absorption occurs at photon energies greater than the band gap energy of the material. However, the coulomb interaction between an electron and hole can form a correlated exciton bound state with an absorption energy less than the band gap energy. The binding energy of this exciton may be increased by confining both the electron and the hole within a two- dimensional quantum well structure [16]. The absorption spectrum that emerges from this quantum well potential profile shows a strong peak just below the band gap energy. This fact has been exploited for use in many modern optoelectronic devices such as modulators and detectors [16, 17, 18]. 32 1.41 1.42 1.43 1.44 1.45 1.46 1.47 Photon Energy, E (eV) 0 2000 4000 6000 8000 10000 12000 14000 Absorption, α (1/cm) 870 860 850 Wavelength, λ (nm) F = 0 kV/cm F = 70 kV/cm F = 140 kV/cm α max -0.50 0.00 0.50 1.00 1.50 2.00 Potential (eV) 0.0 0.5 0.0 0.5 Electron Hole 0 10 20 30 Position, z (nm) -0.50 0.00 0.50 1.00 1.50 2.00 0.0 0.5 0.0 0.5 F = 0 kV/cm F = 70 kV/cm E g = 1.426 eV |ψ h | 2 |ψ h | 2 |ψ e | 2 |ψ e | 2 Figure 2.11: Changes in excitonic absorption spectrum for a symmetric quantum well with increasing electric field. On the left, electron and hole ground state probabilities are shown against the conduction and valence band-edge potential profiles at external field strengths of 0 kV/cm and 70 kV/cm. Decreasing electron-hole overlap and tunneling into the energy gap result in loss of absorption strength and lowered absorption peak energy (shown right) characteristic of the QCSE. Arrows indicate the direction of increasing electric field. Absorption spectra are shown in increments of 10 kV/cm. When a uniform electric field F is applied in the z-direction, perpendicular to the plane of a quantum well, the band-edge potential profile is modified as illustrated in Fig. 2.11. In the figure, the conduction and valence band-edge profiles are shown as a function of position z along with the lowest-energy single-particle probability distribution for electrons, |ψ e | 2 , and for holes, |ψ h | 2 . With increasing applied electric field in the z-direction, electron and hole tunneling results in a decrease of the exciton absorption peak energy. Shifts in the distribution of the electron and hole 33 wave functions lead to a reduction in spatial overlap which in turn reduces the dipole matrix element, diminishing the peak absorption strength. Specializing to the Al x Ga 1-x As material system, the conduction and valence band potential profiles are calculated using the band gap energy g (1.426 1.247 ) eV E x = + (2.7) and a band gap offset ratio between conduction and valence band of 67/33 [26]. The effective electron mass in the conduction band is e 0 (0.067 0.0835 ) m x m = + (2.8) and the effective hole mass in the valence band is h 0 (0.34 0.42 ) m x m = + (2.9) To minimize computation time, the discretization of the quantum well band-edge profile into atomic monolayers is incorporated by using a nearest-neighbor tight- binding Hamiltonian to find the single-particle eigenfunctions in the z-direction. This Hamiltonian for the uncorrelated electrons and holes is given by ∑ ∑ ∑ ∑ + + + + + − = + = > < + + > < + + i i j i i j j i i i j i i j j i c c c c t c c c c t H H H h , h h h h h e , e e e e e h e ) ( ) ( ε ε , (2.10) where t e (t h ) represents the electron (hole) hopping energy, ε e (ε h ) denotes the onsite electron (hole) energy, c e + (c h + ) and c e (c h ) are the electron (hole) creation and annihilation operators, and <i, j> indicates a sum over nearest neighbors only. The single-particle eigenenergies, the electron eigenfunction φ e (z e ), and the hole eigenfunction φ h (z h ) reproduce previously published results [19]. Following Refs. 34 [20] and [19], a variational ansatz is used to obtain the separable 1s exciton wave function ex e h e e h h ( , , ) ( ) ( ) ( ) z z z z ψ ρ φ φ φ ρ = (2.11) in which the in-plane exciton wave function takes the form / 2 2 ( ) e ρ λ φ ρ πλ − = (2.12) where the in-plane coordinate is 2 2 x y ρ = + . The exciton wave function is found by varying the parameter λ to minimize the exciton binding energy using the Hamiltonian 2 h e 2 r 0 2 2 2 h e ex ) ( 1 4 2 z z e H H H − + − ∇ − + = ρ ε πε μ ρ ℏ (2.13) where the first two terms are constant with respect to λ, the third term represents the in-plane kinetic energy of the electron and hole about their center of mass, and the final term is the coulomb potential energy. Here, μ is the reduced mass, r ε is the relative dielectric permittivity, and z e (z h ) is the z position of the center of mass of the electron (hole). This variational method has been shown to accurately model excitonic absorption for quantum well profiles [19]. The energy of the exciton is ex 2 h e 2 r 0 2 ex 2 2 h e ex ) ( 1 4 2 ψ ρ ε πε ψ μλ z z e E E E − + − − + = ℏ (2.14) where the coulomb integral is approximated with a piecewise polynomial found in Ref. [19]. In bulk GaAs, the exciton Rydberg energy, Ry * = E ex is calculated to be 4.3 meV, but due to the confinement of the electron and hole in the z direction, 35 exciton binding energies increase to values near 3×Ry * for structures with well width L = 2.825 nm (10 monolayers), as can be seen in Fig. 2.12. The dependence of the binding energy on the width of the confining potential is approximately logarithmic and extrapolates to a value of 4×Ry * for L = 0.2825 nm (1 monolayer). 0 2 4 6 8 10 12 14 16 18 0 50 100 150 Well Width, L (nm) Binding Energy, E b (meV) Figure 2.12: Exciton binding energy calculated using the tight binding Hamiltonian as a function of quantum well width for Al 0.3 Ga 0.7 As/GaAs quantum well structures. The grey line represents a logarithmic fit to data which approximates the function away from the limits of large and small L. The zero-temperature continuum absorption spectrum may be calculated as a function of photon energy [20] using 36 ( ) 2 0 2 2 2 g g ( ) ( ) | | ( ) ( ) i nm m n E m n i m n i C M E I L E E E E E E E E ω α ω ω Γ = + + + + + + − + Γ ∑∑ ∑ ℏ , (2.15) where L is the width of the quantum well, M(E) is the energy- and polarization- dependant matrix element, |I nm | 2 represents the electron-hole overlap integral, E g is the band gap energy, E i is the energy of the state |i>. C 0 is a constant factor and Γ is a half line width broadening term. Absorption from discrete transitions may be calculated separately with a similar formula [20] and then combined with the continuum result to give the total absorption. This approach accurately captures the main spectral features of other more detailed models [19, 21] of the QCSE while maintaining a low computational cost for calculating the absorption spectrum. As illustrated in Fig. 2.11, for a rectangular quantum well potential profile, the expected loss in absorption strength and shift towards lower energy are obtained in agreement with previous calculations [19]. 2.3.2) Frequency switch: delocalized electron solution For this second example of the adaptive quantum design process, consider the situation where the goal is to design a “frequency switch”, where the absorption peaks at 0 kV/cm and 70 kV/cm are required to have the same strength, but the energy of the absorbed photons is shifted by more than twice the width of the absorption peak (in this case at least 10 meV). In addition, the strength of the absorption peak at each of these two points should be as large as possible. 37 The fitness function for this search is | | 02 . 0 45 1 10 1 | | 2 70 0 70 0 4 70 70 0 E E G − − + + × + − = α α α α α (2.16) where 0 α and 70 α are the strengths of the absorption peak at 0 kV/cm and 70 kV/cm respectively, 0 E and 70 E denote the absorbed photon energies, and all coefficients are set such that on average, each term is of comparable weight. 0.0 0.4 0 2000 4000 6000 8000 10000 12000 Absorption, α (1/cm) F = 0 kV/cm F = 35 kV/cm F = 70 kV/cm F = 140 kV/cm 0.0 0.4 Potential, |ψ e | 2 , |ψ h | 2 (eV) 0.0 0.4 1.42 1.44 1.46 1.48 Photon energy, E (eV) 0 2000 4000 6000 8000 10000 Absorption, α (1/cm) α max F = 0 kV/cm F = 70 kV/cm 0 10 20 30 Position, x (nm) 0.0 0.4 F = 0 kV/cm F = 35 kV/cm F = 70 kV/cm F = 140 kV/cm ΔF = 10 kV/cm Figure 2.13: Broken-symmetry double quantum well obtained from numerical optimization of well width and depth parameters. The target response is an absorption spectrum with maximized, equal-height excitonic peaks at electric fields F = 0 kV/cm and F = 70 kV/cm which are separated in energy by 10 meV. 38 The best solution found by the adaptive quantum design method for this restricted search is shown in Fig. 2.13. It is reminiscent of structures investigated previously [23, 24]. The optimized double well causes the ground state wave function of the hole (broken curve) to develop two maxima whose relative weight is shifted from left to right as the electric field is increased. Simultaneously, the center of the electron wave function (solid curve) moves from left to right, having a maximum spatial overlap with the right peak of the hole wave function at F = 0 kV/cm and with the left peak at F = 70 kV/cm. The resulting exciton peaks in the absorption spectrum, shown on the right, have the desired strength and separation. They are located on two sides of a maximum resonance that is reached at F = 20 kV/cm. In this broken symmetry structure, the maximum of the excitonic absorption peak (α max) initially increases with applied field, and then drops. The corresponding shift of the resonant energy is also nonmonotonic. 2.3.3) Frequency switch: ionized hole solution A second solution found using only two separate regions of constant potential, V(z), is shown in Fig. 2.14. It resembles a conventional quantum well of 38 monolayers, except that near one edge of the well there is a 4 monolayer region of much lower potential. At zero applied field, the hole is localized within this deeper sub-well region, but the electron is delocalized, leading to the initial absorption peak strength. As the field strength is increased, the electron wave function slowly shifts towards the deeper sub-well, while the hole’s position remains pinned. The absorption peak 39 gradually increases until finally the applied electric field is sufficiently large for the hole to ionize out of the sub-well by tunneling to the opposite region of the well. As this ionization process takes place, the electron-hole wave function overlap and hence the absorption peak strength decreases, passing through the same value it held before the field was applied. In this way, the physics behind the trajectory of the absorption peak can be explained in terms of two distinct one-particle stages. 1.36 1.37 1.38 1.39 1.4 1.41 1.42 Photon Energy, E (eV) 0 2000 4000 6000 8000 10000 12000 14000 Absorption, α (1/cm) 910 900 890 880 Wavelength, λ (nm) F = 0 kV/cm F = 70 kV/cm α max -0.25 0.00 0.25 0.50 Potential (eV) 0.00 0.25 0.50 Electron Hole 0 10 20 30 Position, z (nm) -0.25 0.00 0.25 0.50 0.00 0.25 0.50 F = 0 kV/cm F = 70 kV/cm |ψ| 2 |ψ| 2 Figure 2.14: Broken-symmetry quantum well structure obtained from numerical optimization. The target response is an absorption spectrum with strong, equal-strength absorption peaks at 0 kV/cm and 70 kV/cm, separated by 10 meV. The solution involves two single-particle effects: first the electron shifts towards the localized hole, then the hole ionizes out of the deeper sub-well and absorption is suppressed. Electron and hole ground state probabilities are shown on the conduction band-edge profile for visualization purposes. Absorption spectra are shown in increments of 10 kV/cm. 40 Note, that this almost vertical drop of α(E) in a very narrow range of electric fields indicates an exponential sensitivity that is highly desirable for the design of quantum well based modulators. In summary, the adaptive quantum design methodology used in this example reveals design options that have not been explored using conventional ad hoc approaches. The large number of possible solutions along with their exponential sensitivity to small parameter changes render the problem prohibitive to searches by hand [24, 25]. Instead, numerical searches identify optimized broken-symmetry structures which enable desired target functionalities that are useful in the design of quantum-well based photonic switches and modulators. This new approach allows us to engineer many-body excitonic wave functions and thus illustrates a new paradigm in nanoscale design. 2.4) Search: Motivations, terminology, and terrain A typical adaptive quantum design problem begins with a physical model which accurately simulates the physical system under investigation while remaining as computationally efficient as possible. The investigators will decide on a desired output, or “target function”, and then employ advanced optimization methods to iteratively determine the set of input parameters which best approximates the desired result. While the previous examples laid out the procedure involved in taking such 41 an approach, the purpose of this section is to illuminate the nature and difficulty of adaptive quantum design problems in general. 2.4.1) Adaptive quantum design search conditions An adaptive quantum design problem is characterized by a large and complicated search space. That is, the number of possible input parameters must be large enough and the computational cost of each function call must be high enough that a forward approach to the problem is not possible. Also, the complexity of the problem should be such that the results of the physical simulation cannot be approximated by linear interpolation or extraction of trends by random sampling. Having laid out the nature of the problem to be solved, the next point of interest is the spaces and functions involved. 2.4.2) Parameters, spaces, and mapping The parameter space, or search space, S n , is the n-dimensional space of all possible configurations which can be tested or input into the physical model. The size of this space depends on the number of input parameters, n, and the available range of each parameter. This set can either be discreet or continuous, depending on the example [22], but it is almost always bounded due to the finite nature of the problem itself. In example 1, both discreet and continuous search spaces were investigated (on and off- lattice), but in each case the locations of the atoms were bounded by the size of the grid. The physical model acts as a map from the search space, S n , to a subset of the space of all possible observables, P m Õ U m . This mapping is usually not well- 42 behaved, and P m may even be unbounded or disconnected, depending on the nature of the problem. Once S n has been mapped onto P m by the physical model, it must be compared to a human-input target. This target, or objective, T œ U m , specifies a desired value for each of the m observables, and does not necessarily reside in the subset of obtainable solutions, P m . In those cases where m is large and each variable is similar in range, the target can be considered a function of the index of its dimensionality. That is, T = (t 1 , t 2 , t 3 … t m ) can be thought of as T(i), with 1 < i < m. This formulation arises from the fact that for many physical models, the observable variables actually represent discretized values of a function along an axis. In example 1, the target variables t 1 , t 2 , t 3 … t m were the desired state densities at incremental finite energy steps, which together made T(i) = N(E). In order to compare a given set of observables P 0 œ P m to a human-input target T, a cost or fitness function, F(X), must be constructed for X œ U m . The fitness function is a scalar field defined on the space of all possible observables, U m , and not just the space of all outputs from the physical model, P m , because P m cannot be known a priori. At the most basic level, F(X) represents the distance between the two points X and T, and can be visualized as a terrain in U m referred to as the “fitness landscape”[32]. This landscape is a useful concept because even though the physical model may be computationally expensive and impossible to explore exhaustively, 43 the fitness landscape is often much simpler to compute, and can be used to evaluate the effectiveness of both the target T and the fitness function F(X) over the entire range of U m . 2.4.3) Assessing fitness Choosing a fitness function is to a large degree arbitrary, but developing that fitness function into one which will provide the best realizable results during optimization is an involved and deliberate process. The fitness function may need to include various combinations of penalty functions, weights, and even re-casting of variables in order to improve the probabilities of a successful search. Much like the optimization process itself, the fitness function should begin with an initial guess. Assuming we would like to minimize the fitness (rather than maximize), a good initial guess for the fitness function would be ) , D( ) F( i i i t x ∑ = X , (2.17a) where           − −∞ = ∞ = − = otherwise | | if if ) , ( t x t x t x t x D (2.17b) Take, for instance, the fitness function for examples 2 and 3. The observables in this case are P = (ΔE, Δα, α ave ), (2.18) 44 where ΔE represents the change in energy of the exciton peak as the applied field changes from 0 to 70 kV/cm, Δα represents the percent change in absorption at those same points, and α ave represents the average of the absorption strength at those two points. The target in this case is T = (0.02, 0, ¶), (2.19) which implies that a good solution would have a very large change in energy (0.02 eV), minimal change in absorption, and a high average absorption. Applying the guess from Eq. 2.17 to this example, we get F(X) = (0.02-ΔE) + Δα - α ave (2.20) After this, each term should be examined to be sure that the fitness landscape for that particular variable accurately mirrors the desirability of those values. For this example, to improve selective pressure, a square root was used for the ΔE term. For α ave , it was determined that while stronger peaks were more desirable than weaker peaks, the primary goal was to exclude very weak peaks with a large amount of pressure, but not to spend too much computer time on raising peaks which were already of an acceptable height. For this reason, the α ave term was replaced with 1/α ave . The next step is to properly weight each term according to their relative ranges and importance. In this case, each term was considered equally important, so each term was weighted so that on average, when a random result from the physical model, P 0 is plugged in to F(X), all three terms were of the same order of magnitude. This final fitness function used to find the solutions in these examples was 45 E Δ − + × + Δ = 02 . 0 45 1 10 2 2 ) ( F ave 4 α α X . (2.21) Any optimization scheme that uses a fitness function with more than one term, such as the one in Eq. 2.21, is in danger of falling into an undesirable extreme region of the fitness landscape. For example, a solution with the undesirable energy shift of ΔE = 0 but extremely high α ave , and near zero Δα could match fitness with a solution that had acceptable values for each of the three terms. To avoid such unwanted solutions, an acceptable range is assigned to each term in the fitness function. Most iterative search methods gather information from current solutions to generate the next trial set, so deleting solutions with terms beyond the acceptable ranges often results in discarding important information about the location of an optimal solution. Instead, the fitness values for solutions lying outside the accepted range are multiplied by a penalty function. The penalty function is equal to unity over the acceptable range, and increases monotonically outside of that range. The strength and form of the penalty function are chosen to provide sufficient selective pressure to keep the majority of fit solutions within all acceptable ranges. Ultimately, the penalty function(s) must be developed by observing the effects of trial penalty functions on the average relative fitness of each term. For the fitness function in Eq. 2.21, two separate penalty functions are employed for each term of the fitness function. Terms outside of the smaller primary range but inside the greater secondary range are multiplied by a linearly increasing penalty function, while terms outside of the secondary range are multiplied by a quadraticly increasing function. 46 2.4.4) Dimensionality reduction For many physical models, small perturbations in parameter space result in negligible changes in the space of observables, and in some cases this low sensitivity is known to exist over all of parameter space. One such instance is the physical model in examples 2 and 3, where the electron and hole ground state wave functions average out single-monolayer perturbations with energies less than the confinement energy. When this sort of averaging can be assumed, it is possible to adjust the step size of the investigation, or to re-cast the input parameters such that the dimensionality of the parameter space is reduced. For the case of quantum well potential profiles in Al x Ga 1-x As, changes in peak absorption energy greater than half of the peak line width, Γ/2 = 0.007eV and changes greater than 5% in total absorption peak strength require uniform perturbations of 8 monolayers or more in a well 24 monolayers wide. This implies that allowing variation on the level of individual monolayers will not result in any functionality that is no also accessible by varying regions of uniform potential with a width or 8 monolayers or greater. The result of this dimensionality reduction is a dramatic improvement in search results, due to the size of the search space decreasing by many orders of magnitude. If the well width is fixed at 24 monolayers, and aluminum concentration is discretized into steps of 1% with 30% aluminum barriers, there are 30 24 = 2.8x10 35 47 possible configurations in the search space. If, however, changes are made in segments of 8 or more monolayers, there are 30 3 + ((24-8)-8)30 2 = 7200 possible configurations in the search space, where ((24-8)-8) is the total number of possible width combinations for the case of N=2 regions of uniform potential. The possibility of reducing the size of the search space by orders of magnitude makes it advisable to investigate the sensitivity of a given system to perturbations before any large-scale search efforts are launched. 2.4.5) Set-valued target functionals While for most examples, a static target is sufficient to find an acceptable solution, it can be beneficial in many cases to use a more dynamically defined set of targets instead. Consider the search for an electroabsorption modulator, such as the one performed in example 3. A conventional target function would require high absorption values at the beginning of the trajectory, and low absorption values at the end, while varying only slightly in energy. The disadvantage to this static approach is that trajectories which start at low absorption, rise first, and then fall again are excluded, as well as other unanticipated trajectories. The more dynamic approach defines a set of targets based on relative boundaries in absorption and energy space, and then checks to see if any subsection of a given trajectory satisfies the boundary conditions of at least one of the targets in the set. For electroabsorption modulators, a single target is represented by a rectangular region bounded in energy by E and E+ΔE, and bounded in absorption by α max and 48 α max /Δα, and the set of targets used for optimization is the set of these regions for all possible α max and E. Comparing a solution to a target is not performed as a distance calculation here, but rather, the set of targets is examined to determine if the solution’s trajectory satisfies at least one of them. An illustration of this can be seen in Fig. 2.15. The implementation of a dynamic set-valued target system such as this is computationally more taxing than conventional methods, but the fluidity in the method can be worth the effort for cases where desired functionality is less rigid than a point in the space of observables. Figure 2.15: An illustration of a set-valued target functional applied to a specific solution trajectory. In this example, a successful match occurs when the trajectory passes through both upper and lower α boundaries of the rectangular region without crossing through either E boundary. This comparison is performed for all possible locations of the rectangular target region. α E α max α min ΔE ΔF = 10-70kV/cm Δα = α max /α min 49 2.5) Search: Methodology and implementation Choosing the appropriate optimization scheme can be a difficult process of trial and error and often requires meticulous care and thorough sampling of the projection of input parameters onto the fitness landscape. For this reason, the following subsections are meant to serve as a reference guide outlining the implementation and various strengths and weaknesses of many of the more useful algorithms. Assume for all examples below that the fitness function is to be minimized. 2.5.1) Downhill methods The beauty in downhill methods lies in their simplicity. Because they are so easy to implement, it is recommended to try a downhill method before any other methods are considered. A typical downhill search proceeds as follows: 1) Select a random starting point. 2) Make a small random change to a single input parameter. 3) Compare the fitness value to the value before the change was made. If the fitness improved, keep the step, otherwise reject it. 4) Repeat steps 2 and 3 until solution converges. One problem with this method is that if you have a large search space and choose a starting point far away from the minimum, it can take an extremely long time to converge. The primary disadvantage to this method, however, is that it is extremely likely to stop at a local minimum. Most of the following strategies are embellished versions of the downhill method developed specifically to get around the local minima. 50 2.5.2) Simplex methods Not to be confused with the simplex method for linear programming, the downhill simplex method [33, 34, 35] explores n-dimensional parameter space using a simplex, or polytope with n+1 vertices. Here is a sample simplex method algorithm: 1) Evaluate n+1 random points (X 1 , X 2 ,… X n+1 ) in parameter space. 2) Find the point with the highest fitness, and reflect it through the center of the simplex 3) If this new step is now the new minimum point, expand it even further away from the center of the simplex, if it is the new maximum point, contract it closer to the center of the simplex. 4) Repeat steps 2 and 3 until convergence is reached. A good way to visualize this method spatially is to picture a triangle flipping over itself as it tumbles down a hill. While this method has the ability to avoid local minima which are smaller in size than the simplex itself, it is still a local method, and requires a smooth projection of parameters onto the fitness landscape on the order of the size of the simplex. This often results in a simplex too small to avoid local minima when used for nonlinear problems such as those approached for adaptive quantum design. 2.5.3) Gradient methods There exist a wide variety of gradient methods for finding minima [34, 35], most of which become too computationally demanding for problems with more than four 51 input parameters, and many of which provide little, if any improvement over downhill methods in terms of overall speed in finding a solution [22]. These methods (and in particular the Newton method) were designed before computers were used to solve problems iteratively, and therefore focus on reducing iterations in favor of extensive determinate calculation. For this reason, these methods are often favored those who are uncomfortable with stochastic processes, or who do not fully understand the complexity of the problem at hand. Here is an example gradient search algorithm: 1) Pick an initial guess for the solution. 2) Calculate the gradient at that point. 3) Take a step in the direction of steepest descent, using a step size which would lead to the minimum if the function were perfectly quadratic. 4) Repeat steps 2 and 3 until convergence is reached. Two primary faults with this method are that it is often no better at avoiding local minima than the simple downhill approach, and that gradient (and step size) calculation are often too expensive (if possible at all) to calculate in an iterative fashion. Having said this, it should be noted that the gradient information is still extremely useful to any search, so in those cases where this information happens to be inexpensive or when it can be calculated at little extra cost during the calculation of the forward problem [36], it is definitely recommended to use this in the search. 52 2.5.4) Annealing methods Simulated annealing methods [22, 37, 38] are modeled after thermal annealing, a physical process where a substance is heated and then slowly cooled in order to remove crystal defects. The physical process works by allowing ample opportunity for the constituent molecules to find the lowest energy state (a crystal) before the temperature becomes too low to allow for rearrangements. The computational equivalent is some variation on these steps: 1) Pick a starting point randomly 2) Pick a starting temperature, T = T 0 3) Take a random trial step, as in the downhill method. If the step lowers the fitness, keep it. If it increases the fitness, accept the step with probability p = exp(-ΔF(X)/T). 4) Repeat steps 2 and 3, periodically lowering the temperature in small increments, and stopping when convergence is reached. The assumption behind this method is that larger barriers exist further away from the minimum, and that as the minimum is approached, the number and size of these barriers decreases. This assumption is valid for many physical systems in nature, but when simulated annealing methods are applied to those difficult or nonlinear problems for which this assumption is incorrect, little progress will be made towards a nonlocal minimum. Often, the best test for whether or not the annealing assumption applies is to run a simple simulated annealing search and look for any improvement over a downhill search. 53 Many variants exist, most of which find new ways to construct p or lower T. One variation, called “triggered annealing” [22], functions as a downhill method, and only switches to simulated annealing when minima are found, greatly decreasing the search time for problems which deviate from the annealing assumption. The best variant to use depends on the particular problem at hand and should be chosen by running trial searches. 2.5.5) Genetic algorithm methods Genetic algorithms were first used as rudimentary models meant to simulate populations and keep track of changes in phenotype with each generation, but their utility outside the biological sciences was established when John Holland investigated their use for finding the optima of arbitrary mathematical functions [27, 40]. There exist many variations on genetic algorithms, with just as many names, so to disambiguate this section, the term “genetic algorithm” will be used to describe all such strategies. Genetic algorithms selectively combine populations of trial solutions to build new solutions out of successful regions of parameter space, which enables efficient optimization even when local trends do not lead towards a global minimum. This is possible because a genetic algorithm based search is not confined to a neighborhood (however large) around an existing solution, but rather, the search examines the relationship between current parameter values, and builds a new solution from these correlations, entirely independent of proximity in search space. For complicated 54 systems which may exhibit exponential sensitivity to minor changes in input parameters, methods based on neighborhood exploration are seldom able to construct a meaningful downhill trend, leaving neighborhood independent methods such as genetic algorithms as an efficient and reliable alternative. 2.5.5.1) Encoding: Chromosomes The simplest form of the genetic algorithm search is the binary search, which refers to the method used to encode the input parameters. This method encodes a single point in parameter space as a binary array, called a chromosome. The easiest problems to encode in this manor are those with rigid lattice constraints, where each lattice position is represented by a single bit in the array [39]. In these cases, a value of 1 would signify an occupied lattice site, and 0 would signify a vacancy. Most problems, however, are not constrained to a lattice, and have many real valued parameters. The chromosome can then no longer be treated as a binary array and instead contains a single real-valued parameter at each point in the array. The encoding of the parameters is at the same time one of the most important aspects of the search and one of the most overlooked. The primary assumption underlying a genetic algorithm search is that variables or small combinations of variables contribute to the fitness of a solution with some degree of independence from other variables on the array. While the highly popular concept of schemata [27, 40] reinforces this assumption for the case of binary arrays optimizing a one-dimensional 55 parameter space, the schemata explanation becomes misleading for complicated or real-valued problems, and will not be used here in favor of a more general variable- interdependence model which allows for scalability, nonlinearity, and chromosome reordering. 2.5.5.2) Selection schemes The iterative process of optimizing with a genetic algorithm begins with a population of trial solutions, which are given a certain number of copies to contribute to the next generation. These copies are placed in an intermediate generation according to the selection scheme implemented and then operated on by crossover and mutation operators. The purpose of the selection scheme is to ensure that the solutions with the greatest fitness are allowed to have the greatest influence over the next generation, while still avoiding the premature convergence caused by domination of a single solution over the population. Because of this, selection schemes are intimately interdependent on the size of the population, and both should be selected in tandem. The most common method for small binary searches is the roulette wheel method [27, 41]. This method, when used for minimization, randomly selects members of the older generation to be copied into the intermediate population with probability ∑ − ) F( ) F( 1 i j X X , (2.22) 56 where F(X j ) represents the fitness value of the solution being selected, and the sum is over all solutions in the population. Because of this summation, it is obvious that the probability of selecting more than one copy of a successful solution depends directly on the spread of the fitness distribution and inversely on the number of solutions in the population. If this method is to be used, the population size should be selected such that for randomly generated populations, the best solution gets picked multiple times on average when plugged into Eq. 2.22. This size is generally around 20 (depending on the distribution) and generates between 2 and 4 copies of the best solution at each iteration. The roulette wheel selection scheme is primarily used because straight forward to implement and is computationally inexpensive, but when the algorithm begins to converge the fitness value for each solution approaches Δ ± = 0 ) F( F j X , (2.23) where Δ is small compared to F 0 . In this case, the selective power of the scheme diminishes, and every member of the old generation passes on to the intermediate population with equal probability. For this reason, roulette wheel selection is not recommended when the value of the fitness function changes by more than an order of magnitude during the search process. Another commonly used method is the tournament selection method [32, 41], where two solutions are selected at random, and the more desirable solution is placed into the intermediate generation. This selection method has the advantage of maintaining 57 its selective pressure as the fitness distribution converges, but like the roulette wheel method, it provides very few copies of good solutions, and can result in extremely slow convergence for populations of more than 20. The standard deviation method [31, 32] assigns each solution a relative fitness value, f j . If the solution has a fitness value lower than the average, the relative fitness takes on the value 1 ) F( + − = σ μ j j f X , (2.24) where μ is the average fitness value and σ is the standard deviation of the fitness values. This maps the fitness function onto the interval [1, ∞] for solutions which are more desirable than average. If the solution has a fitness value that lies above the average fitness (less desirable), the relative fitness becomes 2 1 ) F( 1         + − = σ μ j j f X , (2.25) which maps these fitness functions onto the interval [0, 1]. Once the relative fitness function has been defined, the intermediate population is filled with one copy of each solution with a relative fitness greater than one. As a copy is placed into the intermediate population, 1.0 is subtracted from the relative fitness value, and the procedure is repeated until there are no longer any solutions with relative fitness greater than one. Each solution is then given a chance to copy itself into the intermediate population with probability f j . If the population size reaches its 58 maximum, the process is halted, and the crossover stage begins. If the process completes without filling all vacancies, any empty slots are filled with copies of the best solution of the generation. The standard deviation selection scheme has the advantage that it maintains selective pressure as the fitness values descend over orders of magnitude, but also carries with it a much higher likelihood of premature convergence due to domination by extremely fit solutions. This sort of domination can occur if one solution is many standard deviations better than any other solution in its generation. In that case, the intermediate population will contain many copies of this single solution, and there is a high probability that subsequent generations will further contribute towards this trend of early convergence. For this reason, standard deviation methods require larger population sizes than other methods. In general, if the ratio of the average fitness value in the first generation to the fitness value of an acceptable final solution is m, and the standard deviation of the first generation is σ, then a population size on the order of 2m/σ will be large enough to avoid domination by poorly optimized solutions. It should also be noted that if there are factors limiting the improvement of the average fitness, such as a physical system where the effects of mutation are so strong as to return most solutions to near-random values, a standard deviation method will seldom avoid premature convergence, and other methods should be considered. 59 There are also many interesting modifications which can be made to a selection scheme to help improve the search. Elitism, or partial replacement, selects a certain portion of the population to remain untouched until the next generation. This practice blurs the lines between generations, and can be a useful tool when the fitness distribution changes drastically from generation to generation. In practice, however, it is often implemented as a result of either an uncomfortability with non- monotonically decreasing fitness functions, or as a crude fix to the problem of having either a population that is too large, or a selection scheme that is providing too little selective pressure. Another important modification is the incorporation of an artificial temperature, much like simulated annealing, and often in the form of a Boltzmann distribution [32]. This can be useful in cases where selective pressure is relatively unimportant at early stages of the search process, for instance, when the search space is large or varies rapidly compared to the number and sampling of a random population, but becomes increasingly more important as a single solution type is chosen. While these options do not by any means exhaust the possibilities for selection, they should provide a broad enough sampling to help in the selection scheme development process. It is important to approach every optimization problem individually, and develop the selection scheme not only around the population size, but the rest of the optimization process as well. 60 2.5.5.3) Crossover Once the intermediate population has been constructed, all intermediate solutions are randomly paired up with each other, and each pair is used to create two new solutions to be placed in the final population. The two final solutions are created by performing crossover and mutation operations on the intermediate pair. The crossover operation randomly assigns elements of each parent chromosome to the two resulting offspring. The most common crossover method is the single point crossover [27], where a point in the chromosome is chosen at random, and the offspring are generated by swapping all chromosome elements before that crossover point. If all elements of the chromosome have only finite-range interdependence on their neighbors, the single point crossover method is a good choice. This method preserves clusters of nearby parameters and passes them on with clusters from the other parent, in hopes that the two successful clusters will combine constructively. A variant of this method is the two point crossover method, where all the information between two randomly chosen points is swapped. This method treats the chromosome as a connected ring, and is ideal for problems with periodic boundary conditions. If strong interdependence exists between elements which are not neighbors in the chromosomal array, either the chromosome should be re-ordered, the crossover should be customized to keep the interdependent variables together, or the parameters should be re-cast in a different basis. Otherwise, clusters of parameters 61 which could add constructively to the fitness of a solution will often be broken, and the crossover operation will be reduced to a randomization of parameters. In the case of binary encoding of occupied states on a two-dimensional lattice, the use of a two-dimensional chromosomal array results in a direct correlation between spatial and chromosomal location. Geometric crossover [39, 42] methods tend to work well in these situations. Geometric crossover methods extend the concept of a swapping point in one dimension to a swapping line in two dimensions or a swapping plane in three dimensions. This extension into higher dimensions, where a chromosomal array with n dimensions is cut by a hyperplane of dimension n-1, demonstrates again the intimate relationship between the encoding of parameter data and the crossover operation. Another useful crossover method is the value-division method, where parameters are assumed to be related by value rather than position in the chromosomal array. This method was used in Example 1 where each element in the chromosome was the x or y coordinate of a particular atom. The largest contributions to the density of states (and the fitness) were due to the interaction between atoms located close to each other, and because the atoms were free to move in 2 dimensions, there was no correlation between spatial and chromosomal proximity. In this case, a line was drawn across the substrate, and all atoms with coordinates below that line were swapped, with exceptions made near the boundary to preserve the total atom number. 62 2.5.5.4) Mutation The mutation step is probably the most flexible step in the genetic algorithm process, and the most understood as well. The purpose of this step is to randomize some members of the population so that a larger portion of parameter space can be incorporated into the search. The mutations rates for various mutation types should be set such that as a solution converges, a significant portion of the population retains a relatively undesirable fitness value. This guideline is based on the assumption that mutations are more likely to damage a solution than improve it, so a significant number of visibly damaged solutions implies a good chance of finding a small number of improved solutions. The ideal rate of mutation operations varies from problem to problem, and should be closely examined to provide a balance between sampling new regions of search space and destroying old solutions. For binary encoded problems, the basic mutation operation is the point, or bit-flip mutation, which flips the value of a particular element from 1 to 0 or from 0 to 1. The real-valued equivalent is perturbing the value of a particular parameter by a random step-size, or alternately, reassigning the parameter with a new random value. Depending on the problem, the effect of this operation on the fitness value can vary greatly, and the mutation rate should be set accordingly, but typical values are on the order of 1% of the total number of bits or parameters. More global mutation operations can be useful when the effects of point-mutation are relatively small. For binary systems, a bit-shift or group-shift mutation will 63 perform a chromosomal translation operation on a 1 (or group of 1s) across an empty (0 valued) array location. Other variations can include rotations of segments in a 2 dimensional array or adding and removing elements such as the atoms used in Example 1. As with all other aspects of the genetic algorithm search method, the individual problem should be examined to determine all relevant and possible mutation types. Once the population has been selected, placed into the intermediate generation, crossed over, and mutated, the newly generated population is then used as the initial population in the next iteration. This process continues until there is no longer any improvement on the best solution, but should always be examined afterwards by downhill exploration to make sure that the solution found is in fact a local minimum. 2.5.6) Hybrid techniques It is often the case that searches of different types (i.e. downhill and genetic algorithm) can be combined to form hybrid techniques [43, 44], which can greatly improve results as well as the rate of convergence. The simplest construction of a hybrid method is to take several downhill steps after the mutation stage of a genetic algorithm, and can greatly improve the trajectory of the fitness distribution for a population. These methods are particularly effective when there are many local minima which are similar in magnitude but not equal, allowing for the minimum of each sampled region to be compared with each other. Depending on the problem, hybrid methods can reduce the randomization generated during mutation, so if a 64 large number of downhill steps are taken, it is advisable to increase mutation rates accordingly, or invert the order and perform the downhill steps first. 2.5.7) Reset techniques All of the above methods are capable of converging to a solution or neighborhood from which no further improvements can be made. Usually, such convergence is considered a sign that a global minimum has been found, but this is not always the case. When convergence is reached, reset techniques [22] can be a helpful explorative tool for determining whether further optimization is still possible. The simplest reset technique to employ is the global reset, which is equivalent to starting the entire search over again from a new random starting point. When the ability to change one parameter is restricted by the values of other parameters, such as in Example 1, where an atom is not able to pass through other atoms on the substrate, one helpful reset technique is particle replacement. This involves identifying the particle(s) which have moved the least in the optimization process and reassigning them to random locations. The optimization will then start over from this new configuration, and approach convergence once again, but with an improved chance of finding a more desirable solution than the solution before the particle was moved. While the methods and techniques presented above do not cover all possible approaches or scenarios, they are to be considered the most likely candidates for a successful search when trying to solve a complex design problem such as those approached in the examples. If there is one resounding theme in these descriptions, 65 it is that no one search parameter or scheme can be picked independently, but rather the search method should be developed in response to the projection of the parameters onto the fitness landscape for the problem at hand. 66 Chapter 3 Optimization without a target function 3.1) Example 3: Quantum well electroabsorption modulator The ionized hole solution to Example 2 met all the target requirements, but displays an interesting behavior in addition to the frequency switching. When the hole ionizes out of the deeper sub-well, the electron-hole overlap diminishes exponentially, which results in a rapid decrease in absorption over a very small range of applied field. This effect shows promise as a high contrast, low chirp, intensity modulator. This section will focus on the development of this discovery into such a device. Optoelectronic intensity modulators made from square wells tend to have low extinction ratios (usually in the range of 4 to 6). In light of this, solutions were sought for a new target function that has the highest possible absorption strength at zero applied field, and then with small applied electric field this absorption strength quickly decreases with negligible shift in frequency and low chirp. For solutions similar to that shown in Fig. 2.14, the rate at which the peak absorption decreases 67 can be approximately measured by examining the fraction of |ψ h | 2 which overlaps the deeper sub well as a function of applied electric field. This function, which can be approximated as a sigmoid, approaches a step function in an ideal switch, and can be characterized by the slope of the best fit line through the center of the step transition. 0.0 0.2 0.4 0.0 0.2 0.4 1.44 1.46 1.48 1.5 1.52 Photon Energy, E (eV) 0.001 0.01 0.1 1 10 100 1000 10000 Absorption, α (1/cm) 860 850 840 830 820 Wavelength, λ (nm) 0.0 0.0 0.2 0.2 0.4 0.4 Potential (eV), |ψ e | 2 , |ψ h | 2 0 10 20 30 40 50 Position, z (nm) 0.0 0.0 0.2 0.2 0.4 0.4 F = 0 - 30 kV/cm F = 31 kV/cm F = -45 kV/cm F = 0 kV/cm F = 35 kV/cm F = 32 kV/cm F = 35 kV/cm |ψ| 2 |ψ| 2 |ψ| 2 Figure 3.1: Broken-symmetry quantum well structure obtained by modifying the structure from Fig. 2.14 to fit a new target functionality: Rapid decrease of the absorption peak with increasing electric field F, small change in absorption peak energy, and stability with respect to crystal growth inaccuracies. The absorption peak diminishes with the application of either positive or negative electric field. Absorption spectra are shown in increments of 1 kV/cm. To aid in the fabrication of devices with this specified potential profile, an additional parameter is added to the search: the zero-field peak absorption energy should be robust with respect to monolayer errors in the potential profile. The conduction band-edge profile, shown in Fig. 3.1, has a sub-well width that initially localizes both 68 the electron and the hole. Robustness with respect to monolayer variations in the potential profile is achieved by making the sub-well shallower and wider than the result illustrated in Fig. 2.14, while still maintaining the ability to localize both particles. This change in form adjusts the sensitivity of absorption peak energy to monolayer fabrication errors from 30 meV to a much more manageable 2 meV. The overall width of the well is also increased so that when the hole ionizes out of the sub-well, there is very little overlap with the electron wave function that remains localized within the sub-well. The ionization of the hole in this solution occurs at 30 kV/cm. The slope at the transition point is θ = 0.27 (kV/cm) -1 , compared to θ = 0.04 (kV/cm) -1 for the solution in Fig. 2.14. Since this solution begins with both the electron and the hole localized within the sub-well, when an electric field is applied in the reverse direction, overlap between the two wave functions again falls rapidly to zero, but in this case it is due to the electron tunneling out of the well instead of the hole. This requires a slightly larger applied field (45 kV/cm compared to 35 kV/cm) but otherwise, the effect is nearly identical. 3.2) Motivations and methodology During the search for an absorption frequency switch the functionality of an intensity modulator emerged without being deliberately sought [30, 31]. This scenario of unintentional discovery prompts the question of how many interesting and useful functionalities exist which have not yet been discovered simply because they have 69 yet to be considered. By harnessing the compute power available today, an intelligent machine-based search may be developed to identify previously undiscovered useful solutions. To investigate this, a new exploration algorithm was implemented to sample solution space randomly and return a set of solutions it found “interesting.” The properties that categorize an interesting solution must be specified prior to the exploration, and should be as general as possible without allowing “uninteresting” solutions into the reported results. The simplest implementation of this exploratory search is the binary test search, which examines random solutions and either saves them for human inspection, or ignores them, depending on whether or not the criteria for an interesting solution have been met. This has the advantage of being simple to implement, but in some cases can result in lists of interesting solutions which are too long to be examined individually by humans, and often will contain many similar or related interesting solutions. Another method for exploratory search is to construct a fitness function based on how interesting a solution is, and optimize to find the most interesting solutions possible. This search should be performed slightly differently from conventional searches like those listed in section 2.5, however because in this case, while the lowest local minima is desired, other local minima should be saved for examination as well, because there may be more than one interesting functionality in the search space. For this reason, the results of many separate optimization efforts should be compounded, and random factors such as temperature or mutation rates should be slightly higher than when a single solution is sought. 70 3.2.1) Defining utility Before any targetless exploration can take place, the qualities which make a solution interesting must be enumerated in some way. While this process is particular to the physical model used, there are some commonalities to the recognition of interesting solutions, which can be addressed here. In most cases, certain ranges of observables will make a solution useless for any functionality, and should not be considered. For the physical model used in examples 2 and 3, any solution which has a low absorption for the entire range of applied field values is effectively transparent and uninteresting. For this reason, in order for a solution to be considered interesting in this scenario, it should have a high absorption for at least one applied field value. When the physical model includes a control parameter in addition to the physical input parameters, such as the applied field in examples 2 and 3, the amplitude of the response to the control parameter is often a good indication of the usefulness of a particular solution. If a solution is unresponsive to the control parameter, it is static and presumable uninteresting, but if there is a large response to a small change in the control parameter, there is a significant chance that the responsive solution could function as a useful device. 71 An initial binary exploration was performed using this physical model where interesting solutions were defined to be those that traversed large distances in either the E or α directions of the E – α(E) plane. A minimum E distance of 0.001 eV (with nonvanishing absorption) or a minimum α distance of 12,000 1/cm would suffice to label the solution “interesting”. Reports of all such solutions (as in Fig. 3.2) were then stored for further examination of these “interesting” features. 1.54 1.56 1.58 1.6 1.62 Photon Energy, E (eV) 0 1000 2000 3000 4000 Absorption, α (1/cm) 800 790 780 770 760 Wavelength, λ (nm) 0 5 10 15 20 25 30 35 Position, z (nm) N=3 Potential Profiles Figure 3.2: Band-edge potential profiles and the corresponding absorption peak paths selected by a machine-based search of configuration space without a human-specified target function. Solutions are sorted from a random set by a binary test for properties which characterize interesting absorption peak paths. In this case, an interesting path is one with a large overall displacement in the E-α plane. Absorption peak paths are tracked as the external applied field increases from 0 kV/cm to 125 kV/cm. 72 Some other characteristics to consider are symmetries, densities, and transitions within the observables. For the physical model in example 1, there is no explicit control parameter, but there is a sequential order to the observables, specifically the 1-dimensional ordering of each finite energy step for which a density of states is calculated. Along this axis, symmetries such as periodicity (in a finite region), and reflection can indicate interesting solutions, as well as areas in energy space with an extremely high density of states, or sharp transitions between areas of high and low densities of states. While the utility implied by various geographical properties such as these vary on a case by case basis, it is hoped that these examples will provide a sufficient starting point for development within other models. 3.2.2) Classifying solutions A key problem faced by both binary and fitness exploratory searches is that of duplicate or near-duplicate solution reporting. This problem occurs when there are a large number of points in parameter space with similar results in the space of observables, allowing for a long list of interesting solutions which differ from each other by only a small amount. Since the purpose of the search is to identify unique candidates for new devices, variations on a single functionality only serve to congest the list of solutions. It is for this reason that the classification of solutions by type is essential to keep the list of solutions manageable. The simplest method of solution classification is to mirror the implementation of the definition of the utility of the solution. In the case of the quantum well structures 73 with an applied electric field, a large change in exciton peak position or in exciton peak amplitude in response to a small change in applied field implied an interesting device. Classification of solutions in this physical model then naturally fell into groups based on the response of these two variables to the applied field, some samples of which can be seen in Fig. 3.3 for structures composed of two equipotential regions (N=2). Figure 3.3: Various categories for exciton absorption peak paths exhaustively sampled for the case of quantum well structures generated from two equipotential regions (N=2). The first arrow in the legend indicates the sign of the first derivative of absorption, α, with respect to applied field, F, and the second arrow indicated the sign of the first derivative of photon energy, E, with respect to F. Comma separated sets of arrows represent consecutive trajectory segments with different first derivatives. For N=1, there is only 1 available category, for N=2, there are 7, and for N=3, there are 11. 74 Another method for categorizing solutions is to isolate similarities and differences in the physical quantities that the observables are constructed from. This is mathematically equivalent to the method illustrated above, but has the advantage of providing intuitive explanations for the differences between functionality categories. Keeping with the absorption peak trajectories, the physical quantities which generate the absorption peak position for a given applied field value are the confinement energies for the electron and hole, the exciton binding energy, and the value of the electron and hole wave functions. An example plot of the electron and hole wave function peak positions for a conventional well structure (N=1) follows in Fig. 3.4. 75 Figure 3.4: The paths of the electron and hole wave function peaks as a function of applied field, F, for a square well 26 monolayers wide. Dark lines represent the peak of the wave function, and thin lines represent the point at which the wave function is equal to half the value at the peak. For N=1, all structures begin with the electron and hole located together and then they steadily separate as in this graph. For all conventional well structures (N=1), the wave function paths are similar. The electron and hole start at the same location, and then separate as the applied field is increased, as in Fig. 3.4. For N=2, there are 7 available categories, 6 of which begin with a segment of increasing absorption rather than decreasing. This difference can be understood by examining the electron and hole paths for structures in these categories. As a result of the asymmetry in the well structure, the hole center of mass is on the lower potential side of the electron center of mass, and as the two 76 peaks approach each other, absorption increases, much like for the structure in Fig. 2.14. Another interesting property that arises from the asymmetry of N=2 structures is tunneling. This appears as a momentary double peak and a sudden change in position on the electron-hole path graph, but is not distinguished in the original first derivative trajectory classification system. Tunneling does not always imply new functionality, though, because if the electron-hole overlap is low, and the tunneling process results in lower overlap, there is little or no effect on the absorption peak trajectory. The identification of effective tunneling processes can be incorporated by checking for a sudden increase or decrease in absorption, and doing so further subdivides the two largest functionality categories in two. An example structure can be seen in Fig. 3.5 which demonstrates asymmetric starting points and tunneling. 77 Figure 3.5: The paths of the electron and hole wave function peaks as a function of applied field, F, for a step well 52 monolayers wide, with a 40 monolayer step of 15% Al. The electron and hole peak positions are initially separated, but as the applied field increases, their two positions approach each other, and then the hole rapidly tunnels out of the sub-well and into the step region. 3.3) Example 4: Double well electroabsorption modulator Drawing from the results of the machine based exploration in Fig. 3.2, an electroabsorption modulator was designed to better separate the electron and hole after ionization, and to better increase overlap before ionization. The solution in Fig. 3.6 consists of two narrow wells of varying potential with a narrow barrier between them. The extinction ratio of this solution is better than that of the well/sub-well solution in Fig. 2.14 by a factor of 100. While this solution looks quite different in 78 form, it works on all the same principles as the well/sub-well solution, and can be generated from the previous solution by adding a small barrier at the well/sub-well interface, and then allowing that barrier to grow in energy to match the potential value of the surrounding material. The narrow wells allow large overlap when the two particles are localized within the same well, and the barrier keeps overlap to a minimum while the two particles are separated. Because of this barrier, the overall structure is smaller in the z-direction than the previous solution. 0.0 0.0 0.2 0.2 0.4 0.4 1.44 1.46 1.48 1.5 1.52 Photon Energy, E (eV) 0.0001 0.001 0.01 0.1 1 10 100 1000 10000 Absorption, α (1/cm) 860 850 840 830 820 Wavelength, λ (nm) 0.0 0.0 0.2 0.2 0.4 0.4 Potential (eV), |ψ e | 2 , |ψ h | 2 0 10 20 30 40 Position, z (nm) 0.0 0.0 0.2 0.2 0.4 0.4 F = 0 - 18.3 kV/cm F = 18.4 kV/cm F = -26 kV/cm F = 0 kV/cm F = 19 kV/cm F = 18.5 kV/cm F = 19 kV/cm |ψ| 2 |ψ| 2 |ψ| 2 Figure 3.6: Asymmetric double quantum well intensity modulator discovered by machine- based exploration without human specification of a target function. Absorption spectra are shown in increments of 0.1 kV/cm. 79 3.4) Example 5: Initially separated double well solution During the categorization of trajectories as in Fig. 3.3 – 3.5, it became apparent that since a high-speed electroabsorption modulator could be built from a trajectory with rapid decrease in absorption due to tunneling, a similar device could be designed from a trajectory which rapidly increases instead. In order to achieve this, the electron and hole wave functions must be initially localized in separate regions of the structure, but once a voltage is applied, the hole would tunnel into the same region as the electron, rapidly increasing the overlap of the two wave functions. This process would be advantageous over previous designs not only because the shifting of the continuum towards longer wavelengths will contribute to a high extinction ratio rather than detract from it, but also because carrier extraction is more efficient at high voltages, making it more beneficial to have high absorption coupled with high applied voltages. Initial separation of the electron and hole wave functions requires at least one narrow sub-well region, typically less than 15 monolayers in an AlGaAs system, in order to exploit the asymmetry between the electron and hole effective masses. In order to isolate the electron and hole wave functions as much as possible before tunneling, a barrier region is also needed. These requirements can both be realized with a structure composed of only 3 equipotential regions (N=3). The first region is a wide well where the electron is localized, the second region is a barrier, and the third 80 region is a narrow well with a lower potential than the first region, which localizes the hole, but which is also narrow enough that the lower potential in the region is much less than the energy required to confine the broad electron wave function. In order to accentuate this effect as much as possible, In 0.27 Ga 0.73 As was used for the 7 monolayer deeper sub-well region, which not only has a lower band gap (1.15 eV) than GaAs (1.42 eV), but is also a strained material, introducing further electron-hole asymmetry. For the 15 monolayer central barrier region, AlAs was used for its high band gap (2.17 eV), and for the wide well region, a conventional 24 monolayer GaAs well was used. The final initially separated double well (ISDW) design can be seen with the electron and hole probability distributions in Fig. 3.7, and the absorption spectrum can be seen in Fig. 3.8. 81 Figure 3.7: Electron and hole probability densities for the ISDW solution. At 0 kV/cm, the electron and hole wave functions are localized in separate wells, due to the difference in electron and hole effective masses. With the application of an applied field, the hole wave function tunnels into the wider well, rapidly increasing the exciton absorption peak strength, which is demonstrated in Fig. 3.8. 82 Figure 3.8: Absorption spectrum for the ISDW solution. At 0 kV/cm, the electron and hole are separated by an AlAs barrier, and there is no discernable exciton peak. As the hole tunnels into the GaAs barrier, the exciton peak rises, and at even higher applied field values, the electron tunnels into the In 0.27 Ga 0.73 As well, again reducing the exciton peak. 83 Chapter 4 Experimental efforts and analysis 4.1) Incorporating fabrication and pragmatics into the model This chapter discusses our attempts to construct and test devices based on the results of the physical model described in examples 2 through 5. While the model used is simple, some pragmatic adjustments were required to better aide in the development of physically realizable devices. For instance, the relative strength of the exciton peak and the continuum varies with the smoothness of the interface between barrier and well layers [45, 46, 47], and that smoothness depends greatly on the growth conditions. For this reason the half-linewidth, Γ/2, and the relative strength of the exciton peak were fit to calibration data rather than calculated from material properties alone. It was also found that the Lorentzian distribution presented with the model in the literature does not fit experimental data as well as a Gaussian, so the model was adjusted to support both peak types. The model was also modified to simulate not only AlGaAs systems, but any combination of the binary compounds GaAs, AlAs, InAs, GaP, and InP, with any 84 binary used as the substrate. Including these new materials requires the use of interpolation of material properties [48, 49, 50, 51], as well as the inclusion of strain effects [52, 53], which are caused by a mismatch in lattice constants. Strain effects are difficult to measure precisely in a laboratory setting, however, and the best interpolative models have been shown to provide inaccurate results in some cases [49], so the interpolations used in our model should not be trusted to provide accurate predictions without calibration. 4.1.1) Poisson solver The well structures in electroabsorption modulators are contained inside the intrinsic region of a PIN diode, which allows for the applied electric field to be adjusted linearly with bias voltage. Provided the background doping levels are low enough, the electric field will remain relatively constant across the multiple quantum well structures. The control parameter for this PIN device is the bias voltage, V bias , which generates the applied field, F, that appears in the physical model. In order to relate our predicted behavior to measured modulation, it is necessary to include a Poisson equation solver for the entire PIN structure. The Poisson solver [54, 55] calculates the band edge of the entire pin structure by iteratively finding a self consistent solution to the equation r 0 d a 2 )) ( ), ( ), ( ( ) ( ε ε ρ z N z N z z Φ − = Φ Δ , (4.1) where Φ(z) is the potential along the z direction, ρ(Φ(z), N a (z), N d (z) ) is the charge 85 density, and N a and N d are the doping concentrations. The charge density, ρ, can be solved numerically using a finite-difference Schrödinger equation solver [54], or by using a Taylor expansion of the exponential representation including the chemical potential at either end of the diode [55, 56], )) ( ) ( ) ( ) ( ( ) ( c d z p z n z N z N e z v a + − − = ρ , (4.2) where e is the electron charge and T k z z e T k z z z e N z p N z n B B max / )) 0 ( ) ( ( a v / )) ( ) ( ( d c e ) ( e ) ( = Φ − Φ − Φ − = Φ − = = (4.3) and z is defined such that the P side of the structure begins at z = 0. The final result of this simulation is that for a given structure, the inherent electric field when V bias = 0 V can be calculated, as well as the linear increase in electric field, F, as a function of V bias . A sample plot of the conduction and valence band edges at V bias = 0 can be seen in Fig. 4.1. 86 Figure 4.1: The conduction and valence band edge profiles for a 6 quantum well AlGaAs/GaAs electroabsorption modulator grown in a PIN layer structure. At V bias = 0V, the electric field across the quantum well structures is -18.6 kV/cm, and for each volt applied, the field increases in magnitude by 12 kV/cm. 4.1.2) Sensitivity analysis for inaccuracies in well thickness While minor perturbations in potential and smoothness at well-barrier interfaces are included by broadening Γ/2, the half-linewidth, inaccuracies in quantum well thickness leads to a shifting in the energy of the exciton absorption peak, and must be considered separately. This occurs because of the shift in electron-hole confinement energies, and the magnitude of this shift decreases for wider well structures. A sample well-thickness analysis for a 24 monolayer Al 0.3 Ga 0.7 As /GaAs conventional well structure can be seen in Fig. 4.2. 87 Figure 4.2: Absorption spectra plotted at F = 0 kV/cm for conventional Al 0.3 Ga 0.7 As/GaAs structures of varying well widths. A well fabricated two monolayers too thick shows a shift in absorption peak energy of ΔE = -7 meV, while a well fabricated two monolayers too thin shifts by ΔE = +10 meV. The well structures grown for experimental verification were all 24 monolayers or larger in size, in order to avoid large unexpected shifts in absorption energy, with the exception of the 7 monolayer narrow well in the ISDW structure of example 5 and Fig. 3.7. At 7 monolayers, growing the structure only one monolayer too thick results in a transition energy shift of -45 meV. The device still functions robustly, however, because ideally there is no absorption from transitions inside the narrow well. The positions of energy levels inside the 7 monolayer well dictate the applied field values necessary to switch the absorption on and off rather than the absorption energy, and appreciable absorption only occurs when both wave functions are localized in the wider GaAs well. In Fig. 4.3, a comparison between transition 88 energies can be seen for both 7 and 8 monolayer well structures, including the GaAs well transition which remains completely unchanged by the error. Figure 4.3: Insensitivity of ISDW structure to monolayer fluctuations, despite narrow well width. The transition energies can be seen for 7 (A) and 8 (B) monolayer structures, but only the robust GaAs transition has a significant overlap integral. The energy levels in the narrow well control the voltage at which the exciton peak switches on, but otherwise have little bearing on the device. InGaAs transition: 1.525 eV (813 nm) Inter-well transition: 1.465 eV (846 nm) GaAs transition: 1.490 eV (832 nm) InGaAs transition: 1.570 eV (790 nm) GaAs transition: 1.490 eV (832 nm) Inter-well transition: 1.482 eV (837 nm) A) B) 89 4.2) Experimental results Experimental efforts included structure growth at three separate facilities, using both Metalorganic Chemical Vapor Deposition (MOCVD) and Molecular beam epitaxy (MBE) growth methods. Structures were grown on GaAs and InP substrates for conventional, step-well, and initially separated double well (ISDW) designs. Measurements of these devices were performed by Wei Hsin, Yuanming Deng, and Donghai Zhu at our own facilities, using a tunable laser to sweep wavelength over the absorption spectrum near the band edge, or using a fixed wavelength integrated laser grown adjacent to the modulator section. The results show that given the accuracy of the growth methods and facilities available, significant improvement on conventional modulation cannot be achieved by control of the exciton peak trajectory alone. 4.2.1) 1550 nm MOCVD Growth InGaAsP integrated devices The first growth to be tested was an InGaAsP conventional device with an integrated laser grown adjacent to the modulator. The absorption spectrum for this structure could not be scanned, but the intensity of the modulated signal was monitored as a function of bias voltage at 1550 nm. The measured absorption coefficient and the simulated curve which was fit to the data are shown in Fig. 4.4. Material composition is In 0.74 Ga 0.26 As 0.73 P 0.27 for the well and In 0.75 Ga 0.25 As 0.49 P 0.51 for the barrier. The well width is 6.5 nm and the band gaps for these materials are 0.786 eV and 0.992 eV respectively, and the half-linewidth is Γ/2 = 0.024 eV. 90 Absorption coefficient measurement for AdTech's EML 0 50 100 150 200 250 300 350 -4 -3.5 -3 -2.5 -2 -1.5 -1 -0.5 0 0.5 V EA (volt) Absorption coefficient (cm-1) I=10 mA I=20 mA I=30 mA I=40 mA I=50 mA I=60 mA I=70 mA I=80 mA I=90 mA I=100 mA Simulation Forward bias Reverse bias Figure 4.4: Measured and simulated absorption coefficients at 1550 nm for a conventional square well device with integrated laser. Material composition is In 0.74 Ga 0.26 As 0.73 P 0.27 for the well and In 0.75 Ga 0.25 As 0.49 P 0.51 for the barrier. The well width is 6.5 nm and the band gaps for these materials are 0.786 eV and 0.992 eV respectively, and the half-linewidth is Γ/2 = 0.024 eV. The next two devices grown were InGaAsP step-well structures similar to the AlGaAs structure in Fig. 3.1. The structure of the first step-well can be seen in Fig. 4.5, while the predicted modulation is plotted in Fig. 4.6. The structure of the step- well and the ionization-induced loss in absorption are similar to the solution in Fig. 3.1, but again, because of the integrated laser, modulation can only be monitored at a single wavelength, rather than observing the entire spectrum. 91 Figure 4.5: Conduction and valence band profiles for InGaAsP Step-Well 1. Material composition is In 0.74 Ga 0.26 As 0.73 P 0.27 for the well, In 0.78 Ga 0.22 As 0.57 P 0.43 for the step, and In 0.75 Ga 0.25 As 0.49 P 0.51 for the barrier. 92 Figure 4.6: Predicted absorption as a function of applied field, F, at 1530 nm for InGaAsP Step-Well 1. The sudden decrease in absorption is due to field induced ionization of the hole, and the subsequent increase in absorption is due to the shifting of the continuum towards longer wavelengths. The predicted absorption at 1530 nm is high at low values of applied electric field, rapidly decreases as the hole ionizes out of the sub-well, and then slowly increases again as the continuum shifts towards longer wavelengths. V bias = 0V corresponds to F = 40 kV/cm, and F increases with V bias by 19 kV/cm for every volt applied. The actual device, when measured, behaves contrary to the initial prediction, as can be seen in Fig 4.7. 93 EAM optical mode absorption coefficient (α α α α=Γ Γ Γ Γ QW N QW α α α α iQW ) vs EAM bias voltage 0 50 100 150 200 250 300 350 -2 -1.5 -1 -0.5 0 0.5 EAM bias (volt) Absorption coefficient (cm -1 ) IDFB=100mA IDFB=90mA IDFB=80mA IDFB=70mA IDFB=60mA IDFB=50mA IDFB=40mA IDFB=30mA IDFB=20mA IDFB=10mA L EAM =200 μm Figure 4.7: Measured absorption as a function of V bias at 1530 nm for InGaAsP Step-Well 1. The modulation is contrary to the predicted response, starting at low absorption and increasing with V bias . The slope of the modulation curve is very steep compared to the conventional device in Fig. 4.4, but because this occurs in forward bias, a high speed (40 GHz) device is not possible. The difference between predicted and measured absorption can be understood in terms of level filling and the location of the first excited state of the hole wave function relative to the ground state. When V bias is near 0V, carriers cannot be efficiently extracted from the PIN structure, and when V bias is in forward bias, carriers are actually being pumped in to the structure, filling low lying energy states, and making the creation of an electron-hole pair at the band edge more difficult. In addition, the first excited state for the hole before ionization lies at the far end of the step, having only negligible overlap with the localized electron. For this reason, the next transition capable of contributing to the absorption before ionization is the second excited state, which has little effect when observing a band-edge wavelength. 94 This explanation supports the measured data, but cannot be quantitatively evaluated because it is not possible to determine the degree to which the low lying states have been filled inside the device. For these reasons, most of our subsequent efforts were spent on devices that switch from low absorption to high absorption and work with the effects of carrier extraction and continuum shifting rather than against them. The next device grown and investigated in the InGaAsP system was an initially separated double well (ISDW) structure similar to the structure in Fig. 3.7. Because of the low barriers inherent in the InGaAsP system and the long minimum layer thicknesses required using MOCVD growth (4.4 nm), strain effects were exploited in addition to effective mass asymmetries to separate the electron and hole wave functions. High compressive strain in the narrow well (1.07%) shifts both the valence and conduction bands towards higher energy, localizing the hole in the narrow well and the electron in the wider well. The conduction band edge and valence band edge profiles are shown in Fig. 4.8. The material composition is In 0.75 Ga 0.25 As 0.49 P 0.51 (-0.14% tensile strain) for the outside barriers, In 0.52 Ga 0.48 As (unstrained) for the wider well, In 0.97 Ga 0.03 As 0.06 P 0.94 (unstrained) for the high inner barrier, and In 0.67 Ga 0.33 As (1.07% compressive strain) for the narrow well. 95 Figure 4.8: Conduction and valence band profiles for InGaAsP ISDW structure. Material composition is In 0.75 Ga 0.25 As 0.49 P 0.51 (-0.14% tensile strain) for the outside barriers, In 0.52 Ga 0.48 As (unstrained) for the wider well, In 0.97 Ga 0.03 As 0.06 P 0.94 (unstrained) for the high inner barrier, and In 0.67 Ga 0.33 As (1.07% compressive strain) for the narrow well. All ISDW structures grown to these specifications resulted in relaxed crystals incapable of modulation. This indicates that the high strain and narrow well width necessary to separate the electron and hole wave functions in the InGaAsP system are beyond current growth capabilities using MOCVD growth methods. For this reason it may be beneficial to reexamine a solution such as this one once significant advancements have been made to improve accuracy and interface smoothness in growth methods. 96 4.2.2) 850 nm MOCVD Growth GaAs device Concurrent with the 1550 MOCVD growth efforts, a step well Al x Ga 1-x As structure similar to the structure in Fig. 3.1 was grown and tested with a laser capable of sweeping over a large wavelength range near the band edge. The barriers for this structure consisted of Al 0.3 Ga 0.7 As, the sub-well was a 7 nm layer of pure GaAs, and the step was a 14.4 nm layer of Al 0.1 Ga 0.9 As. As in the case of the InGaAsP step well in Figs. 4.5-4.7, the step well structure operates against carrier extraction and continuum effects, making modulation ineffective when the exciton peak is small relative to the absorption strength of the continuum. The predicted absorption spectrum for this structure, shown in Fig. 4.9, features a distinct exciton peak at 851 nm, which becomes negligible in absorption strength for V bias less than -1V. 97 Figure 4.9: Predicted absorption spectrum for AlGaAs/GaAs step well structure. At V bias = 0 V, the exciton peak is located at 851 nm. For V bias less than -1 V, this peak disappears. The measured absorption spectrum for this device, plotted in Fig. 4.10, shows two small but distinct peaks at 851 nm and 854 nm, which do not necessarily represent exciton formation. The peak at 851 nm is no longer discernable after Vbias = -1V, and the peak at 854 is no longer discernable after Vbias = -2 V. It should also be noted that the continuum saturates near 849 nm, as opposed to 845 nm as predicted by the model. 98 Low power (0.5 mW) USC step well absorption 0 5000 10000 15000 20000 25000 845 850 855 Photon wavelength (nm) Absorption (1/cm) 0.5 V 0 V -0.5 V -1 V -1.5 V -2 V -2.5 V -3 V -3.5 V -4 V Figure 4.10: Measured absorption spectrum for AlGaAs/GaAs step well structure. At V bias = 0 V, there are two absorption peaks, one at 851 nm, and another at 854 nm. These peaks become indistinguishable from the background at V bias = -1 V and V bias = -2 V, respectively. One possible explanation is that the deeper sub-well was grown two monolayers too wide for several wells in the PIN structure. If the peaks are to be considered exciton peaks, the likely explanation for the existence of two separate peaks is that out of the 6 step wells grown in the PIN structure, several were grown with the proper narrow sub-well width of 7 nm, but several others were grown two monolayers wider than specified, giving a final sub- well width of 7.6 nm. This dual-width growth would explain not only the existence of two distinguishable peaks, but the location of the continuum and the larger bias voltage required for switching the second peak. Further samples could not be grown, however, so comparison to other samples and calibration to a square well could not be performed. 99 4.2.3) 840 nm MBE Growth GaAs devices Highly accurate MBE growth methods were then used to grow conventional square well modulators for calibration purposes. These structures consisted of six 7.3 nm GaAs wells separated by 8.5 nm Al 0.2 Ga 0.8 As barriers, and their predicted absorption spectrum is shown in Fig. 4.11. Figure 4.11: Predicted absorption spectrum for conventional Al 0.2 Ga 0.8 As/GaAs square well structure. At V bias = 0 V, the exciton peak is located at 840 nm, and decreases in absorption strength and shifts towards longer wavelengths as the applied bias increases in magnitude. The measured absorption spectrum did not differ much from the predicted valued at V bias = 0V, but as V bias became increasingly negative, the continuum and exciton peak shifted half the expected distance towards longer wavelength and increased in absorption instead of decreasing as expected. Several separate samples were tested, 100 scanning both wavelength and voltage while holding the other constant, all showing this same increase in absorption strength with applied negative bias. This behavior, shown in Fig. 4.12, is unphysical assuming the device was grown to specification, and resembles the unlikely effects of a voltage applied along the waveguide rather than across it [57]. 6 square wells, sample #2, 0.3mW 0 50 100 150 200 250 300 350 400 450 835 837 839 841 843 845 847 849 Wavelength (nm) Absorption (1/cm) Disconnected +1 V 0 V -1V -2 V -3 V -4 V -5 V -6 V -7 V -7 V (again) Figure 4.12: Measured absorption spectrum for conventional Al 0.2 Ga 0.8 As/GaAs square well structure. The exciton absorption peak and continuum absorption both increase in strength with applied negative bias, contrary to the predictions of the physical model. In this particular sample, the exciton peak can not be distinguished at low bias, but can be separated from continuum effects as it increases in strength. Measurements were taking by sweeping wavelength for each fixed value of V bias . V bias = -7 V was measured first, the voltage was incrementally increased to V bias = +1 V, and then V bias = -7 V was measured again to detect and shift during the measurement process. Two ISDW samples similar to the structure in Fig. 3.7 were grown using the same MBE process and the same PIN structure as the conventional square well above. One sample was grown with a single double-well structure, and the other was grown with three double-wells. The preliminary electroluminescence (EL) data showed 101 peaks at 832 nm, 846 nm, and 814 nm, which allowed for a calibration of the model to better fit the experimental data. The peak at 832 nm corresponds to the transition energy inside the 6.8 nm GaAs well, and the 814 nm and 846 nm peaks correspond respectively to the 2 nm In 0.27 Ga 0.73 As well transition, and the lowest energy cross- well transition. The mass of the heavy hole in the In 0.27 Ga 0.73 As well was known to be inaccurate because the hole dispersion is not parabolic inside strained InGaAs, and while the bandgap of strained InGaAs is well known, the interpolative method of determining band lineup relative to unstrained GaAs is known to be inaccurate [49]. For this reason, these two parameters were varied from their interpolated values until the best possible match of predicted and measured peak locations was achieved. These newly fit values required a shift in the conduction band lineup from 1.3865 eV to 1.3850 eV relative to the InP valence band edge, as well as a change in the heavy hole mass from 0.24m 0 to 0.33m 0 . The absorption spectrum for these devices showed no shift in wavelength with applied voltage, and exhibited only a small increase in the GaAs well transition peak, even for large negative bias values, implying that the effective electric field across the structure remained nearly constant as V bias swept over a large voltage range. This is likely due to the high 4.2 nm high bandgap AlAs barrier restricting the flow of carriers through the device and allowing for a buildup of charge to flatten the band- edge profile even at large negative bias. The undesirable effects of the AlAs barrier on carrier extraction and the band edge profile make this design unsuitable for use as a modulation device, so no further InGaAlAs ISDW structures were grown or tested. 102 Chapter 5 Summary and conclusions The method of computationally solving the inverse problem to design semiconductor devices was presented with the example of adaptive quantum design of atomic clusters, and the stochastic search methods used to investigate the problem yielded not only solutions to the target densities of states, but an understanding of the physics driving those solutions as well. The best solutions could be broken down into individual contributions to the density of states from dimers and small clusters; a simple observation which only became obvious after it was discovered by computer. A further example of these design methodologies was then presented with the adaptive design of excitonic absorption in broken-symmetry quantum wells. A frequency switch was used as a target, and two distinct acceptable solutions were discovered to match that target. The underlying physics, electric field induced ionization, was extracted, and human intuition was again developed following discovery by machine. The physical model used in this example was intentionally simplistic, using the tight binding model in the z-direction, and assuming a 2D hydrogenic exciton in the x-y 103 plane. The calculation of absorption was focused primarily on the lowest energy exciton peak, neglecting many other spectrum characteristics. These simplifications allowed for a relatively fast forward-solve time (typically 10s on a 2.8 GHz machine), at the expense of some accuracy in the absorption spectrum. In practice, the broadening and spread due to fabrication and measuring inaccuracies are large compared to the possible improvements on the physical model, but if growth capabilities improve and more complicated structures become feasible, model improvements may become necessary for accurately predicting the performance of a quantum well device. Specific improvements to consider would be incorporating the band structure, including multiple exciton states and the Sommerfeld enhancement factor into the absorption [20], and developing a more finely tuned material parameter model than the current binary interpolation scheme. Future work on the design and inverse-problem solving approach will be geared towards automation and adaptability. The dimensionality reduction process, where inherent averaging over input parameters is exploited to reduce the search space, could benefit greatly from a more algorithmic or automated approach. The next steps for this effort would likely involve applying the design process to a new physical system, so that the elements of the dimensionality reduction process which are system specific can be distinguished from the elements which are more generalizable. Solution classification can be automized in much the same way. The selection, implementation, and fine-tuning of search methods may also be incorporated into an automated system eventually, but the solution distributions on 104 the fitness terrain, as well as the sensitivity to changes in control and search parameters vary so widely, a large number of physical models would have to be implemented before a generic method-optimizing approach could be developed with confidence. Robustness with respect to fabrication parameters and fabrication errors can also be worked into the fitness function as more problems are approached, but in most cases, this comes at the cost of several extra function calls per iteration. For this reason, it may be best to either limit the incorporation of robustness to solution analysis or include it only periodically within the search process, depending on the distance traveled in parameter space. The parameters, spaces, targets, and costs incorporated in solving an inverse design problem were defined and discussed in order to aide in the future implementation of design problems in new systems. The dimensionality of the parameter space was shown to greatly influence the solvability of the problem, as well as the formulation of both the target and fitness functions. Various search methods were addressed, with particular attention paid to the genetic algorithm and the delicate matter of its efficient implementation. The dependency of search method, forward problem, and fitness terrain was established also for the benefit of future design investigation. As more design problems are approached using these methods, it may be possible to generalize and further develop the notion of the set valued target functional, automate the selection and implementation of the appropriate search methods, and 105 automate the process of dimensionality reduction to ensure the most relevant basis of input parameters. The design for an asymmetric quantum well electroabsorption modulator was developed from the underlying physics discovered by machine during the search for a frequency switch. This unintentional discovery of useful functionality prompted an investigation into the search for interesting or useful device functionality without human specification of a desired target. While further targetless searches did not lead to new functionalities in the case of asymmetric quantum well structures, improved electroabsorption modulator designs were discovered, including the initially separated double well structure, which operates by increasing the electron- hole overlap with applied bias, rather than decreasing the overlap, as in the step well designs. Once these machine based searches without a human specified target have been implemented for multiple physical systems it may become possible to further generalize and automate the search and categorization process, as well as refine the process of defining an interesting solution in any physical model. Finally, efforts towards the physical realization of the electroabsorption modulators designed in the previous sections were presented. Step well structures grown in the InGaAsP system proved to be unsuccessful because the continuum effects were much larger than the effects of the exciton peak, and carrier extraction was more difficult at low bias voltages. ISDW structures proved unsuccessful because the high level of strain required to separate the electron and hole wave functions caused the 106 crystals to relax. The GaAs step well structure faced the same problems with continuum effects, along with difficulties in consistent and accurate growth conditions. The GaAs ISDW structures showed agreement with the physical model at V bias = 0 V, but the high barrier required to separate the electron and hole wave functions also inhibited carrier flow, and resulted in a charge buildup that maintained a near flat band condition. For future work in the area of device fabrications, it would probably be best to focus on improvements in growth and measurement processes, and re-approach these designs after significant advancements have been made. A successful ionizing step well modulator would require complete domination of the absorption spectrum by the exciton peak, which can be achieved by growing narrow wells with atomically smooth layer interfaces, or by developing a low-temperature measurement system. A successful ISDW modulator would require the ability to accurately grow atomically smooth quantum wells at widths less than four monolayers in order to exploit the electron hole effective mass asymmetry, alleviating the need for high barriers or high levels of strain. 107 References [1] T.M. Wallis, N. Nilius, and W. Ho, Phys. Rev. Lett. 89, 236802 (2002). [2] L. Bartels, G. Meyer, and K.-H. Rieder, Phys. Rev. Lett. 79, 697–700 (1997). [3] D. G. Grier, Nature (London) 424, 810-816 (2003). [4] M. Yoshita, H. Akiyama, L. Pfeiffer, and K. West, Appl. Phys. Lett. 81, 49 (2002). [5] J. Thalken, Y. Chen, A. F. J. Levi, and S. Haas, Phys. Rev. B 69, 195410 (2004). [6] W. A. Harrison, “Electronic Structure and the Properties of Solids”, W. H. Freeman and Company (1980). [7] G.V. Nazin, X. H. Oiu, and W. Ho, Science 308, 77 (2003). [8] G.V. Nazin, X. H. Oiu, and W. Ho, Phys. Rev. Lett. 90, 216110 (2003). [9] N. Nilius, T.M. Wallis, M. Persson, and W. Ho, Phys. Rev. Lett. 90, 196103 (2003) [10] N. Nilius, T.M. Wallis, and W. Ho, Science 297, 1853 (2002). [11] T.M. Wallis, N. Nilius, and W. Ho, Phys. Rev. Lett. 89, 236802 (2002). [12] W. H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flammery, “Numerical Recipes: The Art of Scientific Computing”, Cambridge (1986). [13] Y. Chen, R. Yu, W. Li, O. Nohadani, S. Haas, and A.F.J. Levi, J. Appl. Phys. 94, 6065(2003). [14] D. Whitley, Stat. Comput. 4, 65 (1994). [15] I.L. Gheorma, S. Haas, and A.F.J. Levi, J. Appl. Phys. 95, 1420 (2004). [16] D. A. B. Miller, D. S. Chemla, T. C. Damen, A . C. Gossard, W. Wiegmann, T. H. Wood, and C. A. Burrus, Phys. Rev. Lett. 53, 2173 (1984). [17] D. A. B. Miller, D. S. Chemla, T. C. Damen, A. C. Gossard, W. Wiegmann, T. H. Wood, and C. A. Burrus, Appl. Phys. Lett. 45, 13 (1984). 108 [18] X. Huang, A. J. Seeds, J. S. Roberts, and A. P. Knights, IEEE Photonics Technol. Lett. 10, 1697 (1998). [19] P. J. Mares and S. L. Chuang, J. Appl. Phys. 74, 1388 (1993). [20] S. L. Chuang, “Physics of Optoelectronic Devices” (Wiley, New York, 1995). [21] C. Y-P. Chao and S. L. Chuang, Phys. Rev. B 48, 8210 (1993). [22] J. Thalken, Y. Chen, A. F. J. Levi, and S. Haas, Phys. Rev. B 69, 195410 (2004). [23] J. A. Trezza, J. S. Powell, and J. S. Harris, IEEE Photonics Technol. Lett. 9, 330 (1997). [24] N. Susa and T. Nakahara, Appl. Phys. Lett. 60, 2324 (1992). [25] W. Chen and T. G. Andersson, Semicond. Sci. Technol. 7, 828 (1992). [26] G. Ji, D. Huang, U. K. Reddy, H. Unlu, T. S. Henderson, and H. Morkoç, J. Vac. Sci. Technol. B 5, 1346 (1987). [27] David E. Goldberg, “Genetic Algorithms in Search, Optimization, and Machine Learning”, Addison-Wesley, Boston (1989). [28] Adam Prügel–Bennett and Jonathan L. Shapiro, Phys. Rev. Lett. 72, 1305 (1994). [29] C. M. Fonseca and P. J. Fleming, "Genetic Algorithms for Multiobjective Optimization: Formulation, Discussion and Generalization," in Genetic Algorithms: Proceedings of the Fifth Int. Conf., 416 (1993). [30] Jason Thalken, Weifei Li, Stephan Haas, A.F.J. Levi, Appl. Phys. Lett. 85, 121 (2004). [31] Jason Thalken, Stephan Haas, and A. F. J. Levi, J. Appl. Phys. 98, 44508 (2005) [32] Melanie Mitchell, “An Introduction to Genetic Algorithms”, MIT Press, Cambridge, Massachusetts (1996). [33] J. A. Nelder and R. Mead, Comput. J. 7, 308 (1965). [34] Dimitri P. Bertsekas, “Nonlinear Programming”, Athena Scientific, Belmont, Massachusetts (1995). 109 [35] R. Fletcher, “Practical Methods of Optimization”, Wiley, New York (1987). [36] Philip Seliger, et. al., to be published in J. Appl. Phys. (2006) [37] S. Kirkpatrick, C. D. Gelatt, Jr., and M. P. Vecchi, Science 220, 671 (1983). [38] R. A. Rutenbar, IEEE Circuits Dev. Mag., pp. 19-26, Jan. (1989). [39] A. Håkansson and J. Sánchez-Dehesa, Opt. Express 13, 5440-5449 (2005). [40] J. Holland, SIAM J. Comput. 2, pp. 88-105 (1973). [41] D. E. Goldberg and K. Deb, “A comparative analysis of selection schemes used in genetic algorithms,” in Morgan Kaufmann’s “Foundations of Genetic Algorithms”, pp. 69-93, San Mateo, CA (1991). [42] B. R. Moon, Y. S. Lee, C. K. Kim, Journal of Intelligent Manufacturing 9, 401 (1998). [43] J. Renders and S. Flasse, IEEE Trans. Syst. Man Cybern. 26, pp. 243, (1999). [44] Mohamed B. Trabia, J. Mech. Des. 126, 969 (2004). [45] K. Leosson, J. R. Jensen, W. Langbein, and J. M. Hvam, Phys. Rev. B 61, 10322–10329 (2000). [46] M. Sugawara, T. Fujii, S. Yamazaki and K. Nakajima, Appl. Phys. Lett. 54, 1353 (1989). [47] M. Yoshita, H. Akiyama, L. Pfeiffer, and K. West, Appl. Phys. Lett. 81, 49 (2002). [48] T. Ishikawa and J. E. Bowers, IEEE J. Quantum Electron. 30, 562 (1994). [49] J. Minch, S. H. Park, T. Keating, and S. L. Chuang, IEEE J. Quantum Electron. 35, 771 (1999). [50] S. Adachi, J. Appl. Phys. 53, 8775 (1982). [51] P. Lawaetz, Phys. Rev. B 4, 3460 (1971). [52] D. J. Arent, K. Deneffe, C. Van Hoof, J. DeBoeck, and G. Borghs, J. Appl. Phys. 66, 1739 (1989). [53] K. S. Lee, W. S. Han, J. S. Kim, B. Lee, J. Korean Phys. Soc. 37, 587 (2000). 110 [54] I-H. Tan, G.L. Snider, L.D. Chang, and E.L. Hu, J. Appl. Phys. 68, 4071 (1990). [55] F.R. Shapiro, IEEE Transactions on Education, 38, 380 (1995). [56] N. W. Ashcroft and N. D. Mermin, “Solid State Physics”, Thomson Learning, Inc. (1976). [57] D. A. B. Miller, A. C. Gossard, W. Wiegmann, T. H. Wood, and C. A. Burrus, Phys. Rev. B 32, 1043 (1985). 111 Bibliography S. Adachi, J. Appl. Phys. 53, 8775 (1982). D. J. Arent, K. Deneffe, C. Van Hoof, J. DeBoeck, and G. Borghs, J. Appl. Phys. 66, 1739 (1989). N. W. Ashcroft and N. D. Mermin, “Solid State Physics”, Thomson Learning, Inc. (1976). L. Bartels, G. Meyer, and K.-H. Rieder, Phys. Rev. Lett. 79, 697–700 (1997). Dimitri P. Bertsekas, “Nonlinear Programming”, Athena Scientific, Belmont, Massachusetts (1995). C. Y-P. Chao and S. L. Chuang, Phys. Rev. B 48, 8210 (1993). W. Chen and T. G. Andersson, Semicond. Sci. Technol. 7, 828 (1992). Y. Chen, R. Yu, W. Li, O. Nohadani, S. Haas, and A.F.J. Levi, J. Appl. Phys. 94, 6065(2003). S. L. Chuang, “Physics of Optoelectronic Devices” (Wiley, New York, 1995). R. Fletcher, “Practical Methods of Optimization”, Wiley, New York (1987). C. M. Fonseca and P. J. Fleming, "Genetic Algorithms for Multiobjective Optimization: Formulation, Discussion and Generalization," in Genetic Algorithms: Proceedings of the Fifth Int. Conf., 416 (1993). I.L. Gheorma, S. Haas, and A.F.J. Levi, J. Appl. Phys. 95, 1420 (2004). David E. Goldberg, “Genetic Algorithms in Search, Optimization, and Machine Learning”, Addison-Wesley, Boston (1989). D. E. Goldberg and K. Deb, “A comparative analysis of selection schemes used in genetic algorithms,” in Morgan Kaufmann’s “Foundations of Genetic Algorithms”, pp. 69-93, San Mateo, CA (1991). D. G. Grier, Nature (London) 424, 810-816 (2003). A. Håkansson and J. Sánchez-Dehesa, Opt. Express 13, 5440-5449 (2005). 112 W. A. Harrison, “Electronic Structure and the Properties of Solids”, W. H. Freeman and Company (1980). J. Holland, SIAM J. Comput. 2, pp. 88-105 (1973). X. Huang, A. J. Seeds, J. S. Roberts, and A. P. Knights, IEEE Photonics Technol. Lett. 10, 1697 (1998). T. Ishikawa and J. E. Bowers, IEEE J. Quantum Electron. 30, 562 (1994). G. Ji, D. Huang, U. K. Reddy, H. Unlu, T. S. Henderson, and H. Morkoç, J. Vac. Sci. Technol. B 5, 1346 (1987). S. Kirkpatrick, C. D. Gelatt, Jr., and M. P. Vecchi, Science 220, 671 (1983). P. Lawaetz, Phys. Rev. B 4, 3460 (1971). K. S. Lee, W. S. Han, J. S. Kim, B. Lee, J. Korean Phys. Soc. 37, 587 (2000). K. Leosson, J. R. Jensen, W. Langbein, and J. M. Hvam, Phys. Rev. B 61, 10322– 10329 (2000). P. J. Mares and S. L. Chuang, J. Appl. Phys. 74, 1388 (1993). D. A. B. Miller, D. S. Chemla, T. C. Damen, A . C. Gossard, W. Wiegmann, T. H. Wood, and C. A. Burrus, Phys. Rev. Lett. 53, 2173 (1984). D. A. B. Miller, D. S. Chemla, T. C. Damen, A. C. Gossard, W. Wiegmann, T. H. Wood, and C. A. Burrus, Appl. Phys. Lett. 45, 13 (1984). D. A. B. Miller, A. C. Gossard, W. Wiegmann, T. H. Wood, and C. A. Burrus, Phys. Rev. B 32, 1043 (1985). J. Minch, S. H. Park, T. Keating, and S. L. Chuang, IEEE J. Quantum Electron. 35, 771 (1999). Melanie Mitchell, “An Introduction to Genetic Algorithms”, MIT Press, Cambridge, Massachusetts (1996). B. R. Moon, Y. S. Lee, C. K. Kim, Journal of Intelligent Manufacturing 9, 401 (1998). G.V. Nazin, X. H. Oiu, and W. Ho, Science 308, 77 (2003). G.V. Nazin, X. H. Oiu, and W. Ho, Phys. Rev. Lett. 90, 216110 (2003). 113 J. A. Nelder and R. Mead, Comput. J. 7, 308 (1965). N. Nilius, T.M. Wallis, and W. Ho, Science 297, 1853 (2002). N. Nilius, T.M. Wallis, M. Persson, and W. Ho, Phys. Rev. Lett. 90, 196103 (2003) W. H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flammery, “Numerical Recipes: The Art of Scientific Computing”, Cambridge (1986). Adam Prügel–Bennett and Jonathan L. Shapiro, Phys. Rev. Lett. 72, 1305 (1994). J. Renders and S. Flasse, IEEE Trans. Syst. Man Cybern. 26, pp. 243, (1999). R. A. Rutenbar, IEEE Circuits Dev. Mag., pp. 19-26, Jan. (1989). Philip Seliger, et. al., to be published in J. Appl. Phys. (2006) F.R. Shapiro, IEEE Transactions on Education, 38, 380 (1995). M. Sugawara, T. Fujii, S. Yamazaki and K. Nakajima, Appl. Phys. Lett. 54, 1353 (1989). N. Susa and T. Nakahara, Appl. Phys. Lett. 60, 2324 (1992). I-H. Tan, G.L. Snider, L.D. Chang, and E.L. Hu, J. Appl. Phys. 68, 4071 (1990). J. Thalken, Y. Chen, A. F. J. Levi, and S. Haas, Phys. Rev. B 69, 195410 (2004). Jason Thalken, Stephan Haas, and A. F. J. Levi, J. Appl. Phys. 98, 44508 (2005) Jason Thalken, Weifei Li, Stephan Haas, A.F.J. Levi, Appl. Phys. Lett. 85, 121 (2004). Mohamed B. Trabia, J. Mech. Des. 126, 969 (2004). J. A. Trezza, J. S. Powell, and J. S. Harris, IEEE Photonics Technol. Lett. 9, 330 (1997). T.M. Wallis, N. Nilius, and W. Ho, Phys. Rev. Lett. 89, 236802 (2002). D. Whitley, Stat. Comput. 4, 65 (1994). M. Yoshita, H. Akiyama, L. Pfeiffer, and K. West, Appl. Phys. Lett. 81, 49 (2002). 114 Appendix: Quantum well device synthesis program guide This appendix is meant as a reference guide for the suite of functions developed to aid in the optimization and design of quantum well based semiconductor devices. The code itself is commented such that a moderate familiarity with matlab should be sufficient to understand and modify any file which may require improvement for future endeavors. In order to maintain clarity, the following text conventions will be used: matlab code and data file contents >> User input matlab command [optional input] matlab output The functions were all written with ease of use in mind, so this guide should be thought of as a reference rather than required reading. To jump in head first and learn these functions by experience, simply type >> absorption at the matlab command prompt, and experiment with the other functions from there. 115 A-1: Modeling the materials: The first functions to familiarize oneself with are those which recall and calculate the material properties of various semiconducting materials. During the design process, these functions will almost never be called directly, but when working with a particular material, or when fitting data to experimental results, familiarity with the following functions will greatly help to answer questions and solve problems as they arise. The simplest of these functions is Material_data.m, which is called from all other material modeling functions. Material_data.m This function makes no calculations whatsoever, and only serves as an easily accessible table of data for various binary compounds. These values are used in all other functions, so if new experimental data suggests a change in material parameters is in order, the values in this function are the only values which need to be changed. To run the program, type the following command: >> Material_data(binary, [property]) where binary is a string (in single quotes) representing the binary material you are interested in, and property is a string representing the material property you would like to look up. The following binary materials are supported: 116 'InP' 'InAs' 'GaP' 'GaAs' 'AlAs' and the following material parameters can be requested: 'bandgap' Unstrained energy gap, E g (eV) 'epsilon' Dielectric constant, ε r (relative to ε 0 ) 'Me' Electron effective mass 'lambda1' Luttinger Parameter, λ 1 'lambda2' Luttinger Parameter, λ 2 'lambda3' Luttinger Parameter, λ 3 'Mhh' Heavy Hole Effective Mass 'Mlh' Light Hole Effective Mass 'lattice' Monolayer spacing (lattice constant/2) (Å) 'c_band' Position of conduction band (relative to InP valence band) (eV) 'v_band' Position of valence band (relative to InP valence band) (eV) 'C11' Elastic stiffness constant, C 11 'C12' Elastic stiffness constant, C 12 'ac' Hydrostatic deformation potential for conduction band, a’ 117 'av' Hydrostatic deformation potential for valence band, a 'b' Sheer deformation potential for valence band, b 'C_InGaP' Bandgap nonlinearity between InP and GaP 'C_InGaAs' Bandgap nonlinearity between InAs and GaAs 'C_InAlAs' Bandgap nonlinearity between InAs and AlAs 'C_InAsP' Bandgap nonlinearity between InAs and InP 'C_GaAsP' Bandgap nonlinearity between GaAs and GaP 'C_AlGaAs' Bandgap nonlinearity between AlAs and GaAs As an example, here is a request for the energy gap of InP: >> Material_data(‘InP’, ‘bandgap’) ans = 1.3500 AlGaAs.m This function calculates material properties for Al x Ga 1-x As based on the Al concentration, x. All properties are linearly interpolated from the two binary compounds, AlAs and GaAs. This function also has the ability to output a summary or material data for a given concentration, or a row of data which can be copied and pasted directly into material layer data files, which are introduced in the read_layers.m and read_poisson.m function descriptions. To run this function, type >> AlGaAs(x, display) 118 where x is the Al concentration, and display is a string representing either the property to be returned, or the method with which to show a material summary. x may take on any value from 0 to 1, and display may have any value listed under Material_data.m, as well as the following strain and nonlinearity dependant properties: 'strain' Strain 's_C' Strain effect on conduction band (eV) 's_V' Strain effect on valence band (eV) 's_Vl' Strain effect on valence band (light hole) (eV) 's_Mhh' Strained heavy hole mass for strain in the range [0 0.02] 'T_bandgap' Total Energy Gap (including strain effects) (eV) The display string may also be set to one of the following material summary modes: 'view' DEFAULT - display basic properties in command window 'layer' displays info which can be pasted into a material layer data file where the ‘view’ value is assumed if no display string is specified. An example of each possible output type is displayed below. 119 >> AlGaAs(0.1) Al(0.1)Ga(0.9)As Energy Gap: 1.57758 eV Electron Mass: 0.0753 Heavy Hole Mass: 0.41425 Strained Hole Mass: 0.238484 Strain(compressive): 0.000137973 >> AlGaAs(0.1, 'strain') ans = 0.000137973 >> AlGaAs(0.1, 'layer') 1.57758 1.61206 ~~WIDTH~~ 0.0753 0.238484 12.796 InGaAsP.m This function is similar to AlGaAs.m, except that the material modeled is In 1-x Ga x As y P 1-y , so there are four binary compounds to interpolate instead of two. It is also necessary to specify a substrate either by passing it as an argument, or by defining it as a global variable called Substrate. To run the function, type: >> InGaAsP(x, y, [display, Substrate]) 120 where x represents the Ga concentration, y represents the As concentration, display is the same as in AlGaAs.m, and Substrate is a string specifying a binary material defined in Material_data.m. Note that if Substrate is passed as an argument, display must be passed as an argument as well. Several examples are listed below. >> InGaAsP(0.22, 0.57, 'view', 'InP') In(0.78)Ga(0.22)As(0.57)P(0.43) Energy Gap: 0.899146 eV Electron Mass: 0.0687142 Heavy Hole Mass: 0.295518 Strained Hole Mass: 0.239703 Strain(compressive): 0.00302145 >> global Substrate; >> Substrate = 'GaAs'; >> InGaAsP(0.22, 0.57) In(0.78)Ga(0.22)As(0.57)P(0.43) Energy Gap: 1.03469 eV Electron Mass: 0.0687142 Heavy Hole Mass: 0.295518 Strained Hole Mass: 0.235582 Strain(compressive): 0.041256 121 InGaAlAs.m This function parallels the InGaAsP.m function, but for the In 1-x-z Ga x Al z As system. To run the function, type: >> InGaAlAs (x, z, [display, Substrate]) where x represents the Ga concentration, and z represents the Al concentration. The variables display and Substrate are exactly the same as in the InGaAsP.m function. For examples of output, see the InGaAsP.m function listing. A-2: Setting up the quantum well structure: Now that the materials have been modeled, a quantum well structure must be specified so that the absorption can be calculated from it. The quantum well structure is constructed from N rectangular regions concatenated together and bounded on either side by boundary layers. Each rectangular region represents a segment grown from a specific concentration of binary materials and has a width and certain material parameters assigned to it. For example, a conventional Al x Ga 1-x As quantum well structure (N = 1) would consist of one boundary layer (x = 0.3), one well layer (x = 0), and then a final boundary layer (x = 0.3). The boundary layers would have a width of 15 atomic monolayers, the well would have a width of 22 monolayers, and the material parameters for each layer would be calculated based on their composition. The read-in functions of this section are designed to make this 122 quantum well construction process as easy as possible, while also allowing for profiles to be changed and saved with minimal effort. While each read-in function has a different type of input, they all perform the same task: writing a set of global variables to the workspace, to be used later on by absorption.m. This is done so that variables can be easily modified within the workspace after they are read in, and so that the same variables can be used for multiple functions. The global variables defined are: Table A1: Global variables defined before absorption is calculated structure_name string The name of the structure. This will appear in the title of all graphs. Substrate string The name of the binary substrate the device is grown on. Must be defined in Material_data.m width 1×(N+2) array Width of each region (monolayers) E_g 1×(N+2) array Band gap of each region (eV) c_band 1×(N+2) array Position of conduction band relative to InP valence band (eV) M_e 1×(N+2) array Electron effective mass M_h 1×(N+2) array Heavy hole effective mass epsilon 1×(N+2) array Dielectric constant, ε r (relative to ε 0 ) y 1×(N+2) array Defined only for In 1-x Ga x As y P 1-y structures. Specifies As concentration. z 1×(N+2) array Defined only for In 1-x-z Ga x Al z As structures. Specifies Al concentration. x 1×(N+2) array Defined for In 1-x Ga x As y P 1-y , In 1-x-z Ga x Al z As, and Al x Ga 1-x As structures. Specifies Ga concentration for In 1-x Ga x As y P 1-y and In 1-x-z Ga x Al z As structures and specifies Al concentration for Al x Ga 1-x As structures. 123 read_AlGaAs.m This function reads in a text file containing information for each layer of a quantum well structure built from Al x Ga 1-x As. If no input argument is given, the default text file, ‘AlGaAs_data.txt’, is used. The content of the default file, AlGaAs_data.txt, is as follows, and can be copied and edited to create new material data files. Structure_name: Conventional Square Well -----MATERIAL DATA--------------------------------- Substrate: GaAs -----LAYER STRUCTURE------------------------------- X monolayers 0.2 15 0.0 22 0.2 15 The values of the material parameters for each layer are calculated using the AlGaAs function and then assigned to their respective global variables. To run this function, type >> read_AlGaAs ([data_file]) where data_file is a string specifying the location of the data file to be read in. 124 read_InGaAsP.m This function works in exactly the same way as read_AlGaAs.m except that the material system is In 1-x Ga x As y P 1-y instead of Al x Ga 1-x As. To run this function, type >> read_InGaAsP ([data_file]) at the matlab command prompt. If data_file is not specified, the default file, InGaAsP_data.txt, is used. The contents of this file can be seen below: Structure_name: Conventional Square Well -----MATERIAL DATA--------------------------------- Substrate: InP -----LAYER STRUCTURE------------------------------- X Y monolayers 0.2546 0.4909 15 0.2569 0.7323 22 0.2546 0.4909 15 read_InGaAlAs.m This function works in the same fashion as those above for the In 1-x-z Ga x Al z As system. To run this function, type >> read_InGaAlAs ([data_file]) 125 at the command prompt. The contents of the default file, InGaAlAs_data.txt, are included below. Structure_name: Conventional Square Well -----MATERIAL DATA-------------------------------------- Substrate: GaAs -----LAYER STRUCTURE------------------------------------ X Z monolayers 0.8 0.2 15 0.88 0.02 22 0.8 0.2 15 read_layers.m This function reads in layer data for a quantum well structure without specifying a material system, so material parameters must be specified explicitly for each layer of the structure. Because of the large number of input parameters, this function is highly susceptible to human error and is not recommended for any quantum well structure which can be read in by any of the above functions instead. To run this function, type >> read_ layers([data_file]) 126 at the command prompt. The contents of the default file, layer_data.txt, appear below. Structure_name: Conventional Square Well -----MATERIAL DATA--------------------------------- Substrate: InP -----LAYER STRUCTURE------------------------------- E_g c_band monolayers M_e M_h Epsilon 0.9996 1.29 15 0.079 0.24 13.3873 0.7856 1.18 22 0.058 0.24 13.9713 0.9996 1.29 15 0.079 0.24 13.3873 If a single layer of a structure represented in this file is made from a material system described above, this layer may be generated by using the ‘layer’ display option in the AlGaAs.m, InGaAsP.m, or InGaAlAs.m functions. A-3: Calculating the absorption: Once the quantum well structure has been generated, the absorption spectrum can be calculated for that structure as a function of applied field. Again, the function we use generates a set of global variables for use by other programs so that results can be graphed and analyzed at will. This program is complex and flexible, so its 127 description will be broken down into separate sections on input/output and modifications. absorption.m: Input and output This function calculates the wave functions and absorption spectrum for a given quantum well structure at a specific applied electric field, F, with a given exciton absorption peak half line width, Γ. To run absorption.m, type >> absorption([efield, gamma]) where efield is an array of n applied electric field values, F i , in units of kV/cm, and gamma is the half line width, Γ, in eV. If the efield array is not specified, absorption.m looks for a global variable called efield, and uses that variable instead. If that variable cannot be found, it uses the default value of 0:5:160. If gamma is not specified, the global variable gamma is used, and if that variable does not exist, the default value of 0.007 eV is used. Once efield and gamma have been set, the quantum well structure is read in. If the global variables from table 1 have not been defined, the read_layers.m function is called to generate them. Once the wave functions and the absorption spectrums have been calculated for all n efield steps, a set of global variables is created for graphing and analysis. These global variables are: 128 Table A2: Global variables generated during absorption calculation efield 1×n array The input array specifying each value of the applied field, F (kV/cm), for which the absorption spectrum was calculated gamma double The input half line width, Γ (eV), specifying half of the experimentally determined full-width at half-max psi_e2 L×n array Electron probability distribution at each efield step *L denotes the width of the structure (monolayers) psi_h2 L×n array Hole probability distribution at each efield step e_confinement 1×n array Electron confinement energy (eV) h_confinement 1×n array Hole confinement energy (eV) exciton_binding 1×n array Exciton binding energy (eV) exciton_peak 1×n array Location of the exciton absorption peak (eV) spectrum R×(n+1) array The first row contains the x-axis for the absorption spectrum (eV), and the rest contain the absorption, α(E) (cm -1 ), for each efield step. *R is a resolution and window dependant number switching_range 1×2 array The approximated boundary points for the efield range over which switching occurs (kV/cm) It should be noted that the switching_range variable cannot anticipate every functionality type, and may need to be changed by hand after absorption.m has been run in order to be sure that it covers the desired efield range. absorption.m: Modifications There are many other parameters and program variables which are specified within the function itself, for the reason that they should be changed only on rare occasions, and only by a user who has already become very familiar with the program. The comments within the code explain their effect, and will be elaborated upon here. 129 Table A3: Variables, default values, and explanations for absorption.m level_override -1 The absorption continuum is built from the first i wave functions. level_override = -1 means i is calculated automatically. level_override = 0 means there is no continuum, exciton peak only. level_override = i 0 means the first i 0 wave functions are used. lowerbound 0.04 The absorption spectrum is calculated for photon energies in a small window around the exciton absorption peak energy. This specifies how far below the absorption peak energy to begin calculating (eV). upperbound 0.04 This specifies how far above the absorption peak energy to stop calculating (eV). gaussian_peak 1 The exciton absorption peak can be modeled by one of two distribution functions: 1 = Gaussian, 0 = Lorentzian gaussian_continuum 1 The continuum can be constructed from one of two distribution functions: 1 = Gaussian, 0 = Lorentzian cont_max 0.2 The contribution of each transition to the continuum is truncated after a certain point to maintain computational efficiency. If unphysical oscillations are found in the continuum, increasing this value may help. *Warning: changing this may greatly increase computation time. peak_const 0.22 The ratio of the exciton absorption peak strength to the strength of the continuum varies by material, temperature, and other factors, and should be fit to experimental data. The default value was fit for the In 1-x Ga x As y P 1-y system. alpha_const 1.5 The strength of the absorption spectrum is determined up to a multiplicative constant. For any given system, this constant should be fit to experimental data. The default value was fit for the for In 1-x Ga x As y P 1-y system, and also includes the optical confinement factor, Γ. emission_offset Depends on substrate The exciton absorption peak position can be off by a rigid shift in energy due to several physical factors. This variable shifts the peak to match experimental data. InP substrate shift: 0.015 eV GaAs substrate shift: -0.004 eV 130 A-4: Data presentation and analysis The data presentation and analysis functions were all developed to make understanding the physics behind a particular structure easy. They are intended to run in an intuitive fashion, without any need for knowledge of their code or method of operation. This section will cover their use, including more advanced options for in-depth analysis of the quantum well structure and the resulting absorption spectrum. spit_out_data.m This function is called once the absorption.m function has completed its calculations. It prints out a summery of the absorption calculation results, as well as a layer by layer description of the structure itself. An example follows: >> spit_out_data ------------------------------------------------------------ Conventional Square Well ------------------------------------------------------------ Substrate: GaAs Monolayer spacing: 0.282665 nm Switching Range: 0 - 160 kV/cm -------------------------------------------------------------- F = 0 kV/cm: F = 160 kV/cm: Absorption peak (wavelength): 834.391 nm 844.886 nm Absorption peak (energy): 1.48611 eV 1.46765 eV Exciton binding energy: 0.0103335 eV 0.00911237 eV Electron Confinement energy: 0.0579506 eV 0.0508971 eV 131 Hole Confinement energy: 0.0144974 eV 0.00186953 eV -------------------------------------------------------------- LAYER 1: Al(0.2)Ga(0.8)As Energy Gap: 1.73433 eV Electron Mass: 0.0836 Heavy Hole Mass: 0.420521 Strained Hole Mass: 0.238476 Strain(compressive): 0.000275945 Width: 4.23998 nm (15 monolayers) LAYER 2: GaAs Energy Gap: 1.42 eV Electron Mass: 0.067 Heavy Hole Mass: 0.408163 Strain: 0 Width: 6.21863 nm (22 monolayers) LAYER 3: Al(0.2)Ga(0.8)As Energy Gap: 1.73433 eV Electron Mass: 0.0836 Heavy Hole Mass: 0.420521 Strained Hole Mass: 0.238476 Strain(compressive): 0.000275945 Width: 4.23998 nm (15 monolayers) -------------------------------------------------------------- 132 The two applied field values reported in this function are the first value in the efield array (usually 0), and the final value in the switching_range variable. graph_structure.m This function creates a graph of the quantum well structure and includes the ground state energy levels of the electron and heavy hole. A small set or relevant data is also included. This function prints out two sets of data which need to be moved to their ideal aesthetic locations after the graph has been created. The first is the width of each region (except for the boundary layers), and the second is the change in potential at each interface between adjacent layers. To move the data text, click on the ‘edit plot’ tool (the mouse arrow button) and drag them to their appropriate locations. To run this function, type: >> graph_structure An example figure follows for a conventional GaAs quantum well structure after width and interface data has been moved. 133 Fig. A1: Output from the function graph_structure.m after data text has been moved. graph_wavefunctions.m This function plots the probability distributions for the electron and the hole on top of the potential profile at two different electric field steps. To run this function, type: >> graph_wavefunctions([field1, field2]) If no input arguments are given, the first subplot is of the first efield step, and the second subplot is of the last efield step in the switching range. If field1 is specified, it is used for the electric field value for the second subplot. If both field1 and field2 are specified, they are used for the electric field values for each subplot. 134 Fig. A2: Sample figure generated by graph_wavefunctions.m using a conventional GaAs quantum well structure. graph_spectrum.m This function plots the absorption spectrum for the quantum well structure at each efield step. The exciton absorption peak is singled out and labeled with a blue line for points outside the switching range, and a yellow and black line for points inside the switching range. The minimum and maximum values of the x-axis are set by the upperbound and lowerbound variables in the absorption.m function. To run this function, type: >> graph_spectrum([logplot]) 135 where logplot is a control string which sets the y-axis to a log scale if it is set equal to the string ‘log’. The wavelength values are marked along the top of the graph, but the tick marks at the top still correspond to the energy x-axis at the bottom. This is because as of Matlab version 7.0.1, it is not possible to use an alternate x-axis, and the best available work-around for this problem is to overlay separate graphs with separate axes. Fig. A3: Absorption spectrum for a conventional well structure. graph_absorption_at.m This function plots the absorption at a specific wavelength or energy as a function of applied field. To run this function, type: >> graph_absorption_at([position, Vtag]) 136 where position is the wavelength (in nm) or energy (in eV) at which the absorption is plotted, and Vtag, if set to the string ‘v’, indicates that the x-axis should be displayed in volts rather than applied field. This is only possible once poisson.m has been run, and then the applied field values it generates are used to run the absorption.m function. If no input arguments are passed, position is set to the location of the exciton absorption peak at the first efield step. Fig. A4: Absorption at 840 nm as a function of applied electric field for a conventional quantum well structure. emission_table.m This function prints out a table of transition energies between various combinations of electron and hole states. It does not include excitonic effects, and should not be used to fit exciton absorption peaks, but it can help explain minor peaks in PL 137 spectrum when examining preliminary experimental data for possible errors in fabrication. To run this function, type: >> emission_table([nstates, field]) where nstates is the number of eigenstates to calculate, and field is the value of the applied field at which the transitions are to be calculated. By default, the first two eigenstates are calculated, and the applied field is assumed to be zero. An example follows: >> emission_table ------------------------------------------------ Emission table for Conventional Square Well ------------------------------------------------ Transition: eV nm 1e-1h 1.49645 828.629 1e-2h 1.53812 806.177 2e-1h 1.6425 754.945 2e-2h 1.68418 736.263 138 A-5: The Poisson solver Another important piece of the design process is incorporating the entire pin structure into the simulation. This is important for determining the available applied field range as well as the difference in applied field between identical well structures. All of the functions for the Poisson solver operate in the same manor as the absorption functions in order to maintain clarity. read_poisson.m This function is analogous in many ways to the read_layers.m function. It reads in a text file containing material data layer by layer, and creates global variables for the Poisson solver to use. To run this function, type: >> read_poisson ([data_file]) where data_file is the name of the text file to be read in. If no file is specified, the default file, poisson_data.txt, is used. The contents of this file follow as an example. Outer_Structure_name: AlGaAs PIN Structure -----STRUCTURE DATA------------------------------------------- Barrier_width(monolayers): 30 Total_Number_of_Wells: 6 Background_doping(p=1/n=-1/i=0): 1 Background_concentration: 1e15 139 -----LAYER STRUCTURE------------------------------------------ DOPING Con E_g c_band width M_e M_h Epsilon 1 1e19 1.42 1.531 50 0.067 0.40816 13.1 1 1e18 2.35 2.124 800 0.117 0.23844 11.276 0 0 1.89 1.818 80 0.092 0.23847 12.188 -------------------------------------------------------------- | Well Structure Automatically Placed Here | -------------------------------------------------------------- 0 0 1.89 1.818 80 0.092 0.23847 12.188 -1 1e18 2.35 2.124 800 0.117 0.23844 11.276 -1 1e19 1.42 1.531 50 0.067 0.40816 13.1 The layers are much the same as for read_layers.m, except that width is now given in units of nm, DOPING represents the doping type (p=1/i=0/n=-1), and Con represents the doping concentration in cm -3 . If the total number of wells is set to any number greater than zero, it is necessary to set up a quantum well structure. If no well structure has been defined in the workspace, the default function, read_layers.m will be run and the resulting data will be incorporated into the pin structure. The global variables created by read_poisson.m are listed in table A4 below. 140 Table A4: Global variables created by read_poisson.m outer_structure_name string The title of the outer structure. This will appear on graph titles along with the inner structure name. pbg_doping double The background doping type within the structure 1 = P-type doping -1 = N-type doping 0 = no background doping pbg_con double The background doping concentration (cm -3 ) p_doping 1×L array The doping type of each layer 1 = P-type doping -1 = N-type doping 0 = I-type region (background doping only) *L represents the total number of layers p_con 1×L array The doping concentration of each layer (cm -3 ) p_E_g 1×L array The band gap of each layer (eV) p_c_band 1×L array The conduction band position of each layer (eV) p_width 1×L array The width of each layer (nm) p_M_e 1×L array The electron effective mass in each layer p_M_h 1×L array The heavy hole effective mass in each layer p_epsilon 1×L array The relative dielectric constant, ε r , of each layer poisson_solver.m: Input and output This function uses the global variables defined in read_poisson.m to set up and then self-consistently solve Poisson’s equation. The result, much like absorption.m, is to generate a new set of global variables which can be used by various graphing and analysis tools. To run this function, type: >> poisson_solver([volts]) where volts is an array of applied voltages to simulate the structure at. If volts is not specified, and there is no global variable volts either, the default value of 0:0.5:3 is used. The output of this function is a statement of progress printed out for each 141 voltage step, which can be used to evaluate whether the function is working properly, or whether the maximum iterations and tolerance settings should be adjusted. An example follows: >> poisson_solver([0 1 2 3]) Tolerance reached: delta_max = 8.468e-7 lt 1e-6, iter = 20 Tolerance reached: delta_max = 8.974e-7 lt 1e-6, iter = 19 Tolerance reached: delta_max = 7.515e-7 lt 1e-6, iter = 21 Tolerance reached: delta_max = 7.515e-7 lt 1e-6, iter = 23 The phrase “Tolerance reached” means that an acceptable solution was found for that particular voltage step. If the phrase “Max Iterations reached” appears instead, it means that the program has given up without finding an acceptable solution. delta_max represents the largest single difference in potential between successive iterations, and should be less than (lt) the tolerance (here, 1e-6). iter represents the number of iterations required to come within tolerance. If no material data is specified in the workspace, read_poisson.m will be executed using its default values. The potential profile will be solved for each applied voltage indicated, and the global variables denoted in table A5 will be generated for graphical manipulation. 142 Table A5: Global variables created by poisson_solver.m volts 1×n v array The input voltage array (V) *n v represents the number of voltage steps taken phi_c n z ×n v array The potential profile, φ, of the conduction band at each voltage step (eV) *n z represents the number of distance steps taken p_charge_density n z ×n v array The charge density, ρ, as a function of position for each voltage step (C/cm 3 ) p_efield n z ×n v array The electric field, F, as a function of position for each voltage step (kV/cm) p_e_steps 1× n v array The electric field, F, at the midpoint of the structure for each voltage step (kV/cm) p_depletion n d × n v array The length of each depletion region (nm) *n d represents the number of regions with finite charge density poisson_solver.m: Modifications Much like absorption.m, there are many variables defined within the function that effect the way the function runs. These variables are explained within the code itself, and should only be changed by a user who is familiar with the function. Their effects and default values are explained in table A6. 143 Table A6: Variables, default values, and explanations for poisson_solver.m Max_iter 100 The maximum number of self consistent iterations to perform before exiting. Most structures will not require more than 50 iterations, even with high accuracy. tol 1e-6 The maximum allowed difference between potential profiles in successive iterations for a solution to be considered converged (eV). Lower values indicate slower convergence, but higher accuracy. distance_steps 500 The potential profile is reformulated into a set of distance_steps evenly spaced steps instead of a set of layers with finite widths. Higher values increase resolution, but with a large increase in calculation time. N_override 0 If you are using a PN, NP, PIN, or NIP structure, leave this setting to 0. If your structure does not have a clearly defined low-potential side, such as a PIP structure, you can override the automatic n- region calculation by specifying which layers are held at lower potential here. P_override 0 If you are using a PN, NP, PIN, or NIP structure, leave this setting to 0. If your structure does not have a clearly defined high-potential side, such as a NIN structure, you can override the automatic p- region calculation by specifying which layers are held at higher potential here. pgraph_phi.m This function graphs the potential profile of the structure as a function of distance. To run this function, type: >> pgraph_phi([which_v_step, show_valence]) where which_v_step indicates which voltage step to plot, and show_valence indicates whether or not to include the valence band in the plot. If which_v_step is not 144 specified, or if it is set to 0, ‘a’ or ‘all’, the potential profile, φ, will be plotted for each voltage step. Otherwise, only the voltage step indicated by which_v_step will be plotted. When only one voltage step is plotted, the change in potential across the structure is noted on the graph, as well as the slope of the potential profile at the midpoint. Two examples follow with their resulting graphs: >> pgraph_phi() Figure A5: The potential profile for the conduction band at each applied voltage step. 145 >> pgraph_phi(1,1) Figure A6: The potential profile for the conduction and valence bands with zero applied voltage. pgraph_efield.m This function graphs the electric field of the structure as a function of distance. To run this function, type: >> pgraph_efield([which_v_step, show_wells]) 146 where which_v_step indicates which voltage step to plot, and show_wells indicates whether or not superimpose the well structure in the i-region for aesthetic purposes. If which_v_step is not specified, or if it is set to 0, ‘a’ or ‘all’, the potential profile, φ, will be plotted for each voltage step. Otherwise, only the voltage step indicated by which_v_step will be plotted. When only one voltage step is plotted, the slope of the electric field in the i-region is annotated on the graph. Two examples follow with their resulting graphs: >> pgraph_efield() Figure A7: The electric field profile for the structure at each applied voltage step. 147 >> pgraph_efield(1,1) Figure A8: The electric at zero applied volts with the well structure superimposed. pgraph_charge.m This function graphs the charge density of the structure as a function of distance. To run this function, type: >> pgraph_charge([which_v_step]) where which_v_step indicates which voltage step to plot. If which_v_step is not specified, or if it is set to 0, ‘a’ or ‘all’, the potential profile, φ, will be plotted for 148 each voltage step. Otherwise, only the voltage step indicated by which_v_step will be plotted. One example follows with the resulting graph: >> pgraph_efield() Figure A9: The charge density as a function of position for the structure at each applied voltage step.

Abstract (if available)

Abstract Synthesis of semiconductor device design requires access to realistic physical models and adaptive algorithms. To demonstrate that such synthesis is feasible, atomic clusters with specified quasiparticle densities of states are designed using advanced search algorithms, and broken-symmetry quantum-well potential profiles are developed with optical response properties superior to previous ad hoc solutions.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Expanding the library of surface ligands for semiconductor nanocrystal synthesis and photovoltaic applications

PDF

III-V semiconductor heterogeneous integration platform and devices for neuromorphic computing

PDF

Perovskite chalcogenides: emerging semiconductors for visible to infrared optoelectronics

PDF

Topological and non-Hermitian semiconductor lasers

PDF

Modeling graphene: magnetic, transport and optical properties

PDF

Advanced cell design and reconfigurable circuits for single flux quantum technology

PDF

Synthesis, characterization, and device application of two-dimensional materials beyond graphene

PDF

Electronic design automation algorithms for physical design and optimization of single flux quantum logic circuits

PDF

Computationally efficient design of optimal strategies for passive and semiactive damping devices in smart structures

PDF

GaN power devices with innovative structures and great performance

PDF

Designing efficient algorithms and developing suitable software tools to support logic synthesis of superconducting single flux quantum circuits

PDF

Reactive and quantum molecular dynamics study of materials: from oxidation and amorphization to interfacial charge transfer

Asset Metadata

Creator Thalken, Jason (author)

Core Title Search, synthesis, and analysis for semiconductor device design

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Physics

Publication Date 10/06/2006

Defense Date 09/22/2006

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag genetic algorithm,inverse problem,OAI-PMH Harvest,optimization

Language English

Advisor Haas, Stephan W. (committee chair), Bickers, Nelson (committee member), Bozler, Hans M. (committee member), Levi, Anthony F. J. (committee member), Lidar, Daniel (committee member)

Creator Email thalken@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-m87

Unique identifier UC1157594

Identifier etd-Thalken-20061006 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-19832 (legacy record id),usctheses-m87 (legacy record id)

Legacy Identifier etd-Thalken-20061006.pdf

Dmrecord 19832

Document Type Dissertation

Rights Thalken, Jason

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Repository Name Libraries, University of Southern California

Repository Location Los Angeles, California

Repository Email cisadmin@lib.usc.edu