Laser TCAD-level simulation-oriented longitudinal mode self-adaptive selection method

Abstract

The invention discloses a laser TCAD level simulation-oriented self-adaptive longitudinal mode selection method, which specifically comprises the following steps of: effectively extracting photoelectric characteristics of any given structure by utilizing a laser TCAD level photoelectric simulation technology; comprising a transverse optical mode, ABC parameters in a carrier rate equation, a gain spectrum in a parametric analytic expression form, a spontaneous radiation spectrum and refractive index change; and then a special carrier rate equation extracted by using an optical traveling wave model and laser TCAD level electrical simulation is used for joint solving. Starting from an actual device structure and a physical mechanism, compared with other methods, the method does not require a large number of introduced external parameters, and has a great application value.

Description

Technical Field

[0001] The present invention belongs to the field of TCAD-level numerical simulation of lasers, and particularly relates to a method for adaptively selecting longitudinal modes for TCAD-level simulation of lasers. Background Art

[0002] As a key component of the Internet and other optical communication systems, lasers have received great attention throughout the generations. Due to the complexity of the internal physical mechanisms of optoelectronic devices, device designers and researchers highly rely on empirical knowledge obtained from a large number of trial-and-error wafer runs; computer-aided design software TCAD with photon device simulation technology as the core provides another option with lower cost, that is, performing computer numerical simulation of the device before actual device manufacturing to reduce wafer run trial-and-error. The characteristic of such software is that its algorithm starts from the physical model of the device and can provide the device designer with the field information inside the device to guide device design; while general device users do not need to use this software and only need to confirm the device characteristics through the behavioral model.

[0003] The specific challenges in laser simulation are that it is necessary to separate electronics and photonics as much as possible to simplify the problem but also consider the complex interaction between electrons and light; and it also faces the dilemma of selecting theoretical models caused by the diversity of materials, devices, physical mechanisms, and modeling methods. Current research leaves this difficulty to the users themselves to handle [1]: when using simulation software, users not only need to set the basic laser structure, but when dealing with optoelectronic coupling, engineers need to personally and repeatedly debug various parameters extremely tediously (even constructing the initial solution of a complex structure by themselves using simple structure simulation) to obtain a relatively correct result.

[0004] Traditional laser TCAD simulation tools need to separately process the optical model and electrical model of the device, and then use a suitable method to couple these two models. The electrical characteristics of the device are generally described by diffusion drift (abbreviation: DD), which is a set of combined equations, including the Poisson equation:

[0005]

[0006] Among them, ε0 is the vacuum permittivity, and ε dc is the relative permittivity (microwave frequency); q is the unit charge amount; V, n, and p are the electric potential, electron, and hole concentrations. And the continuity (transport) equations of electrons and holes:

[0007]

[0008] N D 、N A 、N tj are the concentrations of shallow-level donors, acceptors, and deep-level traps respectively. f D, f A and f tj are the corresponding energy level occupancies. and are recombination terms related to deep level traps (electrons and holes). J n and J p are the electron density and the hole current density respectively, which can be given by drift and diffusion terms. R au is the Auger recombination term, G imp is the impact ionization term as a reciprocal process. To simulate a laser, a photo - electric coupling term G opt (t), R sp and R st need to be introduced into the continuity equations of electrons and holes. They correspond to the photo - generated carrier generation term, the stimulated recombination term and the spontaneous emission term respectively. During the laser simulation process, they represent the mechanism of optical loss caused by material absorption, the mechanism of converting electrical energy into optical power and the noise source respectively.

[0009] If a quantum region is involved as the active region of the material, the (multi -) quantum well will be regarded as a separate confinement heterojunction SCH. Besides the above equations, since the actual carriers in the quantum well are cold carriers in sub - energy levels constrained by the Schrödinger equation, a cold carrier continuity (transport) equation needs to be added. The original equations (2) and (3) become the transport equations of hot carriers. Their carrier transport equations need to add carrier heating and cooling terms to express the conversion between hot and cold carriers. The Poisson equation (1) also needs to be modified to reflect the contributions of hot and cold carriers. The additional cold carrier electrical equations are:

[0010]

[0011] In the formula, the cold hole concentration and the cold electron concentration in the direction perpendicular to the quantum well for transport phenomena are the cold hole heating time, is the hot hole cooling time, and the relaxation time of electrons is the same. The parameter is the non - radiative recombination term, which includes R au and trap - assisted recombination (SRH recombination).

[0012] The optical model of the device usually uses a spatially and temporally separated standing light field model, that is, a method based on the photon rate equation. For more convenient description, starting from the vector wave equation of light:

[0013]

[0014] The meanings of each parameter are: the relative magnetic permeability μ of the devicer (r, t), relative permittivity (function) ε r (r, t, τ), speed of light c, electric field strength E. The asterisk * represents the convolution within the time τ. The role of convolution is to take into account the memory effect of the material's polarizability and consider the influence of previous excitations on all responses. From the perspective of optoelectronic coupling, the spontaneous emission term in the previous electrical model enters the vector wave equation through the source term F, and the optical gain and refractive index change terms enter the vector wave equation through the relative permittivity (function).

[0015] The wave equation for the coupled electrical problem is basically unsolvable. At this time, it is necessary to separate the time variable and the spatial variable of the vector wave equation. This idea of separating variables can be achieved through mode expansion:

