Probabilistic Markov chain modeling of photonic crystal surface emitting lasers

The photonic crystal surface emitting laser (PCSEL)¹ has emerged as a new class of laser diode. Interest has been sparked by characteristics such as high brightness, and the control of polarization and beam shapes with the PCSEL is the subject of significant current research activity worldwide.² Two-dimensional (2D) distributed feedback results in the PCSEL operating in a single spatial and spectral mode with a very large area, hence power. From simulation and/or the measurement of band structure of the photonic crystal (PC), the coupling constants that represent the coupling of waves in 1D and 2D can be determined.³  Figure 1(a) schematically shows these microscopic scattering parameters for a PC with square symmetry. With these scattering constants, the scattering of photons within the laser can be described. However, device level simulation of the loss mechanisms of PCSELs [shown schematically in Fig. 1(b)] is complicated due to the large, yet finite 3D structure.

Several methods have been developed for PCSEL modeling. The finite difference time domain (FDTD) method is very powerful and utilizes approximate solutions to the associated system of differential equations.^4,5 However, the computational resource requirements can be a significant hurdle for FDTD of very large, yet finite PCSEL devices.⁶ The plane wave expansion (PWE) method describes the field as the superposition of a complete set of plane waves and then uses Floquet–Bloch representation to impose the periodicity of media.⁷ It is computationally less taxing than FDTD but can only consider an infinite PC area, which is a major drawback. Coupled mode theory (CMT) has allowed threshold gain analysis in PCSELs of limited area.⁸ However, it can be very complex in terms of analytical mathematics. Modal index analysis (MIA) views the two-dimensional PC as a multilayer waveguide with periodic longitudinal discontinuities, which makes the MIA method versatile in dealing with device level parameters such as boundary reflectors. However, only square PC atom shapes have so far been considered, limiting applications.⁹ 

In this paper, we describe the application of a probabilistic Markov chain to the simulation of PCSELs. The presented methodology incorporates microscopic input parameters (the in-plane scattering coefficients and those related to radiative and parasitic internal loss)³ and determines macroscopic device level PCSEL parameters such as in-plane optical power loss (allowing threshold gain and slope efficiency to be deduced) and the near-field profile. Following a description of the model, convergence criteria are discussed, and commissioning is completed by demonstrating agreement with state-of-the-art device results in the literature. We then go on to use the model to explore aspects of PCSEL design.¹⁰ The effect of increased in-plane coupling loss and boundary reflectors on the parasitic in-plane loss is described. The effect of spatially varying the reflectivity at the boundary is then explored.

Figure 1(a) schematically represents the microscopic scattering mechanisms for a PC structure with fourfold symmetry (square lattice). Each PC atom is considered as an element of a square. The square symmetry of the PC lattice defines cardinal directions for the transmission of power within the mode. Light traveling in a particular direction (yellow arrow) may remain unscattered or may be subject to 1D scattering resulting in retroreflected light, or in 2D scattering resulting in orthogonal in-plane scattering of optical power. Additionally, out-of-plane scattering (κ_rad) may occur, contributing to output power and being analogous to the mirror loss in a Fabry–Pérot laser. This is the so-called “absorbing state” that terminates a Markov chain.¹¹ Similarly, the internal loss can be defined.

The scattering coefficients (κ_1D and κ_2D) can be calculated from the corresponding coupling constants (κ₁ and κ₃), which were originally used to describe the coupling of waves propagating in the square-lattice DFB structure of PCSELs.³ The scattering coefficients describe how the photon density is distributed within the PC structure. These coefficients typically have units of inverse length (cm⁻¹), but for our model, they are converted to scattering probabilities at each PC atom, with four scattering events possible, being the aforementioned 1D (κ_1D), 2D (κ_2D), out of plane (κ_rad), and internal loss (κ_i) scattering. κ_1D and κ_2D are related to the PC structure. It is possible to engineer κ_1D to almost 0 cm⁻¹.^12,13 κ_2D's value can relate to the asymmetry of PC structure.¹³ The asymmetry of the PC structure also affects κ_rad. κ_i is generally associated with free-carrier loss in the waveguide structure. In this paper, κ_2D is assumed to be equal for scattering into the two possible resultant directions. It is possible to make these two scattering rates inequivalent within the model to accommodate more complex PC structures.¹³ 

The grating area of the PCSEL is described as a matrix of PC atoms with an associated optical power (cardinal directions, internal loss, and emitted power). A time step is defined, which is the time required for light to travel from one PC matrix element to the next. During each time step, new photons are generated within a PC element equally, propagating in each cardinal direction, and those photons in the element in the previous time step transit to a neighboring element and undergo scattering. This procedure is repeated until the system converges.

We note that this model requires a priori knowledge of the scattering constants, and these may be deduced by other simulation,^4–9 experimental methods,³ or assigned arbitrary values to explore new PCSEL device concepts. The use of optical power and no intrinsic solutions of the wave equations result in no consideration of phase difference for power reflected at the PCSEL boundary. While FDTD and our model are based on time-iterative computation, PMC modeling solves neither Maxwell's equations nor any time-dependent difference equations as in FDTD.⁴ We simulate the physical processes of photon power distribution inside the PCSEL device in a time-iterative manner. This requires a-priori determination of these parameters yet reduces computer run-time and memory requirements as compared with FDTD.⁴ 

