Dynamic Modeling of Stress-Induced Defect Expansion in VCSELs

Many failures of semiconductor-based oxide confined vertical cavity surface emitting lasers (VCSELs) are closely related to the generation and expansion of defects in the device structure. However, existing research has predominantly focused on the static study of defect morphology, with little attention given to analyzing the dynamic process of defect expansion, which limited our ability to predict device random failures due to lack of understandings on defect generation and expansion. To address this issue, we present a macroscopic phenomenological evolution model that describes the dynamic expansion of defects in VCSELs, in which the expansion of defects is treated as an anisotropic lattice strain diffusion process. We further exploit a diffusion-limited aggregation (DLA) method in solving the diffusion equation, which describes the random propagation and aggregation of strain in the vicinity of highly strained areas, resulting in defect formation when the stress from accumulated strain surpasses the bonding force of the atoms in the lattice. Our simulation result manages to replicate the dendritic expansion morphology of defects, aligning with experimental observations very well. Our model also predicts an accelerated defect expansion process, which is again consistent with the experimental result. This model finds relevance in applications such as random failure prediction through device aging, burn-in condition setting in device screening, and device structural and/or material design improvement for mitigating defect expansion.

revealed that dark line defects are the primary cause of failure in GaAs-based oxide confined VCSELs [4], [5], [6], [7].These defects are generated due to the excessive mechanical stress on the oxide layer formed by the oxidation of the AlAs or Al x Ga 1-x As (x≥0.9) layer in oxide confined VCSELs.The conversion of Al x Ga 1-x As (1≥x≥0.9) to alumina (Al x O y ) during the oxidation process results in a 6%∼20% volume shrinkage [8], [9], which gives the source of defects [10], [11].Additionally, the mesa columnar structure etched out before the selective wet oxidation process can cause subtle mechanical damage near the mesa edge, as well as a large number of suspended crystal bonds in the semiconductor material, which further contributes to potential sources of defects.Such defects arising from the edge of the mesa can gradually expand into the active region and become one of the most prominent failure modes occurring in oxide confined VCSELs [11], [12], [13].Furthermore, the high Al content of many distributed Bragg reflector (DBR) layers can also result in partial oxidation during the oxidation process, leading to accumulated stress at the mesa edge and contributing to potential defect generation.Although the directly injected current density at the mesa edge is low, the nonequilibrium carriers may not be negligible as they can be generated by the light absorption inside the wells around the mesa edge.These carriers play a crucial role in activating defects around the mesa edge and triggering the strain diffusion inwards.The elevated temperature enhances the lattice vibration, which contributes to the defect activation at the mesa edge and speeds up the strain diffusion and aggregation process.Driven by external forces such as current and/or temperature, these edge defects can again diffuse into the active region.Once they reach the active region, defects would rapidly climb and expand, due to the built-in strain inside the quantum wells, which eventually leads to device failure once the defect density reaches a critical level.Although there are advanced technologies to mitigate the possible defects around the mesa sidewall, particularly near the oxidized ring [12], [14], the traditional processing technology [15], [16] seems still common for fabrication simplicity and/or cost-effectiveness.Moreover, although in advanced VCSELs, the remaining reliability issue is the wear-off degradation which shows a slow and gradual deterioration of the device performance, we still found that the random failure rate in some production batches of devices was not to our satisfactory, even when the mesa and oxidation were processed with the state-of-the-art technology.And such random failure caused by the latent damage or hidden flaw cannot always be screened out through the burn-in process.Our randomly failed VCSELs typically showed the dendritic type of defect patterns, which agreed with the reported result in literature [12], [17], [18], we therefore wanted to find how such featured patterns could possibly be formed.Understanding the dynamic process of defect generation and expansion in VCSELs is essential for improving their reliability.
Among the many research works on VCSEL's reliability and failure analysis, some relevant studies have been reported by reference [4], [5], [6], [7], [10], [11], [12], [13].Although analyses on the root cause and intrinsic mechanism of defect formation and expansion were mainly based on theoretical speculations [4], these studies managed to provide several possible degradation modes of VCSELs and explained the nature of defects, their origin, and potential mechanisms of their expansion.However, only postmortem studies have been carried out on degraded devices for reliability analysis in existing works, in which morphologies were only observed after the completion of defect expansion, making it impossible to obtain a complete picture of the entire defect growth process.Due to the multiple-origin nature of defect generation and complex interactions in defect expansion, the dynamic process that forms the observed random pattern of the morphology and the characteristics indicated by these patterns were not fully understood.A variety of defect simulation models exist, including first-principle physics-based atomic models [19], [20], molecular dynamic models [21], [22], and discrete dislocation models [23], [24].However, the outcome results from these models have limited time and space scales, not to mention that the simulation requires expensive computation resources, making the existing methods unsuitable for macroscopic-scale simulation of VCSEL structures with heterogeneous materials and complex geometry.
To address the aforementioned issues, a dynamic model together with the associated efficient simulation method is needed for describing the generation and expansion of defects at the device structural scale.
In this work, we propose a phenomenological model that describes the defect dynamics in VCSELs.An anisotropic defect diffusion equation with the device operating condition incorporated is first drawn to govern the motion of defects.A diffusion-limited aggregation (DLA) method is then introduced to simulate the random dendritic morphology of defects observed in VCSELs.The defect's dynamic diffusion process in VCSELs can be effectively simulated with its outcomes closely resembling the observed random patterns.Several featured parameters extracted from these patterns also match the observed defect morphology very well.The predicted acceleration, rather than deceleration, of the defect expansion also agrees with the experimental observations.The rest of the paper is organized as follows.Section II describes the derivation of the phenomenological model that governs the defect motion, in which an anisotropic defect diffusion equation is obtained.The solution technique through the DLA approach is then established in Section III, which substantially differs from the traditional method of solving the diffusion equation.Section IV is dedicated to showing the simulation results and their comparisons with experimental observations, with discussions and further quantitative analyses.Finally, this work is summarized by a conclusion.