[0016]

[0017] In the above equation, the complex-valued scalar coefficient a μ (t), real-valued scalar coefficient ω' μ (τ) and eigenvector Ψ μ (r, t) are introduced. The optical field is decomposed into a given number of modes, constituting N eigenvalue problems:

[0018]

[0019] In the equation, the eigenvector Ψ μ (r, t) is the optical mode shape. Generally, the transient mode is not included in the technical solutions of laser TCAD simulation. The time-dependent part of the wave equation after expansion by modes is the rate equation (i.e., the photon rate equation) for describing the time evolution of the average electromagnetic energy S μ (t):

[0020]

[0021] is the electromagnetic power of spontaneous emission coupled to the μ-th mode, ω' μ (t) is the angular velocity of the μ-th mode, ω” μ is the mode change rate of the electromagnetic field of the μ-th mode.

[0022] As mentioned before, the transient mode is not included in the technical solutions of TCAD simulation. Although this method can most simply introduce the meaning of the photon rate equation, the mode Ψ μ (r, t) we hope to obtain is only related to space. Therefore, it is assumed that the light wave in Equation (4) is composed of harmonics with slowly varying envelopes and discrete frequencies:

[0023]

[0024] In the equation, uk and ω k is the amplitude and angular frequency of the kth harmonic. By further assuming the linear characteristics of Equation 4 (i.e., the polarizability has no explicit dependence on the optical electric field intensity) and the polarizability in the device frequency domain It changes much faster in time than the slowly varying envelope function in the frequency domain, and the wave equation 4 should be simplified in the time domain using Fourier transform:

[0025]

[0026] Where ω0 is the reference frequency and c is the speed of light. This equation is for the harmonic frequency ω k＝0 The slowly varying envelope u of the light field k＝0 The constraint equation, J sp It is a spontaneous emission current source. Since there is no quantized light field, the spontaneous emission current source can only be evaluated phenomenologically through the electrical equations (1-4) [2]. is the polarizability at the reference frequency ω0. This variable actually contains implicit time dependence because the applied bias will cause changes in the electrical equations (1 to 3). If the bias V is a time-dependent function, then the dependence of the electrical equations It will also become a function of time and require further separation There are two parts: time-related and time-independent:

[0027]

[0028] Equation 11 can be rewritten by substituting 10 into:

[0029]

[0030] At this time, the opponent 12 adopts a similar pattern to equation 5 to truly separate the space-time variables, and the resulting harmonic u k＝0 The three-dimensional optical mode v l (r) is only the eigenvalue problem of bias V and spatial variables without time variables:

[0031]

[0032] where β is the three-dimensional mode propagation constant of the device at a certain bias voltage V, the integer 1 represents the first eigenvalue and eigenfunction, L is the entire set of eigenvalues, and the harmonic u k＝0 The expression can be written as:

[0033]

[0034] Since the transverse structure of a laser is generally a waveguide that is very good at the operating wavelength, the three-dimensional optical mode of the entire device can be further decomposed into a transverse mode and a longitudinal mode along the light propagation direction. Taking the waveguide structure of an edge-emitting laser as an example, define the z direction as the light transmission direction, the transverse mode is in the xy plane, and define the propagation factor exp(jβ m z) = exp(jk m n eff z), where n eff is the complex effective refractive index of the transverse mode φ m (x, y), and β m and k m are the propagation constant and wave number of the transverse mode φ m (x, y), respectively. Further assume that the device varies very little along the z-axis, and a three-dimensional optical mode v k＝0 of the harmonic u l (r) can be written as:

[0035] v l (r) = xE n (z; φ m (x, y))φ m (x, y)exp(jβ m z) (15)

[0036] where φ m (x, y) is the m-th transverse mode, E n (z; φ m (x, y)) is the n-th longitudinal mode that depends on the specific transverse mode φ m (x, y), and x represents the polarization direction.

[0037] When the waveguide size of an edge-emitting semiconductor laser is equal to or greater than the wavelength, assuming that the single transverse mode and the transverse mode field distribution are independent of frequency within the lasing spectrum, the transverse optical mode φ0 of the laser can be described by the scalar wave equation as follows:

[0038]

[0039] It should be noted that when multiple transverse modes are transmitted in the device, a more complex vector Helmholtz equation needs to be solved.

[0040] The longitudinal mode E n (z; φ m (x, y)) then satisfies an eigenvalue problem determined by the operating state of the device and the device end face conditions under a certain transverse mode φ m＝0 . For simplicity, the expression is directly given here:

[0041]

[0042] where light waves propagate forward f and backward b along z to form a standing wave, and E n (z) is further separated into the longitudinal distribution e f (z) of the forward wave and the longitudinal distribution e b (z) of the backward wave. β is the three-dimensional mode propagation constant of the device under a certain bias voltage V, and g V (z, ω0) is the modal gain (at the reference frequency ω0) of the definite transverse mode φ0 at point z under a certain bias voltage V. Δn V (z, ω0) is the value at the reference frequency ω0 of the change in the refractive index over the entire frequency domain due to the gain / carrier injection process. Γ is the optical field confinement factor, which reflects the degree to which the transverse optical mode of the device is confined in the active region. α eff is the optical modal loss excluding interband absorption. The end face conditions are as follows:

[0043]

[0044] where R r (ω0) is the reflectivity of the right end of the device (at the reference frequency ω0), and R l (ω0) is the reflectivity spectrum of the left end of the device (at the reference frequency ω0).

[0045] The three-dimensional mode of the harmonic u k＝0 should be modified to the following form to reflect the characteristics of the longitudinal mode:

[0046]

[0047] where β m is the propagation constant of a certain definite transverse mode φ m (x, y).

[0048] Assuming single transverse mode, the complete expression of a single harmonic is:

[0049]

[0050] In the laser simulation at the CAD level, the actual photon rate equation used is not Equation 7, but U n (t) in Equation 21 or U l (t) in Equation 13. These two equations are essentially the same. Generally, the photon rate equation of the three-dimensional mode of a certain harmonic is referred to as the photon rate equation of a certain longitudinal mode because the longitudinal mode is solved based on a certain definite transverse mode. The complete form of the photon rate equation is given here without derivation:

[0051]

[0052] where m ∈ L, neff The effective refractive index and wave number of the transverse mode corresponding to the three-dimensional mode of k0, ν g is the group velocity ν g ≡c / ng≈c / n eff , where ng is the group velocity refractive index and n eff ≈ng. Δg Vlm (t) is the cross-coupling coefficient (under bias V) of the three-dimensional modes numbered 1 and m, reflecting the degree of competition mode gain between the two. Δg Vmm (t) is the self-coupling coefficient (under bias V) of the three-dimensional mode, mainly reflecting the gain of a single three-dimensional mode. is the contribution of inhomogeneous spontaneous emission to the m-th three-dimensional mode.

[0053] Many details are omitted here, but the general model for laser simulation is already available. Next, it will be combined with the appendix Figure 1 as a technical background description.

[0054] Appendix Figure 1 The (002) technique in it is a method for solving the eigenvalue problem [3], used to solve the eigenvalue problems of optical transverse and longitudinal modes.

[0055] Appendix Figure 1 The (003) technique in it is to expand the electrical equations 1-4 using discrete methods such as finite volume (for the purpose of improving the solution convergence), and then add the scalar photon rate equation 21 of each longitudinal optical mode and solve it by strong coupling using the Newton method [4]. The flow chart shown in Appendix Figure 1 also includes the laser simulation scheme at the TCAD level: at a certain bias, iteratively and step by step, use the (002) technique to calculate the transverse / longitudinal optical mode problems and use the (003) technique to calculate the laser output problems related to electricity until convergence, rather than solving all the above equations simultaneously in one calculation. The reason for adopting this iterative method is that on the physical model, the photon rate equation and the electrical equations 1-3 are strongly coupled, while the mode distribution is weakly coupled to the laser output: the average electromagnetic energy U m (t) directly comes from the optical material gain and spontaneous emission terms contained in the electrical equations 1-4, and near the lasing threshold of the laser, even a small bias change applied to the device will cause a huge change in U m (t); in contrast, even if the refractive index changes due to the gain or carriers inside the device, techniques such as gradually increasing the bias as shown in Appendix Figure 1 can always be tried to ensure that the optical modes do not change significantly in one iteration process.

[0056] Appendix Figure 1The biggest problem faced by the laser simulation scheme in the shown flowchart comes from the physical model itself. Since there is energy leakage in the laser, the total energy is not conserved within the cavity of the device. At least the traveling wave of the optical field along the z-direction cannot be expressed as the sum of a finite number of standing waves, that is, the harmonic wave composed of an infinite number of standing wave sums reflected in Equation 20. However, due to the requirement of finite numbers in numerical calculations [5], truncation will definitely be introduced into Equation 20 (which will lead to loss of precision), resulting in a well-known difficult problem of how to select longitudinal modes.

[0057] In the Figure 1 shown technical process, the initial guess of the longitudinal mode in Technique (001) is very important: the longitudinal mode distribution must be a problem coupled with the electrical equation. However, the longitudinal mode distribution is required before solving the electrical equation, and solving the longitudinal mode distribution requires solving the electrical equation. Therefore, it is necessary to find an initial guess of the longitudinal mode distribution to initialize the entire laser simulation scheme. Due to the characteristics of the laser, a small number of longitudinal modes will definitely have the vast majority of the power. It is very important to retain relatively accurately the longitudinal modes that may lase during the simulation process. Currently, the main methods for selecting longitudinal modes in research are the cold cavity model [6], the floating bias model [7], and the Green's function method [8]. The core idea of the cold cavity model is to retain as many longitudinal modes as possible to reduce the number of times of solving longitudinal modes in the entire simulation process, while the core idea of the floating bias model is to only retain the most important modes and frequently calculate Equation 21 to ensure that no important modes are missed.

[0058] References:

[0059] [1] US201314025902 [P]

[0060] [2] M. Lax and W. Louisell, "Quantum noise IX: Quantum Fokker - Planck solution for laser noise," in IEEE Journal of Quantum Electronics, vol. 3, no. 2, pp. 47 - 58, February 1967.

[0061] [3] THE DAVIDSON TYPE METHOD AND ITS IMPLEMENTATION. Journal on Numerical Methods and Computer Applications, 2006, 27(3): 218 - 240

