Development of a Multi-Objective Evolutionary Algorithm for Strain-Enhanced Quantum Cascade Lasers

An automated design approach using an evolutionary algorithm for the development of quantum cascade lasers (QCLs) is presented. Our algorithmic approach merges computational intelligence techniques with the physics of device structures, representing a design methodology that reduces experimental effort and costs. The algorithm was developed to produce QCLs with a three-well, diagonal-transition active region and a five-well injector region. Specifically, we applied this technique to AlxGa1−xAs/InyGa1−yAs strained active region designs. The algorithmic approach is a non-dominated sorting method using four aggregate objectives: target wavelength, population inversion via longitudinal-optical (LO) phonon extraction, injector level coupling, and an optical gain metric. Analysis indicates that the most plausible device candidates are a result of the optical gain metric and a total aggregate of all objectives. However, design limitations exist in many of the resulting candidates, indicating need for additional objective criteria and parameter limits to improve the application of this and other evolutionary algorithm methods.


Introduction
In the last two decades, quantum cascade lasers (QCLs) have become the standard for high-efficiency emission sources in the mid-infrared (mid-IR) and THz spectrums [1,2].QCL radiative transitions occur between conduction-band subbands formed in quantum wells, making QCLs unipolar devices-relying only on electrons as carriers.These devices offer reduction in efficiency loss, for example by reducing Auger recombination, and have higher tunability than their interband counterparts.However, traditional QCLs use material combinations in lattice matched or strain-compensated regimes that limit their useful emission range to wavelengths greater than 3 µm.Specifically, GaAs based devices are generally limited to emission greater than 5-7 µm because of the low conduction band offset, ∆E c , available when matched with AlGaAs.
Applications such as infrared counter measures, trace gas detection, and free-space communication require high-efficiency and high-power emission and would benefit from operation in the first atmospheric transparency window between 3-5 µm.However, few QCL devices emit in this range.Strain-compensated InAlAs/InGaAs devices grown on InP substrates, traditionally used for these applications, have limited tunability due to a fixed ∆E c .One method we have proposed is the use strained active region designs [3][4][5] with indium incorporation in the wells.While these devices can be grown on the traditional (100) surface, indium incorporation is limited to only a few percent due to strain-induced three-dimensional growth modes [4].However, pseudomorphic growth of (Ga)InAs on GaAs (111)B is possible [4,5], opening up a new class of devices.These devices benefit from strain-induced piezoelectric effect: reducing threshold voltage, extended wavelength tunability, and enhanced non-linear susceptibility [3,5].
The large range of design possibilities using this new material system necessitates an approach that ensures efficient, high-power capable devices.Our approach, as reported here, is the development of a multi-objective evolutionary algorithm for optimized and intelligent search of this design space.This algorithm uses a self-consistent Schrödinger-Poisson solver (nextnano 3 ) for simulation of the QCL device candidates [6].In this work, nextnano 3 was configured to compute device structures using an effective-mass approximation of the envelope function with plots of the relevant moduli squared of the wavefunctions and Dirichlet boundary.Also incorporated in this simulation are calculations of direct and indirect energy gaps, their temperature dependencies, and conduction and valence band deformation potentials accounting for strain effects in pseudomorphic thin layers [7]. Figure 1 shows a general flow diagram of this algorithm.The nextnano 3 simulation tool provides electron energy levels, wavefunction probability (Ψ 2 ), carrier lifetimes, dipole matrix element data to the algorithm.This data is used to evaluate candidates and calculate fitness values based on a multi-objective scheme rooted in the physics of QCL device operation.