A. Anisotropic Defect Random Walking
Firstly, we assume that a defect moves in the lattice following a random walking process described by: where R N+1 and R N are vectors indicating the position of the defect at step N+1 and N, respectively; U N is a unitary matrix that defines the walking direction.In a Cartesian coordinate system, once the walking direction is specified by a pair of polar angles (θ, ϕ), the unitary matrix will take the form of: Next, we will need to find the likely directions that the defect prefers to move in the random walking process.Assuming r the walking step, R the position of the defect at step N+1, p(r) the probability of walking from position R − r to R (which is independent of R), and P N +1 (R) the probability of the defect at R at step N+1, we can establish the following recursive relationship: Expanding the right-hand side of (3) into the Taylor series, we obtain: If it takes a period of Δt for the defect to move from step N to N+1, we find [25]: where ρ = P N (R), ∂ρ/∂t = [P N +1 (R) − P N (R)]/Δt, and the diffusion coefficient D = r • r /6Δt.It describes how the density of defects evolves over time, with the Laplacian term representing the spatial spread of defects and the diffusion coefficient D characterizing the rate of diffusion.The formula is based on the assumption of isotropic diffusion, meaning that the motion of the defect follows a standard isotropic diffusion process.However, experimental observations showed that the defect expansion in InGaAs-AlGaAs-GaAs systems is not isotropic.Rather, it tends to grow in favor of the [110] direction [4], [26].Fig. 1 shows the defect pattern we observed by plan-view transmission electron microscopy (PV-TEM) in failure VCSELs.In all 6 samples, we found that indeed [110] is the dominant direction of defect growth.Therefore, the isotropic diffusion model described by (3) must be modified to capture this feature.Although qualitatively we know that the interplanar crystal Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.spacing, as well as the valence bond density and strength, play an important role in determining the defect growth direction, we show that the anisotropy of defects' motion in zinc blende crystalline structure possessed by InGaAs-AlGaAs-GaAs systems can be quantitively described with a strain-stress model which matches with experimental observations.
The absolute stress tensor σ endured by the lattice in a crystal is related to the lattice strain tensor ε by Hooke's law [27]: where C is the elastic modulus of the material.The lattice strain tensor ε is produced by the lattice relative displacement vector u: Therefore, the net (i.e., the relative) stress vector F of the lattice is: Phenomenologically, we assume that, in any space direction i, (a) the defect growth probability ρ i is proportional to the lattice displacement vector u i , and (b) the defect growth rate ρ i /t is proportional to the lattice stress component F i , which turns (8) into: with δ introduced as a coefficient to match the dimension on the two sides of ( 9).These assumptions are reasonable because the defect is a direct result of lattice displacement and the stress is the driving force of defect generation, as reported by previous experimental studies [4], [5], [6], [7].It is worth noting that unlike the scalar quantity ρ in ( 5) that defines the probability of the defect at a specific location at certain step, the vectorial quantity ρ = [ρ i ] 1×3 defines the probability of the defect's random walking direction from one step to the next.
In InGaAs-AlGaAs-GaAs systems with the zinc blende crystalline structure, the elastic modulus bears the cubic symmetry.As a result, there are only three independent components of the elastic modulus: C 11 , C 12 , and C 44 .Eq. ( 9) can therefore be simplified to: where x, y, and z, indicate the [100], [010], and [001] directions, respectively.
Apparently, (10), shown at the bottom of the next page, takes the form of an anisotropic diffusion equation, with a different probability of defect random walking in each different direction reflected by the different diffusion coefficients in the corresponding direction.