[0062] [4]R.E.Bank, D.J. Rose and W.Fichtner, "Numerical methods for semiconductor device simulation," in IEEE Transactions on Electron Devices, vol. 30, no. 9, pp. 1031-1041, Sept. 1983.

[0063] [5]Brüning, E., & Blanchard, P. (1992). Variational Methods in Mathematical Physics: A Unified Approach.

[0064] [6]Y. Xi, X. Li and W.-P. Huang, "Time-Domain Standing-Wave Approach Based on Cold Cavity Modes for Simulation of DFB Lasers," in IEEE Journal of Quantum Electronics, vol. 44, no. 10, pp. 931-937, Oct. 2008.

[0065] [7]Xun Li, A.D. Sadovnikov, W..-P. Huang and T. Makino, "A physics-based three-dimensional model for distributed feedback laser diodes," in IEEE Journal of Quantum Electronics, vol. 34, no. 9, pp. 1545-1553, Sept. 1998.

[0066] [8]B. Tromborg, H. Olesen and X. Pan, "Theory of linewidth for multielectrode laser diodes with spatially distributed noise sources," in IEEE Journal of Quantum Electronics, vol. 27, no. 2, pp. 178-192, Feb. 1991. Summary of the Invention

[0067] In view of the above problems, the present invention provides a method for adaptively selecting longitudinal modes for laser TCAD-level simulation.

[0068] A method for adaptively selecting longitudinal modes for laser TCAD-level simulation according to the present invention effectively extracts the optoelectronic characteristics of any given structure by using optoelectronic simulation technology at the laser TCAD level, including transverse optical modes, ABC parameters in the carrier rate equation, gain spectra in parametric analytical form, spontaneous emission spectra, and refractive index changes; subsequently, a combined solution is carried out by using an optical traveling wave model and a special carrier rate equation extracted from laser TCAD-level electrical simulation; specifically, it includes the following steps:

[0069] Step 1: Conduct electrical simulation on the laser structure by solving the diffusion drift Poisson equations at each bias voltage V, and extract the analytical form of the gain spectrum.

[0070] Step 2: Solve the transverse mode of the laser, and calculate the modal spontaneous emission spectrum and modal gain spectrum.

[0071] Step 3: Extract a phenomenological carrier rate equation model from the results obtained by electrical equation simulation.

[0072] Step 4: Solve the coupled equations of the carrier rate equation and the traveling wave optical model.

[0073] Further, step 1 is specifically as follows:

[0074] Step 1.1: Assume that when the bias is below the lasing threshold, the photon density term U m (t) in the photon rate equation is 0. Since the photon density is small, the stimulated emission term is also small, and the stimulated emission term R in the photon rate equation and the electrical equation is deleted. st ; Study a slice along the z-axis, and the electrical equation to be solved is two-dimensional or one-dimensional.

[0075] Step 1.2: Expand the equation spatially by the finite volume method and then solve it by the Newton method, and sequentially solve by scanning the current bias from 0 to a large current bias until the gain spectrum becomes transparent (just neither absorbing nor amplifying) at a certain frequency, and then continue to scan to 100 times the current magnitude of this transparent current bias to avoid missing the gain at a specific frequency; this step requires the use of a gain model, and all the parameters required by the gain model can be obtained by solving the electrical equation.

[0076] Step 1.3: According to the specific type of active region, different gain models are adopted for bulk material active regions and multiple quantum well active regions:

[0077] A. For the gain model at a single r position (in k space) of the bulk material active region, it is:

[0078]

[0079] In the formula, q is the electron charge, h is the Planck constant, k is the electron wave number, is the square of the free electron mass, ε0 is the permittivity of vacuum, n τ is the refractive index of the optical frequency, c is the speed of light, v is the photon frequency, E ij = hv; f i and f j is E i and E j are the Fermi-Dirac distribution functions of the energy levels, D r (E ij ) is the joint density of states between two transition energy levels, M is the momentum (transition) matrix. Under the condition of momentum conservation, the gain at all possible ij energy levels needs to be integrated and summed in the momentum space to obtain the gain spectrum g(r, hv) at a single point r in space.

[0080] B. For the quantum well active region, the gain model needs to correct the joint density of states D r (E ij ) to be the two-dimensional density of states:

[0081]

[0082] In the formula, m r is the reduced effective mass, d z is the quantum well width, is the reduced Planck constant.

[0083] Step 1.4: Adopt a unified spontaneous emission model for the active region and the passive region. The general form of the spontaneous emission model r sp (located in the k space) is:

[0084]

[0085] In the formula, D opt is the photon density of states. Similarly, r sp is also transformed into r sp (r, hv) in the same way as the gain spectrum.

[0086] Step 1.5: Calculate the gain, spontaneous emission, and refractive index change by electrical simulation from the momentum space through a rigorous physical model to parameterize the stimulated emission term and the refractive index change term caused by carrier injection:

[0087]

[0088] In the formula, is the gain spectrum in the analytical expression form, is the photon energy, N is the cold carrier concentration in the active region. It is assumed that the electron-hole concentrations in the active region are equal and stimulated emission only occurs in the active region. Other parameters are fitting parameters; is the refractive index change in an analytical expression form. Other parameters are fitting parameters except for the photon energy and the cold carrier concentration in the active region.