The model is run with power being continually added to the mode and lost through the edges, through surface emission, and through internal loss. The simulation stops at a time step where convergence is obtained, which is based on the accuracy required for the in-plane loss described in the following. The in-plane power lost at the edge of the PCSEL can be gathered from the edges of the matrix, and the out-of-plane power can be calculated from every atom in the matrix. The matrix is updated after calculating the power transfer of every atom in every direction. In most current work discussing gain measurements of laser didoes, threshold gain, mirror loss, and internal loss are usually discussed in terms of whole reciprocal centimeters (cm⁻¹).^14–22 We, therefore, set our convergence limits to be to within ±0.1 cm⁻¹. However, this can be changed according to requirements.

Figure 2 shows how the calculated α_// varies with elapsed time steps for a range of PCSEL sizes and input parameters as discussed earlier.²³ The sizes are the number of PC atoms along a side of the PCSEL (e.g., 100 × 100 has 100 PC atoms on a side). Also plotted are normalized 2D maps of the surface emitted power from the 700 × 700 PCSEL at a range of time steps. This result should, therefore, depict the near-field of the lasing mode once convergence is achieved. We determine values of dα/dt that provide values of α_// to be within 0.1 cm⁻¹ of those at “infinite time,” which we assume to be at 100 000 time steps. From Fig. 2, we also observe the following: larger size PCSEL devices have lower initial and final values of α_//, and that smaller PCSEL devices reach convergence within fewer time steps. Additionally, by eye, we see very little difference between near-field patterns in insets C and D, due to a variation in surface emitted optical power within each atom changing by no more than 1 part in 10⁴.

Figure 3 shows the calculated in-plane loss as a function of PCSEL size using input data for a triangular PC atom on a square lattice.²³ This report gives κ_rad = 38 cm⁻¹ and κ_i = 5 cm⁻¹, and we extract scattering coefficients³ and κ_1D = 1334 cm⁻¹ and κ_2D = 386 cm⁻¹ from the reported band structure.²³ The inset is the calculated near-field pattern for a 700 × 700 PCSEL. The 3D mode shape is in qualitative agreement with other reports.²⁴ The trend of the curve is in qualitative agreement with MIA simulation results of α_// as a function of PCSEL size.⁹ 

For this PCSEL of size 700 × 700, threshold current and slope efficiency have also been reported in addition to the band structure.²² Our simulation indicates α_// = 17.9 cm⁻¹ predicting a slope efficiency of 0.74 WA⁻¹. This is in excellent quantitative agreement with the experimentally obtained slope efficiency of 0.73 WA⁻¹. Our source of errors in this calculation are in our convergence limits (0.1 cm⁻¹) and in the practicalities of extracting input data from the literature.

The model now allows a range of effects to be explored. We first investigate the effect of modifying the microscopic scattering coefficients on in-plane power loss. In Fig. 4, we consider a PCSEL, where the κ_1D and κ_2D input scattering coefficients are both simultaneously multiplied by a scalar M (0.5, 1, 2, and 3).

We see that multiplying the scattering coefficients does not significantly change the curve shape, and significant in-plane loss is still observed for a 700 × 700 PC PCSEL (α_// = 10 cm⁻¹ for M = 3 as compared to α_//4 = 17 cm⁻¹ for M = 1). This result agrees with the general trend within the canon of increasing PCSEL area to achieve high slope efficiencies and powers.²⁵ 

While our probabilistic Markov chain simulator cannot deal with different reflected phase conditions, we can consider the effect of reflecting light back into the PCSEL with the same phase as the emitted light. We consider a range of power reflectivity, R, from 0% to 95% at the lasing wavelength. Again, the general shape of the curve (shown in Fig. 5) is unchanged, but the results indicate that with modest reflectivity at the boundary of the PCSEL, the in-plane loss can be significantly reduced.

This model also allows the spatial variation in reflectivity around the perimeter of the PCSEL. Figure 6(b) plots the in-plane loss for four cases. Type A is a normal PCSEL device with no perimeter reflector; type B is a 95% reflector evenly spreading around the edges of a PCSEL device; for types C and D, each edge of the device has been evenly divided into three parts: type C with 95% reflector in the center of a side, type D with 95% reflector at each corner. In Fig. 6(a), it can be seen that the reflector can reduce in-plane loss (compare A and B). It can also be seen that a reflector at the center of the perimeter side has a greater impact upon reducing α_// than the corners of the perimeter (compare C and D).

The near-field pattern for all cases is also plotted in Fig. 6(c). The scaling of four nearfield modes in Fig. 6(c) is the same. It can be seen that a flatter, more uniform near-field is achieved using reflectors (compare A and B). We also note that the mode pattern is significantly modified by the use of more complex feedback schemes as shown for C and D.

In this paper, probabilistic Markov chain modeling of PCSELs has been discussed. The ability to input microscopic PC scattering parameters and describe device-level PCSEL performance has been demonstrated. The effect of scaling the in-plane scattering parameters, using reflective boundaries on in-plane optical power loss, and the use of spatially varying boundaries on the near-field have been explored.