B. Reduced Two-Dimensional Equation
Reliability studies of VCSELs have found that defects are mainly confined inside the active region, which may likely be due to the built-in strain in the quantum well and/or barrier and the heat generated by the high carrier density as they collide with the lattice and recombine non-radiatively [28], [29].Our experimental observations also support this finding.In Fig. 2(a), the top plane-view (PV) transmission electron microscope (TEM) observation of a failure VCSEL is presented, in which the active region underneath the oxide aperture is not removed.A spreading defect pattern is clearly visible.However, the defect pattern disappears after the active region is removed by etching, as depicted in Fig. 2(b).Fig. 2(c) further shows the cross-section (XS) TEM image of the same device sample, in which the defects are evidently concentrated inside the active region.
Fig. 3 depicts the spatial morphology of the defects typically observed.It reveals that the defects originated from the curved inner edge of the oxide layer and spread downwards into the active region.As the active region has the built-in strain and is heated relatively heavily, additional stress coming from outside can trigger the generation of a large number of local defects due to lattice debonding and dislocation.Since the generation of the defects releases the strain, their spreading is confined inside the strained layers and their close neighborhood.They cannot climb far vertically inside the unstrained region for their momentum is lost as the strain disappears.The difference between the defect spreading scales inside the active region (001) plane and along the vertical [001] direction is approximately 2∼3 orders of magnitude, with the featured expansion dimension of defects in the (001) plane significantly larger than that along the vertical direction.This fact suggests that a two-dimensional (2D) diffusion model is sufficient in describing the defect evolution.Hence, we further reduce our (8) from three-dimensional (3D) to 2D.
The 3D defect diffusion (10) can therefore be reduced to its 2D form shown below: Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Due to the symmetry of the cubic crystalline structure, the defect diffusion equation in the (001) plane can further be simplified as: with ρ(x, y) defined as the probability of the defect's random walking direction at position (x, y) in the (001) plane.Although ρ returns to a scalar quantity, its diffusion rate along the different directions in the (001) plane is different, as evidenced by the different diffusion coefficients in the x and y directions.
We set ρ = 1 for t = 0 at the origin (x = 0, y = 0) and solve equation (11b) for t > 0, the solution of ρ(x, y) gives a distribution shown in Fig. 4.Those ellipses represent the solution ρ(x, y) at different time steps.As time elapses, the ellipse expands outwards.Since the diffusion pattern is a set of ellipses rather than cycles, it is evident that the diffusion is not isotropic, but with a dominant direction along the long axis of the ellipse, which is the [110] direction, consistent with the experimental observations.
Since ρ is the probability of the direction that the defect takes in its random walking process, after a deterministic distribution of ρ(x, y) is obtained by solving (11b) in the neighborhood of (x, y)ࢠC, with C denoting a circular area in radius R, a random number generate (RNG) will be exploited to select an azimuthal angle ϕ by forcing the RNG to follow the distribution of: Equation ( 12) enforces a constraint to the distribution of the random number generator which gives the probability of the walking direction.For example, if ρ(x, y) = ρ(r cos ϕ, r sin ϕ) gives a non-uniform distribution along angle ϕ, D(ϕ) will also be non-uniform and will be in proportion to ρ.The random number generator will pick up the direction along ϕ as frequently as D(ϕ) and therefore as frequently as ρ.As such, the defect propagation will have the highest chance to take the direction along which ρ is maximum, but also with a probability to deviate from that direction.