[0089] Furthermore, step 2 is specifically as follows:

[0090] Step 2.1: Ignore the refractive index change caused by gain and carrier injection, and use the cold cavity model to solve the optical transverse mode. The frequency scanning range is set to the gain spectrum range obtained from step 1.5.

[0091] Step 2.2: Normalize the parameters obtained from the electrical solution in step 1.2 for each bias V to modal parameters according to the transverse mode:

[0092]

[0093] In the formula, Σ is the xy-plane regions along the z-axis where the transverse modes of the entire device are located, and Σ ar is the area of the active region. By evaluating the modal parameters at a specific point z0, φ 2 (x, y) is the optical transverse mode distribution, Γ is the mode confinement factor, and α eff (z) is the modal loss, is the polarizability in the frequency domain, Δn(z, ω0) is the refractive index change, and g(z, ω0) is the gain spectrum.

[0094] Furthermore, step 3 is specifically as follows:

[0095] Step 3.1: Starting from the strict spontaneous emission term used in the actual electrical equation and the strict expression in the steady state is:

[0096]

[0097] In the formula, E g is the bandgap between the ij energy bands, and its empirical expression is:

[0098]

[0099] In the formula, N and P are the electron and hole concentrations in the active region, N0 and P0 are the electron and hole concentrations in the active region at the equilibrium state, and Br is the spontaneous recombination coefficient.

[0100] Step 3.2: Extract the Br spontaneous recombination coefficient through the weighted average method. Here, it is first emphasized that is an implicit function of the electron N and hole concentration P in the active region, and what is used in the solution is Spatial discretization form of variables

[0101]

[0102] In the formula, the part with [m] refers to the m-th discrete variable of the spatial discretization form of the variables in the active region, and M is the total number of discrete points in the active region.

[0103] Step 3.3: Extract the steady-state non-radiative recombination terms actually used Coefficients A and C of, and divide into Auger recombination term and SRH recombination term

[0104]

[0105] In the formula, A e , A h are the SRH decay rates of electrons and holes respectively, τ e , τ h are the SRH time constants of electrons and holes respectively, C e , C h are the Auger recombination coefficients of electrons and holes respectively.

[0106] Step 3.4: Obtain the expressions of coefficients A and C by weighted average:

[0107]

[0108] Step 3.5: Introduce the simplified carrier rate equation established only along the z-axis in the electrical equations solved by coupling the diffusion-drift-Poisson photon rate equation with the Newton method:

[0109]

[0110] In the formula, N(z,t) is the cold carrier concentration in the active region without distinguishing electrons and holes, is the slowly varying envelope of the optical electric field intensity along the z-axis, J(z,t) is the injection current density, η is the injection efficiency, n eff is the effective refractive index of the transverse mode, Γ is the mode field confinement factor, d is the thickness of the active region, ω0 is the reference frequency, Σ ar is the area of the active region; g als (ω; N(z,t)) is the analytical expression of the gain extracted in Step 1.5.

[0111] The beneficial effects of the present invention compared with the prior art are:

[0112] First, it has good self - adaptability. Starting from the actual device structure and physical mechanism, the present invention does not require a large number of externally introduced parameters compared with other methods. For a specific structure, it can adaptively calculate its parameters to solve the longitudinal modes.

[0113] Second, it utilizes advanced TCAD technology compared with ordinary methods and actually constitutes a key component of the TCAD simulation of lasers. Combined with other components, it can form a complete commercial TCAD laser simulation software. BRIEF DESCRIPTION OF THE DRAWINGS

[0114] Figure 1 It is a flowchart of a typical TCAD - level laser simulation.

[0115] Figure 2 The transverse optical mode (TM mode) of a quantum cascade laser obtained by solving.

[0116] Figure 3 It is the energy band diagram obtained by using the (003) technology to solve the electrical structure of a quantum cascade laser.

[0117] Figure 4 It is the modal gain spectrum obtained by using the (003) technology to solve the electrical structure of a quantum cascade laser.

[0118] Figure 5 It is the modal spontaneous emission spectrum obtained by using the (003) technology to solve the electrical structure of a quantum cascade laser.

[0119] Figure 6 It is the spectrogram of the QCL longitudinal mode obtained by our technology. DETAILED DESCRIPTION OF THE INVENTION

[0120] The present invention will be further described in detail below with reference to the drawings and specific embodiments.

[0121] Step 1: Electrically simulate the laser structure by solving the drift - diffusion Poisson equations at each bias V, and extract the analytical form of the gain spectrum.

[0122] Step 1.1: Assume that when the bias is below the lasing threshold, the photon density term U m (t) in the photon rate equation is 0. Since the photon density is small, the stimulated emission term is also small, and the stimulated emission term R in the photon rate equation and the electrical equation is deleted st ; Study a slice along the z - axis, and the electrical equation to be solved is two - dimensional or one - dimensional.

[0123] Step 1.2: After spatially expanding the equation by the finite volume method, solve it by Newton's method, and sequentially solve to gradually scan the current bias from 0 to a large current bias until the gain spectrum becomes transparent (neither absorbing nor amplifying) at a certain frequency. Then continue to scan to 100 times the magnitude of this transparent current bias to avoid missing the gain at a specific frequency. This step requires the use of a gain model, and all the parameters required for the gain model can be obtained by solving the electrical equations.