Background of Evolutionary Algorithms and QCLs
Evolutionary computation is a set of optimization algorithms that use the basic processes of Darwinian nature: selection, breeding, mutation, inheritance, etc. to make decisions about a population based on objectives.This process is iterative and the intended outcome is progress, based on the objectives, in the population.
In this work, the techniques of evolutionary computation, specifically multi-objective evolutionary algorithms (MOEAs), are used to search the massive design space of the diagonal transition QCL.Our current focus is on strain-enhanced (Al x Ga 1−x As/In y Ga 1−y As) active region designs; however, this algorithmic design approach is much more broadly applicable.As described in greater detail in the following section, the gene representation is real-valued, classifying our approach as an evolutionary algorithm rather than a traditional genetic algorithm that uses binary gene representation.Otherwise, the techniques are most similar to the second generation non-dominated sorting genetic algorithm (NSGA-II) [8].
Recently, some research groups have applied the techniques of global and direct optimization to the design of QCLs [9][10][11][12][13][14][15][16][17][18][19][20][21][22][23].These include use of iterative or direct optimization [9][10][11], simulated annealing algorithms [12,13], single-objective genetic algorithms [14][15][16][17][18][19][20][21], particle swarm optimization [22], and multi-objective genetic algorithms [23].QCL structures are complex multi-layer systems and optimization requires a scheme capable of handling a large number of variables with complex interactions.Direct and iterative optimization methods lack the ability to efficiently search the parameter space.Likewise, simulated annealing methods, while capable of global search, are inefficient for a large number of variables and do not take advantage of parallel computing.Evolutionary algorithms such as genetic algorithms and particle swarm optimization are much more efficient methods for the large parameter space of a QCL design and are well suited for parallel computation.
The key to a successful optimization algorithm is the choice of the objective or merit function(s): the trait(s) being evaluated that classify the fitness of a design.Almost all of the current studies used a form of gain [9][10][11][12][13]15], or a measure related to gain such as wall-plug efficiency [16], population inversion [17,21,23], oscillator strength [17,21], and period length [17].Some groups used an objective function specific to the design, such as non-linear susceptibility (for second-harmonic generation) [18,19] and third harmonic power [22].Only one study is known to have used a multi-objective, non-aggregated, approach [23].However, this group's apparent interests lie in the comparison of evolutionary algorithm techniques rather than the resulting QCL.No evolutionary programming based studies are known to date that incorporate the full extent of design principles from a multi-objective approach.
Many of the groups limited their optimization algorithms to the active region of a QCL design [9,11,12,15,[17][18][19]21].There is a delicate interplay between the injector and the active region and these groups must either set strict limits on the position of the upper state energy level E U and the lower state energy level E L to suit a given injector design or will need to perform additional optimization to find a suitable injector.The coupling between the lower state and the injector is especially critical in our chosen material system.It is possible that the lower energy state will lie below the conduction band of the adjacent injector region well.This necessitates that the injector region be a integral part of the design and optimization methodology.
Various methods for calculating the relevant physics have been used including rate equation modeling [10,12,13,15,[17][18][19], density matrix transport model [11,16], Monte Carlo methods based on Boltzmann transport [21], Pauli master equation [22], and Schrödinger(-Poisson) solver [9,14].Each of these methods have merits for use when considering the delicate balance between accuracy when compared to experiment and computational efficiency.The main objective of this study is to evaluate the applicability and effectiveness of the various objective functions to produce feasible QCL designs.This goal balances toward the use of a computationally efficient method; however, our target material system, which is a strained system, was heavily considered.
The use of global search, especially evolutionary algorithms, is well suited for any choice or even variation of the QCL material system.However, in our case, the use of a strained material system necessitates the use of additional objective evaluation criteria, as is discussed in greater detail in the following sections.

Algorithm Description and Design Approach
The flow sequence of the algorithm is shown in Figure 1.The general parameters of this algorithm are shown in Table 1.The application specific parameters are shown in Table 2.The design of interest is the traditional three-well, diagonal-transition active region.Five (5) injector well/barrier pairs and three active region pairs makes for a total of 16 layers that may be varied in the device.The other two variables are the alloy concentrations of each well and barrier.Only the active region may contain indium.Injector or active region doping is not considered in this work.
The data consists of a set of runs with an initial (training) population and a randomly generated initial population.This makes for a potential for 6400 individual candidates tested in each run.However, this algorithm includes elitism via carryover (the copying of the most fit candidates unchanged into the next generation); therefore, the number of unique candidates is lower.

Initial Population
An initial population of QCL candidates is generated from a set of eight parents selected from known working QCLs produced by groups including Faist, Razeghi, Sirtori, etc. [24][25][26][27][28][29][30] or as a randomly generated set of parents.The initial stage of the algorithm converts this initial population into a gene sequence, illustrated in Figure 2.Each parent is crossbred, using a real-valued crossover operation, generating a population of 64 candidates.These resulting candidates undergo random mutation, at a rate of 2%. Figure 3 gives a flow sequence of the conversion of the initial population.
The initial population gene sequence is given as values of the alloy concentration percentage and thickness of each well and barrier in nanometers.The gene sequence in this algorithm is a real-valued sequence of integers.A conversion of the lengths takes into account the alloy percentages of the well and barrier and converts this into a lattice constant a using Vegard's Law, resulting in monolayers (ML).

Simulation and Data Collection
After generating the initial population, each candidate's gene sequence is converted into a nextnano 3 input file.nextnano 3 is a Schödinger-Poisson solver and, for our purposes, is configured to solve the single-band approximation with effective mass.It does this self-consistently with Poisson's equation.The software generates the first 20 eigenvalue solutions to the wave equation.The software allows specification for the data to be output in text files-used to collect data from the simulation.
The simulation output data collected includes wavefunction data, ψ 2 i , energy levels, E i , dipole matrix elements, |z ij |, position and subband dependent effective mass, m * , and LO phonon energies, E LO .Relevant energy levels, ones with the highest probability, are determined for each well.Once the E inj , E 3 , E 2 , and E 1 levels are found, relevant electron lifetimes, τ 3 , τ 2 and τ 32 are imported from the simulation data.The dipole matrix elements are used to calculate the oscillator strength, f osc,ij of each possible energy transition.The oscillator strength equation is given as where m 0 is the free-space electron mass and h is the modified Planck constant.

Fitness Evaluation and Candidate Selection
The algorithm is designed to produce diagonal-transition three-well active region designs with a five-well/barrier pair superlattice injector region.Current iterations of this algorithm rank candidates based on the following objective criteria: 3-5 µm emission and E 3 confinement, Equation (2); LO phonon resonance between E 2 and E 1 and alignment of E 1 and injector well, Equation (3); injector energy levels alignment/coupling, Equation (4); oscillator strength of the E 32 transition and a population inversion/gain metric, Equation ( 5).
where E ph is the expected photon energy, E c,max and E c,min are the maximum and minimum conduction band energies in the active region, respectively, ∆E inj is the thickness in eV of the injector miniband, and ∆N is the change in electrons-a measure of carrier gain.Parents are selected via rank selection.Each candidate is ranked based on the fitness of each objective.The top two of each are selected for carryover into the next generation.Additionally, top candidates from the aggregate ranks (mean of two objective fitness values) are selected for carry over.Finally, the top two from the complete aggregate fitness values are selected as the final carryover population.This is a very elitist approach, which has the potential to find solutions very quickly, but at the risk of being stuck in a locally optimized regime.To mitigate this, random genes are inserted into the population by mating (crossover) the top candidates.The rest of the population is generated by breeding (crossover and mutation) the carryover population.

Fitness Objective Calculation
The desired traits of candidates are given in Equations ( 2)-( 5).Normalized fitness evaluation is performed as given in Equations ( 6)-( 14).Equation ( 6) evaluates if transition energy, E 32 , is in the desired range.However, the E 3 energy level must be adequately contained in the first active region well: evaluated by Equation (7).Equation ( 8) evaluates if the E 21 transition is near the LO phonon scattering energy.It is possible, due to the lowering of the well with indium incorporation, that the E 1 level will lie below the conduction band minimum of the adjacent injector well.This undesirable condition is mitigated using Equation ( 9).In the injector region, electrons should transition from well to well via energy relaxation.The algorithm quantifies this transition by evaluating the oscillator strength of adjacent energy levels via Equation (10).Additionally, the last injector energy level should align with the E 3 level of the active region as evaluated by Equation (11).In this equation, evpdf refers to the extreme value probability density function, also known as the Gumbel distribution.Ideally the injector energy levels should form a mini-band with minimum energy difference; we use Equation (12) for this.Finally, a measure of the potential gain characteristics of the candidate device is evaluated by looking at carrier lifetimes, Equation (13), and the strength of the E 32 transition, Equation (14).
f 3b (∆E 3,inj ) = evpd f (∆E 3,inj , 0.01, 0.03) (11) Each of the fitness functions are combined into a set of four (4) aggregate fitness measures, Equations ( 15)- (18).These fitness evaluations are used to rank (maximize fitness value) and select the top candidates for carryover and breeding to generate the next generation:

Results and Discussion
To evaluate the quality of the fitness functions, they were tested against a reference design [27] from the training population, Figure 4.The resulting normalized fitness values are given in Table 3.The lower score for the injector coupling objective is indicative of the algorithm not considering all the possible transitions or tunneling, but rather just a high coupling factor to the highest probable wavefunction in the adjacent well.Per the sum rule, the oscillator strength coupling to all energy levels will be one.This particular result indicates that further study is needed on how to evaluate an effective relaxation/injector region.Currently, injectors for long wavelength devices are designed as a graded super lattice structure so as to produce a miniband.This approach produces many energy levels to which the electron may scatter.Since the gain metric is normalized using a Sigmoid function, Equations ( 13) and ( 14), the resulting fitness will have a maximum value of 1.However, the range of possible values for the ∆N metric is unknown.Therefore, the Sigmoid function parameters were chosen based on the range of values found in the reference designs.The goal of the algorithm is to maximize these values.Of note in this design is the fitness values for the injector coupling metric.The injector in this design does contain a miniband, but the energy levels are above the E 1 level, which for the coupling factors necessitates energy absorption-a factor this algorithm tries to reduce.

Average Fitness Values over Time
Figures 5 and 6 show the average fitness values for each objective over all generations, with an initial training population and with randomly generated initial candidates, respectively.With the initial population derived from a training set, it is clear that when compared to the random set, the initial average fitness values are higher, as expected.For the randomly generated population, average fitness values initially increase.However, the fitness values decrease then level after around 30-40 generations.With the insertion of random genes into the population, the lowering of fitness is not unexpected, but the large amount of elitism present in this algorithm should prevent the search from straying significantly.The randomly generated initial population is an interesting case study to see if by optimizing the objectives, the algorithm would produce recognizable or plausible devices.For the LO phonon and injector coupling based objectives, the algorithm shows a slow average overall improvement for around 70 generations, after which the improvement seems to level off.This is likely an indication of the random nature of inserting new genes into the population.Given many more generations in both the initial and random populations, it would be interesting to see if the algorithm would improve average fitness values.The jagged but relatively level average values of fitness are indicative of both the random gene insertion and the elitism present in the algorithm with no advanced methods of changing the selection mechanism over time.To improve, the algorithm needs to adapt a dynamic method of parental selection more advanced than simple rank selection.The multi-objective nature of QCL design would be a good candidate for divide and conquer methods where individual objectives are optimized in a separate process and design traits are shared between processes.

Bandstructure Diagrams of Top Candidates
QCL candidate band-structure graphs with corresponding wavefunctions are given for the top candidates in each objective ranking, as well as aggregate rankings.These graphs serve as representations of what is considered the 'top' output of the algorithm; however, the practicality of these candidates is discussed.Additionally, some commentary is given on the performance of a given ranking scheme and how additional ranking metrics could improve the outcome of the algorithms.

Initial Training Population
Given an initial training population, the algorithm is expected to be able to produce plausible optimized devices.Not all of the initial training population has optimized values for our target wavelength.However, they contain elements from working designs.Figure 7 shows the candidate with the highest combined fitness after 100 generations of a simulation run.Table 4 shows the normalized fitness values of this design candidate.The design, especially the active region, has traditional elements with a nicely aligned injector level and well confined E 3 level.The E 21 transition is on the order of an LO phonon and the target wavelength is achieved.However, as indicated by the lower fitness value, the injector region is not optimal.There is not a clear miniband being formed, and the third and fifth injector barriers are very wide reducing tunneling probability (a factor not currently considered by the algorithm).As will be shown shortly, modifications to the injector optimization algorithm are necessary to account for barrier thickness-using the oscillator strength alone does not seem to be enough.
Figure 8 shows the candidate with the best fitness value for the wavelength objective (1) after 100 generations of a simulation.One of the reasons for looking at the strained Al x Ga 1−x As/In y Ga 1−y As material system is the potential for increased conduction band offset with increasing indium incorporation.This would have the benefit of giving a large E 32 energy as well as good E 3 confinement.The candidate has alloy values of x = 0.22 and y = 0.25.Here, the algorithm has utilized the advantages of indium incorporation, but it is surprising that the aluminum percentage is lower.Given the built-in field due to the Piezo-electric properties of the (111) orientation, it is likely not necessary to have such a large percentage of aluminum.However, the algorithm does not take into account any negative effects of the large strain, nor the difficulty of pseudomorphic growth of highly strained layers.
Figure 9 shows the candidate with the best fitness value for the LO phonon objective (2).It is expected that this portion of the algorithm would be affected by the second and third active region well widths.The only challenge here is the accounting for the potential of the energy level to be too low for effective extraction into the injector region, which is affected by the well alloy ratio.The algorithm does not account for coupling strength, via oscillator or other mechanism, between the levels.This is likely why the algorithm produced a candidate with large barriers between the wells.Improvements should include a method for accounting for the coupling of wells and energy levels and/or tunneling probability to make for more realistic device design.One of the more difficult implementations in this algorithm was the optimization of the injector region.Figure 10 shows the candidate with the best fitness value for the injector region objective (3).The design shows a tight miniband of energy levels.However, as previously noticed in the overall best design, the injector region contains a large barrier.Again, this should be accounted for to improve the design.Optimization of the gain function alone tends to produce devices with many of the characteristics expected of a working QCL. Figure 11 shows the candidate with the best fitness value for the gain objective (4).The gain metric accounts for the positioning of the three active region levels; therefore, it is not surprising to see such good results from only this objective.However, this metric alone does not account for wavelength of the device.The candidate shown below has an expected wavelength of 9.66 µm.Therefore, this metric cannot be used alone.In addition, there is no direct optimization of the injector; however, the oscillator strength of the transition between E 3 and E 2 is being considered.It is interesting to see that the gain metric does not seem to be affected by lack of a minigap in this design.

Random Initial Population
To further test the ability of the algorithm design plausible QCL structures, we removed the training population by producing a randomly generated initial population.If the objectives and fitness measures are sufficient, given enough time, the algorithm should produce workable devices.Figure 12 shows two candidates with high aggregate fitness values.While these candidates show some elements of known working QCLs, they have other features that are not traditionally characteristic of a working device.The candidate shown on the left has very large barriers in the active region, effectively uncoupling the wells and has large wells in the injector region.The candidate on the right shows a more traditional coupled well active region, but the target wavelength is missed.Additionally, the injector region contains both large barriers and wells.Improvement of the algorithm is necessary both in the parental selection and the fitness metrics being used.However, it is remarkable that after just 100 generations given no training set, the results do show some promising signs.Figure 13 shows top candidates for each of the four objectives after 100 generations of a simulation run with only a random initial population.
These candidates, like the ones produced with a training population, show similar elements, but also similar weaknesses.Objectives (1) and (2) produce candidates with large barriers to decouple the wells, objective (3) produces large barriers, and objective 4 alone produces traditionally looking active regions but does not account for the desired wavelength.
Adding objectives/fitness metrics to the scope of the algorithm may produce better results (working devices) but necessitates more sophisticated algorithms to assure speedy convergence on optimal solutions.Our results here show how the different objectives led to specific results and when combined can produce plausible devices.There is, however, room for improvement.

Conclusions
Our past work explored the possibility of quantum cascade laser designs based on the non-traditional GaAs(111) surface.We found that growing pseudomorphic InAs on GaAs(111)B is possible [4], which opens up a new material system on which to develop devices.Ternary In x Ga 1−x As at reasonable alloy content percentages benefit from lower threshold voltages and higher efficiency while extending the operating range [3].
To explore this concept, a multi-objective evolutionary algorithm was developed and used to search the design space made possible by the promise of this new material system.The algorithm was successful in searching the design space, but was ultimately limited in many of its conclusions.The algorithm found designs, which it ranked high in all or most of the objective categories, but would be unlikely to work in real life.Further enhancement of the algorithm is possible by adding more objectives that will regulate the fitness of the candidate solutions more effectively.Future work should at least incorporate objectives for mini-band and mini-gap formation.Additionally, the interplay between objectives in this multi-objective approach should be considered carefully to improve the resultant designs.Ultimately, there is hope that, with further refinement, the algorithm would be successful in finding working devices more effectively than traditional design approaches, such as the design toolbox method of [1].This algorithm uses rank selection for its parental selection and mutation operations for enhancement of multimodal search range.Future algorithm development should expand the use of more advanced parental selection methods to provide quicker and more effective search.More advanced evolutionary algorithm techniques may be in order to help narrow the variable range.This could be accomplished using sequential objective evaluation approach, where one objective is optimized and used to narrow the variable parameter range, then another objective, and so on.There are many potential approaches, but, overall, the results of this work show promising possibilities for the viability of an intelligent search approach for solving the problem of designing device structures.

Figure 1 .
Figure 1.General flow diagram of the evolutionary algorithm used in this study.

Figure 2 .
Figure 2. Generic gene sequence with value ranges.

Figure 3 .
Figure 3. Schematic flow for the initial population.

Figure 4 .
Figure 4. Al x Ga 1−x As/In y Ga 1−y As QC laser reference design, adopted from [27] with alloy proportions changed to x = 0.25 and y = 0.1 and simulated on GaAs (111).

Figure 5 .
Figure 5. Plots of normalized average fitness per generation for (a) wavelength, (b) LO phonon, (c) injector coupling, and (d) gain objectives when initial population was generated from a training set.

Figure 6 .
Figure 6.Plots of normalized average fitness per generation for (a) wavelength, (b) LO phonon, (c) injector coupling, and (d) gain objectives when initial population was generated randomly.

Figure 8 .
Figure 8. Al x Ga 1−x As/In y Ga 1−y As QC laser candidate simulated on GaAs (111) with top objective (1) fitness rank given a training initial population.Relevant wave functions, E 3 and E 2 are in bold and colored.

Figure 9 .
Figure 9. Al x Ga 1−x As/In y Ga 1−y As QC laser candidate simulated on GaAs (111) with top objective (2) fitness rank given a training initial population.Relevant wave functions, E 2 and E 1 are in bold and colored.

Figure 10 .
Figure 10.Al x Ga 1−x As/In y Ga 1−y As QC laser candidate simulated on GaAs (111) with top objective (3) fitness rank given a training initial population.Relevant wave functions, E inj,ij are in bold and colored.

Figure 11 .
Figure 11.Al x Ga 1−x As/In y Ga 1−y As QC laser candidate simulated on GaAs (111) with top objective (4) fitness rank given a training initial population.

Figure 12 .
Figure 12.Al x Ga 1−x As/In y Ga 1−y As QC laser candidates simulated on GaAs (111).(a) candidate with top aggregate fitness rank and (b) candidate showing diversity of designs.

Figure 13 .
Figure 13.Al x Ga 1−x As/In y Ga 1−y As QC laser candidates simulated on GaAs (111) with top fitness rank for (a) wavelength, (b) LO phonon, (c) injector coupling, and (d) gain objectives given a random initial population.

Table 1 .
General simulation parameters for the multi-objective evolutionary algorithm (MOEA).

Table 2 .
Application specific simulation parameters for the MOEA.

Table 3 .
Reference design fitness evaluation.

Table 4 .
Initial population best candidate design fitness evaluation.