C. Driving Source
The homogeneous (11b) must be modified to include a driving term (the defect source) to obtain nontrivial solutions.Moreover, the ambient temperature and the bias current have a strong impact on defect evolution in semiconductor lasers [30], [31].We Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE I SIMULATION PARAMETERS AND VALUES
thereby incorporate these effects into our equation by an aging factor A(I, T ), shown in (13), which is introduced by following the conventional Arrhenius and polynomial models to describe the temperature and the current dependence, respectively [18], [31], [32], [33].(13) where I -the bias current T -the ambient temperature I ref -the reference current T ref -the reference temperature k -the Boltzmann's constant E a -thermal activation energy n-current power index Our final equation can be expressed as: where s(x, y) represents the distribution of the source defect generation rate.Depending on different device structures, the initial location of the source defect can be placed at any given place, which makes ( 14) general enough for studying the defect evolution in a variety of VCSELs with different geometry and/or in different sizes.To avoid over-complication, we have also assumed that the temperature and current aging effect on the diffusion coefficient and the source are the same.As a result of our fitting to experimental data, we found that to use the same aging factor was indeed sufficient, as a single set of extracted activation energy and current aging index used for both aging factors could achieve a good agreement.The activation energy and the power index need to be extracted through experiments [18], [32], [33].In this work, these values were extracted through separate experiments and listed in Table I in Section IV-B, which also matches very well with those reported in the literature [18], [32], [33].