[0124] Step 1.3: According to the specific type of active region, different gain models are adopted for the bulk active region and the multiple quantum well active region:

[0125] A. The gain model for a single r position (in k-space) in the bulk active region is:

[0126]

[0127] In the formula, q is the electron charge, h is the Planck constant, k is the electron wave number, is the square of the free electron mass, ε0 is the permittivity of free space, n τ is the refractive index of the optical frequency, c is the speed of light, v is the photon frequency, E ij = hv; f i and f j are the Fermi-Dirac distribution functions of the E i and E j energy levels, D r (E ij ) is the joint density of states between two transition energy levels, M is the momentum (transition) matrix. Under the condition of momentum conservation, the gain of all possible ij energy levels needs to be integrated and summed in momentum space to obtain the gain spectrum g(r, hv) at a single point r in space.

[0128] B. For the multiple quantum well active region, the gain model needs to correct the joint density of states D r (E ij ) to a two-dimensional density of states:

[0129]

[0130] In the formula, m r is the reduced effective mass, d z is the quantum well width, is the reduced Planck constant.

[0131] Step 1.4: Adopt a unified spontaneous emission model for the active region and the passive region. The general form of the spontaneous emission model r sp (in k-space) is:

[0132]

[0133] where D opt is the photon density of states, and r is the same as above sp is also transformed into r in the same way as the gain spectrum sp (r, hv).

[0134] Step 1.5: Calculate the gain, spontaneous emission, and refractive index change to parameterize the stimulated emission term and the refractive index change term caused by carrier injection from the momentum space through electrical simulation using a rigorous physical model:

[0135]

[0136] where is the gain spectrum in analytical expression form, is the photon energy, N is the cold carrier concentration in the active region. It has been assumed that the electron-hole concentrations in the active region are equal and stimulated emission only occurs in the active region. Other parameters are fitting parameters; is the refractive index change in analytical expression form. Other parameters are fitting parameters except for the photon energy and the cold carrier concentration in the active region.

[0137] Step 2: Solve the transverse mode of the laser and calculate the modal spontaneous emission spectrum and modal gain spectrum.

[0138] Step 2.1: Ignore the refractive index change caused by gain and carrier injection, and use the cold cavity model to solve the optical transverse mode. The frequency scanning range is set to the gain spectrum range obtained in Step 1.5.

[0139] Step 2.2: Normalize the parameters obtained from the electrical solution in Step 1.2 for each bias V to modal parameters according to the transverse mode:

[0140]

[0141] where Σ is the xy-plane regions along the z-axis where the transverse modes of the entire device are located, Σ ar is the area of the active region. By evaluating the modal parameters at a specific point z0, φ 2 (x, y) is the optical transverse mode distribution, Γ is the mode confinement factor, α eff (z) is the modal loss, is the polarizability in the frequency domain, Δn(z, ω0) is the refractive index change, and g(z, ω0) is the gain spectrum.

[0142] Step 3: Extract a phenomenological carrier rate equation model from the results obtained from the electrical equation simulation.

[0143] Step 3.1: Starting from the rigorous spontaneous emission term used in the actual electrical equation, the rigorous expression in the steady state is as follows:

[0144]

[0145] where E g is the band gap between the i - j energy bands, and its empirical expression is:

[0146]

[0147] where N and P are the electron and hole concentrations in the active region, N0 and P0 are the electron and hole concentrations at equilibrium in the active region, and Br is the spontaneous recombination coefficient.

[0148] Step 3.2: Extract the Br spontaneous recombination coefficient by the method of weighted average. Here, it is first emphasized that is an implicit function of the electron N and hole concentration P in the active region, and the spatial discrete form of the variables is used in the solution

[0149]

[0150] where the part with [m] refers to the m - th discrete variable of the spatial discrete form of the active region variables, and M is the total number of discrete points in the active region.

[0151] Step 3.3: Extract the A and C coefficients of the steady - state non - radiative recombination term actually used, and divide into Auger recombination term and SRH recombination term

[0152]

[0153] where A e and A h are the SRH decay rates of electrons and holes respectively, τ e and τ h are the SRH time constants of electrons and holes respectively, and C e and C h are the Auger recombination coefficients of electrons and holes respectively.

[0154] Step 3.4: Obtain the expressions of the A and C coefficients by weighted average:

[0155]

[0156] Step 3.5: Introduce the carrier rate equation simplified from the electrical equations solved by the coupled Newton method of the diffusion - drift - Poisson - photon rate equations, which is established only along the z - axis:

[0157]

[0158] In the formula, N(z, t) is the cold carrier concentration in the active region without distinguishing between electrons and holes, is the slowly varying envelope of the optical electric field intensity along the z-axis, J(z, t) is the injection current density, η is the injection efficiency, n eff is the effective refractive index of the transverse mode, Γ is the mode field confinement factor, d is the thickness of the active region, ω0 is the reference frequency, Σ ar is the area of the active region; g als (ω; N(z, t)) is the analytical expression of the gain extracted in step 1.5.

[0159] Step 4: Solve the coupled equations of the carrier rate equation and the traveling-wave optical model.

[0160] Example:

[0161] Take the solution process of the initial longitudinal mode of a quantum cascade laser (QCL) with a distributed feedback Bragg grating as an example. In a QCL, the active region material consists of a superlattice composed of engineering layers of different (usually two alternating) materials, where the structure emits light under an electrical bias. The region where optical transitions occur is called the "active region", and the transport layer connecting subsequent active regions is called the "injection region". When population inversion is achieved, the material becomes a positive gain material, which means it can amplify the passing light.

[0162] This example is an AlInAs / InGaAs / InP laser (λ = 4.7 μm), strain-compensated Al 0.638 In 0.362 As / In 0.669 Ga 0.331 As / InP laser structure, the active region consists of 35 periods. Designed for four-well two-phonon resonance. In nanometers, starting from the injection barrier, the layer sequence of one period of the structure is: 3.8, 1.2, 1.3, 4.3, 1.3, 3.8, 1.4, 3.6, 2.2, 2.8, 1.7, 2.5, 1.8, 2.2 , 1.9 , 2.1 , 2.1 , 2.0, 2.1, 1.8, 2.7, 1.8. The bold layers are AlInAs layers, and the total thickness of one period is 50.4 nm. The underlined layers are n-doped to 2.5×10 11 cm -2 .

[0163] As described in step 1, in combination with the attached Figure 3, the energy band diagram of the quantum cascade laser is obtained, by which the gain spectrum, spontaneous emission spectrum and refractive index change of each local point can be calculated and their analytical expressions can be extracted. It is analyzed that the specific optical working conditions of the device are that the optical wavelength should be between 4.65 μm and 4.75 μm and it is the TM mode (determined by the gain).

[0164] As described in step 2, the optical transverse mode of the device is solved, and the optical transverse mode is drawn in the appendix Figure 2 and the optical refractive index distribution map is marked. In the appendix Figure 4 shows the modal gain spectrum of the device under different injection conditions in the electrical simulation, and in the appendix Figure 5 shows the modal spontaneous emission spectrum of the device under different injection conditions in the electrical simulation.

[0165] As described in step 3, each parameter required for step 4 is extracted as shown in the following table.

[0166] Table 1 Parameter Information

[0167]

[0168]

[0169] As described in step 4: The carrier rate equation and the optical traveling wave model can be jointly solved. For the distributed feedback Bragg grating quantum cascade laser, the optical traveling wave equations are in the following form:

[0170]

[0171] ν g is the group velocity of the light wave. Assuming a pure refractive index coupling grating, the first-order diffraction coefficient κ1 of the first-order Bragg grating is κ -1 = 30. δ is the phase shift factor relative to the Bragg wavelength, and its expression is:

[0172]

[0173] where c is the speed of light and Λ is the grating period value equal to 110.3125 nm.

[0174] The remaining several equations that need to be coupled and solved are as follows:

[0175]

[0176] where the carrier rate equation constitutes the constraint conditions for the optical traveling wave e b (z,t) and e f (z,t). Since almost all the parameters are obtained in steps 1, 2 and 3, only specific boundary conditions are needed to solve.

[0177] The boundary conditions are as follows:

[0178] e f (0) = R l (ω0)e b (0)

[0179] e b (L) = R r (ω0)e f (L)

[0180] In the formula, it is assumed that the laser ranges from 0 to L along the z-axis direction, where L is the length of the laser, which is set to 600 microns here. R l is the reflectivity of the left end face of the device, with a value of 0.9. R r is the reflectivity of the right end face of the device, with a value of 0.1.

[0181] The longitudinal mode results obtained by solving in Step 4 are as Figure 6 shown, and it can be seen that the detuning is 0.02 microns.

[0182] In summary, the advantages of the present invention are mainly reflected in the following two points: First, it has good self-adaptability. Starting from the actual device structure and physical mechanism, the present invention does not require a large number of externally introduced parameters compared with other methods. For a specific structure, it can adaptively calculate its parameters to solve the longitudinal mode. Second, it utilizes advanced TCAD technology compared with ordinary methods and actually constitutes a key component of the laser TCAD simulation. Combined with other components, it can form a complete commercial TCAD laser simulation software.

Claims

1. An adaptive longitudinal mode selection method for laser TCAD-level simulation, characterized in that The optoelectronic characteristics of any given structure are effectively extracted using optoelectronic simulation technology at the laser TCAD level, including the transverse optical mode, the ABC parameters in the carrier rate equation, the gain spectrum, spontaneous emission spectrum, and refractive index change in the form of a parametric analytical formula; subsequently, the optical traveling wave model and the special carrier rate equation extracted from the electrical simulation at the laser TCAD level are jointly solved; specifically, the following steps are included: Step 1: Electrically simulate the laser structure by solving the diffusion-drift Poisson equations at each bias V and extract the analytical form of the gain spectrum; Step 2: Solve the transverse mode of the laser, calculate the modal spontaneous emission spectrum and modal gain spectrum; Step 3: Extract the phenomenological carrier rate equation model from the results obtained by the electrical equation simulation; Step 4: Solve the coupled equations of the carrier rate equation and the traveling wave optical model.