III. SOLUTION TO DEFECT EXPANSION MODEL
There are two different approaches in modeling defect expansion through the set of proposed governing (1), ( 2), (12), and (14).Distinguished by different methods to set the distribution of the defect generation rate s(x, y) in ( 14) and to identify the defect trajectories, these approaches are commonly known as: (a) the diffusion-limited growth (DLG) method, and (b) the diffusion-limited aggregation (DLA) method.
In the DLG method [34], the defect pattern is directly formed by the trajectory of its diffusion process.The defect expansion process follows a sequence which the diffusion starts from a defect's initial position (i.e., the position of the seed defect) and randomly walks outward.The direction of the random walking at each step in (1) is determined by solving (14) with its source given at the position of the previous step.Its solution is subsequently substituted into (12).The RNG then selects an angle based on the distribution given by the result of (12), as described in the last paragraph of Section II-B.Through (2) (in its 2D version by setting θ = 0), this angle finally determines the position of the defect at the followed step.The above process is repeated until a maximum number of random walking steps is reached.The trajectory of the defects is given by the defect generation rate [i.e., s(x, y)] throughout the entire process.The pattern of the grown defects is single crack-like with no brunches, similar to the path of classical random walking, as shown in Fig. 6(a).We did see such defect patterns in some other device structures, but not in VCSELs.At least they are not the majority defect patterns in VCSELs from top-view observation in the (001) plane.
In the DLA method [35], [36], [37], however, the defect pattern is not directly formed by its trajectory of diffusion path, rather, it is generated by defect aggregation through the diffusion process.The defect pattern begins in a specific region some distance away from the defect's starting point, it describes the gathering of the defects after random walking.Initially, an aggregation point is randomly formed at a location separate from the defect source.(In this work, the defect source is around the oxide aperture, whereas the aggregation area in which the defect pattern appears is the active region made of strained quantum wells.)The defects then randomly walk away from their origin by following the aforementioned DLG process.The defect's walking stops if it either touches the aggregation points (i.e., the defect goes to a position adjacent to the aggregation points) or moves out of the computation window.As such, the aggregation point expands and forms the defect pattern.As evidenced by the images presented in Fig. 6(b), the defect patterns obtained by the DLA method clearly show the featured dendritic structures, which agrees with what has been observed in VCSELs.The DLA process can be illustrated by Fig. 5.A step-by-step procedure of the DLA algorithm adopted in this work is given in Appendix A.
As being referred to in explaining the DLG and DLA methods, Fig. 6 shows the numerically simulated defect patterns obtained by the two different approaches, respectively.These simulation examples vividly illustrate the distinctions between the two methods.Notably, only the patterns generated by the DLA method exhibit the dendritic structures, which is considered to be a crucial feature that needs to be captured in modeling VCSEL's defect evolution.Consequently, in the subsequent analyses, we will exclusively employ the DLA method, rather than the traditional DLG method, for addressing the defect evolution model presented in Section II.
Fig. 7 provides a qualitative comparison between the simulation results obtained using the 2D DLA method and the experimental data.In numerical simulations, we initiated the defect origin at different locations by following the observation of the measured defect patterns in different TEM images.In Figs.7(a), (b), and (c), the origin points of defects are the mesa edge at the lower right corner, mesa edge at the lower right

A. Sample Preparation and Defect Pattern Observation
The experimental samples are commercial oxide confined VCSELs with an emission wavelength at 850 nm.The crosssectional schematic structure of the sample device is illustrated in Fig. 8.The epitaxial structure was grown on the N-doped GaAs substrate on the (001) crystal plane by MOCVD.The active region consisted of 3 pairs of GaAs/AlGaAs strained-layer quantum well and barrier.The DBR sections were made of a number of pairs of alternating Al 0.15 Ga 0.85 As and Al 0.9 Ga 0.1 As layers.The cylindrical mesa with a diameter of 30 μm was defined by photolithography plus dry etching.The oxide aperture in a diameter of 8 μm was finally formed by steam-oxidation of the pre-inserted AlAs layers in the P-side DBR section.
The samples were packaged in unsealed TO-46 headers and were lit under high temperature and current conditions in Chroma 58604 aging equipment.The typical threshold current and current density values for the tested samples are 0.49 mA and 0.98 kA/cm 2 .Prior to the aging process, the samples went through a burn-in screening under an ambient temperature of 120 °C and a constant bias current of 10 mA (∼20 kA/cm 2 ) for 24 hours.Devices showing early degradation were removed.The remaining samples were then exposed to an accelerated aging test at 85 °C and 12 mA (∼24 kA/cm 2 ).We conducted defect characterization on device samples that stopped emitting light.The samples were cut by dual focused-ion beams (FIB) and then PV-or XS-observed by the TEM, through FFI Helios DualBeam and Talos F200X, respectively.
We conducted a detailed simulation analysis and discussion on the failure defects of three typical samples (referred to as samples A, B, and C) in the top row of Fig. 1.Defects originated at the mesa edge and expand towards the device's interior region, driven by the accelerated temperature condition and current bias.The expansion pattern shows clear random branching, and the expansion trend is primarily along the [110] direction.

B. Simulation Results and Discussions
The simulation parameters and their values for the proposed model are outlined in Table I.
Fig. 9 shows the simulated defect growing patterns of the three failure samples in the top row of Fig. 1.The initial defect distribution is marked in red and placed in the first column of the picture of each sample simulation result.Each horizontal row of images from left to right shows the dynamic evolution of the defects from the beginning to the end for each sample.Although the simulated pattern cannot precisely replicate the exact morphology due to the stochastic nature of defect expansion, it closely resembles the real observed image, indicating that the model is promising for describing the defect expansion phenomena in the experiment.To validate this model, we further studied two hidden features in these images as discussed below.
1) Agreement on Branch Spacing: One of the hidden features of the dendritic pattern is the average spacing between branches.Measured in unit of their width, the first order branches (i.e., the branches directly attached to the main trunk) are gapped by 6∼10 units, the second order branches (i.e., the branches attached to the first order branches) are gapped by 2∼4 units, and the third order branches (i.e., the branches attached to the second order branches) are gapped by 0.8∼1.2units, respectively.In each branch order, by "width" we meant the width of the branch in that order.The absolute measure on the average width of the defect branches is around 87 nm.We traced 40 sets of such branches and extracted the gap distance between the first, second, and third order branching from the measured and simulated patterns shown in Figs. 1 and 9, respectively.Fig. 10 shows the comparison result on these spacings extracted from the patterns simulated and from the images experimentally observed.The branch spacings shown in the simulation agree very well with those in experimental observation, as evidenced by the calculated correlation coefficients, 99%, 96%, and 97% for the first, second, and third branch spacing, respectively.This result clearly shows that, on the branch spacing feature, the simulated pattern is almost 100% similar to the measured one.The simulation also shows that the branch spacings remain stable and are not significantly affected by the source density, location, and distribution, or the initial aggregation point.Limited by our phenomenological mode, we are unable to give a physical explanation of this finding.It only reveals that this feature might be determined by the material's microscopic crystalline structure, for it only has a non-appreciable dependence on macroscopic external conditions.
2) Defect Growth Velocity: Our simulation through the DLA model shows that the growth of defects is accelerating rather than decelerating, which is substantially different from a classical DLG based random walking process with its average growth rate given by the famous expression [39]: where t denotes the time and b a constant.In the DLA process, however, the growth rate can be found as: with the derivation given in Appendix B, where D represents the average distance from the defect source to the initial aggregation point, l the length that the defect grows from 0 to t.In Fig. 11(a), we plotted the degradation of the laser output power as a function of time for the three samples A, B, and C shown in top row of Fig. 1.Although the degradation started at different times, the trend is all similar -in an accelerating mode.Since it is impossible to directly observe the defect pattern in its growing process, we have to use the power degradation rate as the indicator of the defect growth rate, for it is reasonable to assume a linear dependence of the power degradation rate on the defect growth rate.
We calculated the defect growth rate given by ( 16), and have it compared with the power degradation rate.Fig. 11(b) shows the comparison of the normalized degradation rate as a function of the normalized time between the experimental data and the calculated result.The experimental degradation rates of the three samples are extracted by taking numerical derivatives of the degradation curves as shown in Fig. 11(a), whereas the degradation rate predicted by the DLA theory is directly calculated by solving (16).These degradation rates, as well as the time variable, are all normalized in order to show the quality of the matching without knowing the actual ratio between the power degradation percentage and the defect growing length and without calibrating the actual time duration represented by each time step used in calculation.As can be seen in Fig. 11(b), it is clear that the degradation rate predicted by the DLA theory indeed matches with the experimental result after normalization, which indicates the self-consistence of our model.

V. CONCLUSION
With the DLA solution technique, our proposed phenomenological defect evolution model can describe the dynamic expansion process of the defects typically found in semiconductor-based oxide confined VCSEL failure samples.This model successfully captures the main features of the defect pattern as experimentally observed on: (a) the dominant growing direction (i.e., along [110]), (b) the dendritic morphology, (c) the branch spacing, and (d) the accelerating growth rate.The simulation result quantitatively matches with the experimental one on the branch spacings in every order.Although this model and the associated solution technique (DLA) do not reveal much physical insights, it offers a rather accurate tool for modeling the dynamic process of defect expansion.It is therefore appealing to be used for random failure prediction, burn-in condition setting for early degradation screening, and aging condition setting for gradual degradation statistics.Moreover, the comparison between the DLG and DLA simulation reveals that the dendritic defect pattern is a result of inward strain diffusion plus aggregation in a vulnerable region, rather than direct outward growth, which gives us a hint to mitigate the risk of such defect expansion by introducing an isolation ring to separate the vulnerable region (i.e., the active layers) and the area with high density of strain and/or point defects (i.e., the oxide aperture and the mesa sidewall).