2. The method for adaptively selecting longitudinal modes for laser TCAD-level simulation according to claim 1, wherein The specific content of Step 1 is as follows: Step 1.1: Assume that when the bias is below the lasing threshold, the photon density term U in the photon rate equation m (t) is 0. Since the photon density is small, the stimulated emission term is also small. Delete the stimulated emission term R in the photon rate equation and the electrical equation. st ; Study a slice along the z-axis, and the electrical equation to be solved is two-dimensional or one-dimensional. Step 1.2: After spatially expanding the equation by the finite volume method, solve it by the Newton method, and sequentially solve by scanning the current bias from 0 to a large current bias until the gain spectrum becomes transparent at a certain frequency, that is, neither absorbing nor amplifying, and then continue to scan to 100 times the size of this transparent current bias to avoid missing the gain at a specific frequency; this step requires the use of a gain model, and all the parameters required for the gain model can be obtained by solving the electrical equation; Step 1.3: According to the specific type of active region, different gain models are adopted for the bulk material active region and the multiple quantum well active region: A. The gain model for a single r position in the bulk material active region is: where q is the electron charge, h is the Planck constant, k is the electron wave number, is the square of the free electron mass, ε0 is the permittivity of free space, n τ is the refractive index of the optical frequency, c is the speed of light, v is the photon frequency, E ij = hv; f i and f j is the Fermi-Dirac distribution function of E i and E j energy levels, D r (E ij ) is the joint density of states between two transition energy levels, M is the momentum matrix. Under the condition of momentum conservation, the gain for all possible ij energy levels needs to be integrated and summed in momentum space to obtain the gain spectrum g(r, hv) at a single point r in space; B. For the quantum well active region, the gain model needs to correct the joint density of states D in the above equation r (E ij ) is the two-dimensional density of states: where m r is the reduced effective mass, d z is the quantum well width, is the reduced Planck constant; Step 1.4: Adopt a unified spontaneous emission model for the active region and the passive region. The general form of the spontaneous emission model r sp is as follows: where D opt is the photon density of states, and r sp is also transformed into r in the same way as the gain spectrum sp (r, hν); Step 1.5: Calculate the gain, spontaneous emission, and refractive index change parametric stimulated emission term and the refractive index change term caused by carrier injection from the momentum space through electrical simulation using a strict physical model: In the formula, is the gain spectrum in an analytical expression form, is the photon energy, N is the cold carrier concentration in the active region. It is assumed that the electron-hole concentrations in the active region are equal and stimulated emission only exists in the active region. Other parameters are fitting parameters; is the refractive index change in an analytical expression form. Other parameters are fitting parameters except for the photon energy and the cold carrier concentration in the active region.

3. The method for adaptively selecting longitudinal modes for TCAD-level simulation of a laser according to claim 2, wherein The specific content of Step 2 is as follows: Step 2.1: Neglect the refractive index change caused by gain and carrier injection, use the cold cavity model to solve the optical transverse mode, and set the frequency scanning range to the gain spectrum range obtained in Step 1.5; Step 2.2: Normalize the parameters obtained by electrical solution at each bias V in Step 1.2 to modal parameters according to the transverse mode weighting; where Σ is the regions of each xy plane along the z-axis where the transverse modes of the entire device lie, and Σ ar is the area of the active region, by evaluating the respective modal parameters at a specific point z0, φ 2 (x, y) is the optical transverse mode distribution, Γ is the mode confinement factor, α eff )z) is the modal loss, is the polarizability in the frequency domain, Δn(z, ω0) is the refractive index change, and g(z, ω0) is the gain spectrum.

4. The method for adaptively selecting longitudinal modes for laser TCAD-level simulation according to claim 3, characterized in that The specific content of Step 3 is as follows: Step 3.1: Starting from the rigorous spontaneous emission term used in the actual electrical equation and the rigorous expression under steady state conditions is as follows: where E g is the band gap between the i and j energy bands, and its empirical expression is: In the formula, N and P are the electron and hole concentrations in the active region, N0 and P0 are the electron and hole concentrations in the active region at the equilibrium state, and Br is the spontaneous recombination coefficient; Step 3.2: Extract the Br spontaneous recombination coefficient by the method of weighted average. First, it is emphasized here that is an implicit function of the electron N and hole concentration P in the active region. When solving, the spatial discrete form of the variable is used The part with [m] in the formula refers to the m-th discrete variable of the spatial discrete form of the active region variable, and M is the total number of discrete points in the active region; Step 3.3: Extract the non-radiative recombination terms of the actual steady state for the A and C coefficients of, and divide into Auger recombination terms and SRH recombination terms where A e and A h are the SRH decay rates of electrons and holes respectively, τ e and τ h are the SRH time constants of electrons and holes respectively, C e and C h are the Auger recombination coefficients of electrons and holes respectively; Step 3.4: Obtain the expressions of the A and C coefficients by the method of weighted average; Step 3.5: Introduce a simplified carrier rate equation established only along the z-axis in the electrical equations for the coupled solution of the diffusion-drift Poisson photon rate equation by the Newton method; Wherein, N(z,t) is the cold carrier concentration in the active region without distinguishing electrons and holes, is the slowly varying envelope of the optical electric field strength along the z-axis, J(z,t) is the injection current density, η is the injection efficiency, n eff is the effective refractive index of the transverse mode, Γ is the mode field confinement factor, d is the thickness of the active region, ω0 is the reference frequency, Σ ar is the area of the active region; g als (ω; N(z,t)) is the analytical expression of the gain extracted in Step 1.5.

Images