APPENDIX A DLA-BASED SOLUTION TECHNIQUE
A step-by-step description of the DLA-based solution method is given below: 1) Define the region of the defect source distribution (e.g., around the oxide aperture) and specify the initial aggregation point (i.e., the defect growth starting point).2) Choose a defect position in its source distribution region to start with.The defect randomly walks away from its source distribution region following (1).This process is computed by following sub-steps: a) Solve (14) for ρ(x, y) with its source positioned at R N ; b) Substitute ρ(x, y) into (12) to find D(ϕ); c) Use an RNG to select ϕ according to D(ϕ); d) Substitute ϕ into (2) (where θ is set to zero for 2D) to obtain U N ; e) Insert U N into (1) to find the next position R N +1 at which the defect arrives; f) Repeat (a)-(f).Steps (a)-(f) form the classical DLG solution approach.
3) If the defect moves to a position adjacent to an aggregation point, or moves out of the computation window, the process (DLG) in the above step 2 stops.Once a defect touches an existing aggregation point, its location becomes an aggregated point.Hence the defect grows as the aggregated point set expands.Record the aggregated points and proceed to the above step 2 again by choosing a new defect position in its source distribution region to start with.4) Repeat the above steps 2 and 3 until the defects started from all different positions in their source distribution region are exhausted.

APPENDIX B DERIVATION OF DEFECT GROWTH RATE
Following the classical diffusion (DLG) process, the average diffusion length L is proportional to the square root of the diffusion time t: where b denotes a constant.
In the DLA process, by assuming that the defect growth length in a time period from 0 to t is l, we find: where D represents the average distance from the defect source to the initial aggregation point.This is because for any point 0 ≤ x ≤ l, the time required for a diffusion process from D to x is given by the integrand in (18) according to (17).The total required time for x to move from 0 to l is therefore given by (18).The integration in (18) gives: By taking the derivative of l with respect to t, we find: )

Fig. 1 .
Fig. 1.PV-TEM images of the failure VCSEL samples reveal that [110] is the dominant direction of defect growth.

Fig. 2 .
Fig. 2. TEM images of a failure VCSEL sample: (a) PV-TEM observation with the active region underneath the oxide layer, (b) PV-TEM observation after removal of the active region, and (c) XS-TEM observation of the defect region.

Fig. 3 .
Fig. 3. Typical defect morphology in VCSELs: (a) plan view (from the top) and (b) cross-sectional view (from the side).

Fig. 4 .
Fig. 4. Solution of the 2D diffusion (11b) shows that the dominant direction of defect random walking is along [110].

Fig. 6 .
Fig. 6.Simulation result of the two different approaches: (a) the DLG method, and (b) the DLA method.The red marked point is the origin point, and the black patterns are the aggregation points.

Fig. 7 .
Fig. 7. Comparison between experimental and simulation results on VCESL's defect pattern.The different patterns are generated by different distributions of the defect source.

Fig. 10 .
Fig. 10.Comparison between the simulated and experimentally observed branch spacings.(r indicates the Pearson correlation coefficient, defined as the covariance of two variables divided by the product of their standard deviations.).

Fig. 11 .
Fig. 11.(a) Degradation of the laser output power as a function of the aging time for three different samples A, B, and C. The optical output power at different aging time is measured at a fixed bias current (6 mA).(b) Normalized degradation rate v. s. normalized time, which shows that the degradation rate predicted by the DLA theory as given in (16) fits the experimental result well.