Numerical model for small-signal modulation response in vertical-cavity surface-emitting lasers

We present a numerical model allowing for simulations of small-signal modulation (SSM) response of vertical-cavity surface-emitting lasers (VCSELs). The model of SSM response utilizes only the data provided by a static model of continuous-wave operation for a given bias voltage. Thus the fitting of dynamic measurement parameters is not needed nor used. The validity of this model has been verified by comparing experimental SSM characteristics of a VCSEL with the results of simulations. A good agreement between experiment and simulations has been observed. Based on the results obtained in the simulations of the existing laser, the impact of the number of quantum wells in the active region on the modulation properties has been calculated and analyzed.


Introduction
Optical data transfer is a growing field of application for semiconductor lasers. In the existing and possible future short-range systems, vertical-cavity surface-emitting lasers (VCSELs) are natural candidates for the emitters in such systems. State-of-the-art optical links based on VCSELs reach bandwidths over 30 Gbps [1][2][3][4][5][6], approaching even 100 Gbps [7] (the four-level pulse amplitude modulation (PAM4) scheme was used in the last reference). A good review of this subject can be found in [6] and [8]. Further improvement in the optical communication link capacity is expected because there is a constant need for increasing the amount of data sent between computers and nodes in data centers and other systems relying on optical data transfer [9]. Future optical links will need emitters that can be modulated with higher frequencies. Further development of such emitters could be intensified if useful models were available. Despite the development in the technology, the existing models for the modulation performance of VCSELs and semiconductor lasers in general are still ineffective. The existing models of modulation response are based on many parameters whose values can be chosen arbitrarily to a certain extent, or have to be fitted to experiment [10][11][12][13][14][15][16]. These models have shown that the approach based on rate equations can give good agreement with the small-signal modulation (SSM) response experiment. This approach allows for determination where the limitations of the obtained modulation bandwidth come from, but not always how to overcome the problem. The greatest challenge is to calculate, not fit, the parameters in these equations in such a way that the results of SSM experiments can be predicted for different operating conditions.
In this paper we present how SSM response of a VCSEL can be calculated using numerical simulations. In order to find the parameters describing the laser's response to the modulation, a reliable model of above-threshold continuous-wave (CW) operation of the VCSEL is needed. In this paper we show that based on such a model developed by the Photonics Group at Lodz University of Technology (LUT), we can develop a model for predicting the SSM characteristics. Although the Photonics Group has been modeling various kinds VCSELs for many years, their work has been focused on static parameters, especially at the laser threshold, although above-threshold characteristic have also been calculated [17][18][19]. Recently, a model of electric dynamic parameters (SSM reflection) was presented [20]. The model presented for the first time here allows for calculating SSM response which is more general much more complicated problem. This model is verified by comparing the modelled and measured characteristics of a GaAs-based high-speed VCSEL.

Model
The entire numerical model of VCSEL operation consists of two parts: (a) a numerical model of CW above-threshold operation; and (b) an analytic model of SSM response under the conditions determined by the CW operation model.
The CW (static) model we used is an already existing model developed by the Photonics Group at LUT. Such models are also developed by other groups, and used to model VCSELs [21,22]. The dynamic model uses the standard rate equations. The novelty of our approach is the fact that we determine the parameters for the rate equations based on the CW simulations only. As a result, we do not need any frequency-dependent measurements to simulate the laser's modulation response. The CW characteristics of the laser depend on many parameters that are not universal material parameters, such as the resistivity of the electric contacts; or can, in real devices, be slightly different than designed (e.g. layers' doping, QW thickness and composition, etc). Because of this, usually, some adjustments of the CW model parameters are necessary in order to get a good quantitative agreement with experiment. In the dynamic model, however, we use only such parameters whose values can be derived from the CW model, without relying on modulation experiments. In what follows, we describe each of these parts in detail. Section 2.1 describes our CW model which is a modification of an already published model [23]. Section 2.2 describes the text-book rate-equations model, improved by us by taking into account the impact of the non-zero photon lifetime that is usually neglected in the literature. Section 2.3 contains original research on how to connect both previous models. We believe that our approach can be used also with different CW models. To the best of our knowledge this the first demonstration of such an enhancement of a VCSEL model.

CW operation model
The numerical model of CW above-threshold operation is based on models of the following physical phenomena: The model presented here is a modification of the model presented in reference [23]. The main modification we made was intended to limit the number of times the optical model is invoked, because in our case, the optical model, which is a 3D fully-vectorial model [24] is the most time-and memoryconsuming piece.
A scheme of the calculations performed by the CW model is shown in the flowchart in figure 1.
First, distributions of the electric potential and the temperature in the structure are determined, for a given voltage, by the thermal-electric model [25,26]. In the next step, optical modes in the cavity are found. The temperature distribution from the previous step is used to take into account its impact on the refractive index in the cavity. The optical gain in the laser's active region is set to 0, so as to calculate the total mode losses. The mode is a solution of Maxwell's equations where the mode is a monochromatic (i.e. having a well-defined frequency) electromagnetic wave. The mode's frequency is a complex number: In the case of a lossy cavity, the imaginary part ω i > 0. This means, that the mode's energy in the cavity E M will decay according to the following formula: so the number ω i describes all the mode losses. In our case, where the gain in the active-region is 0, it includes optical absorption and emission of photons from the cavity. Since Maxwell's equations are linear, the amplitude of the mode (and hence its energy or the number of photons) can be arbitrary. This means that the mode's electric or magnetic field distribution can be multiplied by an arbitrary number. In order to find the amplitude of the mode, it is necessary to have another equation. In the model presented here, this equation is given by the condition that the number of photons generated per unit of time in the active region (G) is equal to the number of photons lost, either due to absorption or emission: The number of photons generated in the active region is given by the following equation: where Act denotes the active region, r = (x, y, z), c is the speed of light, ε 0 is the dielectric constant, n r is the refractive index (it also can depend on r), g is the material gain, ℏω r is the photon energy, and E is the mode's electric field in the standard complex-number notation. The impact of the stimulated emission on the carrier (and hence gain) distribution, manifesting as the spatial hole burning, is taken into account by the presence of the stimulation-emission losses in the carrier diffusion equation: where u is the carrier density in the active region, D is the ambipolar diffusion coefficient, A, B, C are, in general temperaturedependent, mono-molecular, bi-molecular and Auger recombination coefficients [27]; j ⊥ is the perpendicular component of the current density injected into the active region, d is the total thickness of the active-region QWs, and e is the elementary charge. In our implementation, this equation is solved in 2D, and the mode's distribution is averaged in the vertical direction in the QWs. The optical gain distribution is calculated based on the carrier, temperature and mode distributions. First, optical gaing is calculated assuming thermal (Fermi-Dirac) electronic states occupancy. Our simulations showed that such a gain model is not sufficient to describe the observed rollingover of the VCSEL's light-current (LI) curves. To our model we added the effect of spectral hole burning (SHB) using the following relation, similar to one used in the literature [10,16]: The factor ε0n 2 r 2ℏωr ∥E∥ 2 is the density of photons that interact with the carriers. Because only the electric field of the electromagnetic wave contributes to stimulated emission, this factor is simply the energy density of the electric field divided by the photon energy. The factor α has the unit of m 3 and describes the strength of SHB in this system. In our calculations, we fitted α to measured CW LI curves.
The total energy of the mode is given by the following formula: where Res denotes the entire resonator, η 0 is the free space impedance and H is the mode's magnetic field. Although locally the energy densities of the electric and the magnetic fields are usually different (in a standing wave the nodes of the magnetic field coincide with the anti-nodes of the electric field and vice versa), the total energy of both fields in the whole cavity is equal. Then the mode energy can be also expressed in the following ways: The latter form, using the magnetic field, is generally more convenient, as the integrated function is continuous (in the case of non-magnetic materials), contrary to the energy density of the electric field. Because the optical model (Maxwell's equations) makes no distinctions between solutions that differ in the amplitude only, we can use in the calculations electromagnetic field distributions normalized in such a way that the total number of photons in the cavity is equal to 1 and the actual number of photons as a separate parameter. Let us denote the corresponding electric and magnetic field scalar distributions (i.e. the distributions of the squared norm of each of the fields) as M E and M H . Then: Any possible distribution of the mode's electromagnetic field intensities can now be expressed using an amplitude P, being simply the number of photons in the cavity: where F 2 stands for ∥F(r)∥ 2 for F being either E or H. With this notation, we can write equation (3), using equation (2), in the following form: Although the amplitude P cancels out on both sides, the lefthand side does depend on P, because the carrier density and hence material gain g do. The above equation can be written as:ˆA ct cε 0 n r g(r, P) There is a fundamental question, whether this equation has solutions. The optical gain g is a decreasing function of P when g > 0 (see equations (5), (6)). We are interested only in the cases where the left-hand side integral is positive (because the right-hand side is positive), and we can assume, that the left-hand side is maximal for P = 0. If Act cε 0 n r g(r, 0) 2ℏω r M E (r) dxdydz < 2ω i (13) then the mode is below its threshold. In the opposite case, the mode is above the threshold (or at the threshold if there is an equality), and there is a single value of the amplitude that satisfies equation (12). This solution is unique, because the left-hand side is decreasing while the right-hand side is constant. Knowing the amplitude we can calculate the laser output power as the flux of the electromagnetic field through the emission surface of the laser. Using the equations presented in this section, we can write an equation for the average threshold gain: Or, using the light speed in the active-region medium v = c/n r : We can write this equation in a form widely used in the literature [10,16]: where τ = 1/(2ω i ) is the photon cavity lifetime and Γ is called the confinement factor. However, in the literature the confinement factor is often defined as [16,28]: The difference, i.e. the lack of the medium permittivity ε = n 2 r usually does not make a huge difference, but equation (17), unlike (18), has a clear physical sense-this is the ratio of the number of photons interacting with the gain medium to the total number of photons in the resonator.

SSM response model
Our model of SSM response takes into account three phenomena: (a) electric capacitance-related phenomena; (b) interaction between the carriers and photons in the active region; and (c) the non-zero photon cavity lifetime.
We assumed that the changes in the temperature in the laser caused by the current modulation are negligible, because thermal processes are so slow that there is no significant modulation of the temperature when the current modulation frequencies are of the order of 1 GHz or higher. As a result, all the temperature-dependent parameters used in this model are calculated at the temperature calculated in the CW model for a given current.
The first element of the model is calculated with the aid of our model of capacitance [20]. The second module is based on well-known rate equations, however our model is different from the standard text-book approach in important details. The rate equations have been used to describe relaxation oscillation or modulation response in semiconductor lasers since the 1970s [29][30][31]. In general, the equations are non-linear and their solutions have to be found numerically. In the considered case of small-signal modulation it is very convenient to use a linear approximation of the equations, because their analytical solution is known. Although the equations and their solutions have been known for many years, we added to this solution a description of the fact that the instantaneous number of emitted photons is not exactly proportional to the instantaneous number of photons in the active region. To the best of our knowledge, this has not been done before. Additionally, we take into account the impact of the spectral hole burning (SHB) in a different way than is generally presented in the literature.
We present our analysis for a single-mode case, starting from the model of carrier-photon interactions. The starting point is a set of two differential equations, the same as those used in the literature [10,16,32]: where N is the number of carriers in the active region, I is the current injected into the active region, R are the carrier losses in all the processes other than stimulated emission. P is the number of photons in the whole cavity, and L stands for the photon losses. As defined in the previous section, G is the rate of stimulated emission. In this convention, all the rates on the right-hand side are positive and have a unit of s −1 . In the analysis of small signal modulation, we assume that the driving current is a DC value I 0 with a sinusoidal modulation δI(t) = B I sin(Ωt) with a small amplitude. Then, both N(t) and P(t) can be presented as an average value modified by a small modulation: We assume that the current modulation amplitude is small enough that we can assume that all the rates in the differential equations are linear functions of N and P. Then, equations (19) lead to a system of equations where there are only the amplitudes of the oscillating functions: This is a non-homogeneous system of two linear differential equations that can be written in a matrix form: Because the current oscillation δI is sinusoidal, so are oscillations δN and δP. They have the same frequency, but different amplitudes and phases. We are mostly interested in the amplitude of the oscillation of the number of photons in the cavity P. This system of equations can be easily solved analytically, and the amplitude oscillation of the number of photons is as follows: where Ω is the angular frequency of the modulation and γ is the square matrix from equation (24). So far, our calculations are equivalent to those in the textbooks. However, in the literature, the amplitude of the number of photons is usually identified with the amplitude of the output power [10,13,16,32]. This assumption can be false (but it depends on the modulation frequencies one wants to consider), because the cavity in a VCSEL may have a high quality factor. This effect can be easily taken into account. The photons in the cavity decay exponentially, so the density of the probability that a photon generated at time 0 will be in the cavity at time t is described by the following formula: The number of photons emitted P em (or simply the instantaneous output power) is then given by the following convolution: Because we assume that P(t) is a sinusoidal function, this convolution can be calculated analytically. As a result, we get a simple formula for the amplitude of the output power oscillation, assuming a constant amplitude of the current modulation: The second factor in this equation, describing the impact of the non-zero photon lifetime, has the same form as the factor presented in [16,32] as the impact of the laser's parasitics. In our model, the parasitics are taken into account by means of a model for capacitance-related phenomena developed recently by our group. Details of this model can be found in [20]. Using this model, we can find the active (i.e. not apparent) current flowing through the active region as a function of the modulation frequency. In other words, assuming that the voltage modulation amplitude is the same for all the frequencies (as it is the case in SSM experiments), we can calculate the function B I (f ) which is the amplitude of the modulation of the current injected into the active region. Then, the full formula for the amplitude of the emitted power in the SSM experiment is as follows: where f is the modulation frequency. The formula we obtained is not significantly different from what can be found in the literature, however we take into account the reduction in the oscillation amplitude related to the cavity lifetime and use a separate model for capacitancerelated phenomena. Our main goal, however, was not to rely on fitting the parameters of this model to dynamic measurements. In order to achieve this goal, we have to calculate the elements of matrix γ from equation (24). It is worth emphasizing that thanks to this fundamental difference it is possible to predict the performance of a laser already at the designing stage.

Determination of matrix γ
The physical sense of the elements in matrix γ is given by equations (23) and (24), namely: All the derivatives are calculated at N 0 , P 0 , which are the number of carriers in the active region and the number of photons in the cavity in a laser driven by a constant current I 0 . These values are found using our CW above-threshold model described before.
Here we present how the derivatives are calculated. Derivative ∂L ∂P is pretty straightforward. The losses of the mode are proportional to the number of photons it contains, so this derivative is the proportionality constant, in this case (see equation (2)): In order to calculate the derivatives of the recombination loses R we need to define this function in our model. It describes all the carrier losses other than those caused by stimulated emission and is calculated using the following formula: And the number of carriers N is related to the carrier concentration u in the following way: Because none of the terms in equation (34) depends explicitly on P, we see that: The derivative ∂R ∂N is not properly defined in this picture, because R as presented in equation (34) is not a function of N because R depends on the distribution of carriers, not only on their number. In our case, however, when the possible carrier distributions under CW operation are restricted to those given by the solution of equation (5), we know how to calculate R(N 0 ). We approximate the derivative by a difference quotient: In order to calculate R(N 0 + δN, P 0 ) we need to know the distribution of δN, i.e. a function δu(r) such that Then Distribution δu has to be associated with the shape of the current density injected into the active region. In the CW regime the shape of the carrier concentration is calculated using the diffusion equation (5). Now, however, we are interested in the change in the shape caused by a quickly changing current of a small amplitude. There are two extreme cases that we can consider: 1) the diffusion processes are fast enough to spread the carriers injected by the modulation current; or 2) the diffusion is too slow to introduce any noticeable changes to the injection profile. The range of the carrier spreading can be estimated from above by the spreading of a Dirac-delta-like impulse i.e. by a half of the FWHM of the Green's function for the diffusion equation, equal to: where D is the diffusion constant. In our simulations we solve diffusion equation (5) for D = 10 cm 2 s . For t = 0.05 ns, being the current rise time for modulation frequency of 10 GHz, the spreading radius is equal to only around 0.26 μm. In our simulations we consider frequencies up to 40 GHz, so the second assumption presented above is much better than the first one. We also verified this conclusion by performing S 21 simulations in both cases and we got a much better comparison with experiment in the second case.
We use the same δu distribution to calculate ∂G ∂N as an appropriate difference quotient, where G(N 0 , P 0 ) and G(N 0 + δN, P 0 ) are calculated using equation (4). In this equation, the only function that depends on u is the gain g.
The last remaining derivative is ∂G ∂P . As shown in the description of the CW model, function G can be expressed as: The dependence of the gain on P is a result of spectral hole burning. This phenomenon is decribed in the CW model by the following formula (see equation (6) and (10)) There is a question whether this equation can be used when P varies with frequencies of tens of GHz. If the time needed to deplete the lasing electronic states is significantly lower thañ 10 ps (half of the period for frequency of 50 GHz), then the above equation can be used to calculate ∂g/∂P. However, if after a time comparable with a half of the oscillation period the depletion is only partial, the constant α in equation (42) should be reduced accordingly to a value α mod . This reduction should increase with the modulation frequency. In the presented dynamic model we wanted to avoid any parameters that would have to be fitted. Like in the case of the shape of function δu discussed above, there are two extreme cases: 1) α mod = α; or 2) α mod = 0. Both assumptions allows us to avoid using of a parameter whose value we are not able to determine theoretically, but it is rather unlikely that any of them can be true in the whole spectrum of considered frequencies. In order to decide if any of them allows us to obtain reasonable results, we compared simulated S 21 (f ) curves, calculated using each of the two considered assumptions, with experiment. The functions obtained using the value of α mod = α CW (i.e. using the values we used in modelling of CW LI curves) did not resemble measured S 21 curves, so we decided to use the approximation α mod = 0, which gave us good results as it is shown further in this paper. With this assumption ∂g/∂P = 0, and: Collecting the results presented above, we can write the elements of matrix γ in a form reduced when compared to equations (31) and (32): where the two partial derivatives are calculated numerically in the way described above, and their values do not depend on any parameters other than those used in the CW model. It is worth mentioning that the calculations related to the dynamic model are not time-and memory-consuming compared to the static calculations. In the implementation used for simulations presented in this paper, the 3D vectorial optical model needed by far most of the resources (however less sophisticated, much faster, models can be used). Next were the iterations of the diffusion equation (the loop in the flow-chart in figure 1). Because the SSM model is based on rate equations which have analytical solutions, this part does not increase the computation time in any noticeable way.

Limitations of the model
The model presented here can be applied to the case of small signal modulation. Otherwise the rate equations are no longer linear and their solution cannot be obtained, in general, analytically and their parameters will not be given as the partial derivatives we used here. In the form presented here the model describes only sinusoidal voltage modulation, but with the aid of a Fourier transform could be used to simulate eye diagrams as long as the assumption of small amplitude of the signal is valid. This is usually not the case in eye-diagram experiments, but nevertheless such simulations can be very useful. This subject, however, is beyond the scope of this work.
Another important constraint of the version of the model presented in this paper is the assumption of single-opticalmode operation. This is not a fundamental limitation. The rate equations can be used to take into account several modes. The CW model can also be used to simulate multi-mode operation, but then the numerical complexity increases greatly.
The presented model can be used only if one can assume stable modal operation. For the most common type of semiconductor laser, the EEL, this is not the case, except for sophisticated constructions such as the DFB laser. In VCSELs, however, this assumption is generally true. Additionally their dynamic properties are of great importance. For this reason, this model is presented as a model for VCSELs, although after slight modifications resulting from the different light propagation direction in EELs it should be suitable to DFB lasers. Because the rate equations have been successfully used to model ordinary EELs it could be possible to find effective parameters describing the whole group of the modes the laser switches between. This subject, however is not within the scope of this paper.

Results and comparison with experiment
In order to verify our model, we compared the results of simulations with measured characteristics of gallium-arsenidebased VCSELs emitting at~980 nm. The processing of the VCSELs was performed at the Technical University Berlin using planar procedures developed there. First, top metal pcontact rings were patterned on the surface of the sample using photolitography. Next the Au/Zn/Au contacts were laid down using thermal evaporation. In the following step a photolitography was used to define mesas. The mesas were then dryetched in an inductively-coupled plasma reactive-ion etching (ICP-RIE) reactor using Cl/BCl 3 plasma. The high-Al-content layers placed close to the active region were then oxidized in an oxidation oven in 420 • C at 50 mbar in a water vapor atmosphere in order to form the oxide apertures. The photolithography was then used again to define second, large diameter mesas. The second mesas were dry-etched all the way down through the bottom DBR to the ∼ 1.5 μm thick (n+)-doped GaAs buffer layer. In the next step the bottom Ni/AuGe/Au ncontacts were patterned using photolitography and deposited onto the 1.5 µ m-thick (n+) GaAs layer using thermal evaporation. In order to lower the contact resistance of the contacts, the sample was annealed at 420 • C for 1 minute using a rapid thermal annealing furnace. The sample was then planarized using a photo-sensitive BCB spin-on film, which was selectively removed in order to expose the top mesas (and the top metal contacts) and the bottom contacts. In the final step, high-frequency ground-signal-ground (GSG) Ti/Au pads were patterned using photolitography and deposited using thermal evaporation.
The static measurements involved light-voltage-current measurements performed with use of an integration sphere with an InGaAs photodetector and optical emission spectra versus bias current measurements performed using an optical spectrum analyzer with a 0.02 nm resolution. The high-speed SSM measurements were performed with use of a Hewlett-Packard 8722 C vector network analyzer (VNA). Each VCSEL was biased with a DC bias current from a current source which was fed to the VNA via a bias-tee and mixed with a small amplitude sinusoidal AC current. The AC part was modulated with frequencies from 0.05 to 40 GHz. The signal was sent from the VNA port 2 through a transmission line to a high speed GSG probe which was connected to the VCSEL. The light from the VCSEL was collected by a 62.5 μm core-diameter OM1 multimode lensed optical fiber and fed to to New Focus 25 GHz photodetector module whose output is connected to the VNA's port 1. Both the static and dynamic measurements were performed at 25 • C heat-sink temperature.
In the simulations, as a first step, we modelled the laser's LIV curves. The gain model is calibrated using photoluminescence (PL) spectra from a test wafer, where the spectrum is not modified by the presence of an optical cavity. Our model is able to calculate spectra of both stimulated and spontaneous emission. By adjusting the parameters of the active region we obtained a simulated PL spectrum with the peak width and position matching the experiment. The same parameters were then used to simulate optical gain in the laser. The resistance of the electric contacts and the lateral electric conductivity of the DBRs were measured at LUT using the transmission line method (TLM). All other experimental results were obtained at TU Berlin.
One of the parameters of the structure that significantly influences both the static and dynamic parameters is the oxide aperture diameter. In the model, its impact is taken into account in may ways. The aperture diameter impacts the VCSEL's: resistance, capacitance (as it is related to the surface of the oxide insulation), and the shape and the loses of the optical mode. The aperture diameter is not precisely known. Our estimation is based on the oxidation rate measured on calibration structures and the fact that the simulated lasers operated on a single optical mode. We estimate the uncertainty to be around ±1 μm. In the simulation we used the value of ϕ = 3.6 μm. Another important parameter is the Auger recombination coefficient C in equation (5). Its actual value can depend on, for instance, the parameters of the QWs, so there is no universally correct value of this parameter. In the simulations we assumed C = 3 · 10 −29 cm 6 /s at 300 K with a linear dependence on temperature dC/dT = 3 · 10 −31 cm 6 /(sK). The value assumed lies within the range of the values for active regions emitting at 980 nm reported in the literature [33,34]. The value of the output power depends also on material absorption in the layers of the laser. For some layers, for instance the gradient layers in the DBRs, there is little information on the absorption, so the values used in the simulations, based on interpolation, can be inaccurate. Besides, there can be some imperfections in the fabrication that are not taken into account in the simulations. We think that these can be the reasons why the simulated LI curve in figure 2 is almost twice as high as the measured one. At this stage, we fitted the value of the parameter α in equation (6), which in our model, is responsible for the LI roll-over. For α = 3.6 · 10 −16 cm 3 we got the shape (not the actual values, though) of the simulated curve very similar to the experiment, as it can be seen in figure 2. The value of parameter α is difficult to determine theoretically. The strength of SHB effect depends on factors such as (among others) carrierphonon interaction. In principle, it is possible to expand the rate equations to model SHB [35], but then new parameters would be necessary, and the complexity of the whole model would increase greatly. Instead, we decided to describe this process in a simplified way through the parameter α whose value is fitted to CW characteristics. Then, we can avoid fitting any parameters in the rate equations.
The static simulations are also used to determine all the parameters for the dynamic model as described in section 2.2, so no fitting to the experiment is performed when the static model parameters are set. In figure 3 the simulated and measured SMM response (S 21 ) for the modulation frequency up to 40 GHz are presented. The frequency response curves are obtained for above-threshold currents from 1 mA up to 5 mA.
The model correctly reproduces the following important features (see also figure 4): • With increasing current the resonance peaks shift toward higher frequencies but this effect is saturated at the current around 5 mA. • The height of the peaks decreases with current, while their width increases. Except for the lowest currents, the simulated widths and heights are in good agreement with experiment. • The slope of the decreasing part of the resonance curve is modelled correctly.
The most apparent difference between the simulated and measured curves is the fact that the resonant peak frequency f r and the −3 dB frequency f 3 dB are not the same. Figure 4 shows more clearly where this difference comes from. In the experiment, f 3 dB is higher at low currents, but the slope df 3 dB /dI is very similar to the modelled one. Near the current of 2 mA, the experimental f 3 dB (I) curve bends, unlike in the model. As a result, both curves intersect at a current between 2 mA and 3 mA. Both curves reach their maximal value of similar currents between 5 mA and 6 mA.

Summary
We have presented a numerical model of SSM response curves for VCSELs. The main part of this model is based on the well known rate equations, but the parameters are not fitted, but rather they are calculated based on the results of CW simulations of the laser. This significant improvement allows for determination of the impact of different modifications of the structure even at the designing stage.
A comparison with experimental data showed a good qualitative agreement, although the model predicts slightly higher values of the maximal frequency f 3 dB . We have shown the usefulness of this model by simulating the impact of the number of the active-region QWs on the modulation properties. We have shown that in the analyzed laser design having more QWs in the active region would be beneficial.
We have also analyzed what information on the nonradiative losses and stimulated emission rate can be obtained from an analysis of the parameters of the resonant peak as a function of the driving current.

Acknowledgment
Funding: this work has been partially supported by the Polish National Science Centre, grant no. 2016/21/B/ST7/03532.

Appendix: Analysis of the parameters of the model
Having shown that our model simulates correctly the most important features observed in the SSM experiment, we can focus on understanding how the parameters of the model change with changing conditions. In figure A1 recombination rates and the number of carriers in the whole active region (QWs) are presented as functions of the driving current. The number of carriers in the active region of a lasing laser is such that the resulting gain balances all the losses, exactly like at threshold. The surplus injected carriers are turned into photons. The recombination rate denoted as 'losses' is given by the integral of in equation (34), while the stimulated emission rate is given by equation (4). The rates sum up to a linear function of the current, namely the number of carries injected into the active region per a unit of time. The shape of the stimulated-emission curve is nearly identical with the shape of the LI curve (see figure 2). The possible differences may come from the fact that the absorption coefficients for the laser materials and the mirror reflectivity are temperature and wavelength-dependent, so they depend on the driving current. The losses are initially (up to around 3.5 mA) almost constant (and even slightly decreasing), because the number of carriers is not changing significantly. The initial small increase in the number of carriers is caused by the generation of carriers in the distant (from the aperture) region of the active region due to absorption of the mode occurring there. At low currents, where the temperature rise is small, the mode is not strongly confined (as it can be seen in figure A2), and overlaps to a certain degree with the high absorption outside the aperture.
The drop in the number of carriers observed further is a result of a slight decrease in the modal losses caused by the observed thermal focusing and the fact that the gain spectrum, that at RT is blue-shifted relative to the modes wavelength, red-shifts with temperature faster than the emission. All these effects, however, are relatively subtle. The situation changes dramatically when the SHB effect, increasing with the number of photons, becomes significant. It reduces the gain and in order to overcome this decrease, more carriers are necessary. At higher carrier concentrations and temperatures, Auger recombination (with the rate of Cu 3 ) becomes dominant, what explains the superlinear increase in the losses. Figure A2 shows also the effect of spatial hole burning, namely the  decrease in the carrier density at the center, where the mode intensity is strongest, relative to the density near the aperture edge, where the current injected is highest and the mode intensity is smaller than at the center.
The simulated SSM response spectrum S 21 (f ) is a product of three functions as shown in equation (29). In figure A3 the S 21 for the driving current of 3 mA and the three factors are shown.
The general shape of the response spectrum is determined by the photon-carrier interactions described by function B P . The parasitic effects, i.e. function B I do have an important impact on the height of the resonance peak and the frequency f 3 dB . In the case of this laser, where the top DBR has only 14.5 pairs, the impact of the cavity-lifetime (function B C ) is very weak in this range of frequencies. This, however, may not be true in the case of VCSELs having high-Q-factor cavities. For example, calculations for VCSELs with 22.5 and 26.5 top (coupling) DBR pairs, otherwise identical to the VCSEL simulated here, show that the reduction in the frequency response caused by the cavity-lifetime is 5.5% and 8.5% respectively, at 40 GHz.
In order to analyze the behaviour of function B P we plotted, in figure A4, the derivatives that define this function as described in equations (25) and (44).
The term ω i is almost constant and its value is over 10 times higher than the other two. This high value is a result of the low quality-factor of the cavity of the analyzed laser. Derivative ∂G ∂N is strongly related to the averge number of photons in the cavity (see equation (41)), so the shapes of this function and the LI (figure 2) curve are similar. The last derivative, ∂R ∂N is much smaller than ∂G ∂N for currents other than those near the threshold and the roll-over, because in this region losses R are much smaller than the photon generation G as presented in figure A1.
The position of the resonance peak, i.e. the frequency f r is approximately equal to the position of the maximum of function B P , because the other two factors in equation (29) are slowly varying functions. By finding the minimum of the denominator in equation (25) we can get a fairly simple formula: The right-hand side of this formula can be negative, which means that the resonant peak is not present. This may happen in two cases: (a) the derivative ∂G ∂N is close to 0; and (b) the sum ∂G ∂N + ∂G ∂N is big enough compared with ω i . This can be approximated by the condition that ∂G ∂N is at least twice as big as ω i when ∂G ∂N ≫ ∂R ∂N which should be fulfilled whenever a VCSEL operates far enough from its threshold. The first condition must be satisfied near the threshold, but when, like in the analyzed VCSEL, ω i ≫ ∂R ∂N the peak appears almost immediately above the threshold. The second condition can or cannot be fulfilled, depending on the cavity Q-factor. VCSELs with Q-factors high enough (with ω i low enough) may have response curves without the resonance peak. A transition from S 21 curves with a peak to curves without a peak caused by changes in the Q-factor has been shown experimentally in [36], confirming the theoretical analysis presented here. In the laser considered in this paper, however, the Q-factor is so low that all the S 21 have a resonant peak.
The height of the resonance peak can be roughly approximated (neglecting the impact of functions B I and B C ) by the maximum (if present, as discussed above) of the function B 2 P (Ω)/B 2 P (0). This is described by the following formula: Function H has a maximum when the following condition: is satisfied. Derivative ∂G ∂N starts from 0 at the threshold and rises with the current. Derivative ∂R ∂N is initially approximately constant (see figure A4) and its value is related to the recombination itself. Because the recombination rate is a superlinear function of carrier concentration, a higher value of the derivative means a higher value of the recombination rate. As a result, the current I H for which the S 21 has the highest peak is an indication of the level of unwanted carrier recombination in the active region. The higher the current is, the lower this recombination is. Because the condition (A3) does not contain ω i , the value of I H (compared to the threshold current) can be a useful qualitative measure of the recombination losses of the carrier in the active region. In figure A5 the results of the simulations of the height of the resonant peak (including the whole formula (29)) and its spectral position are shown. When comparing the dashed line in this figure with figure 3 one can notice that our simulations probably underestimated the rate of the non-radiative recombination, because in the measurements the highest resonant peak appears at a current somewhere in the vicinity of the current of 2 mA, while in the simulations it is slightly above 1 mA. The formula (A1) for the resonant frequency can be radically simplified under the condition ω i ≫ ∂G ∂N ≫ ∂R ∂N : .
Although this formula is an approximation and its validity is limited by the conditions discussed above, it gives an idea how f r and hence f 3 dB can be increased. Frequency f 3 dB itself, unlike f r , cannot be reliably estimated using equation (25), because the other two factors in equation (29) do impact the result. It seems that, in order to increase f r , one should increase ω i by reducing the number of DBR pairs (reducing the Qfactor) and increase ∂G ∂N by, for example, increasing the number of photons in the cavity. However, if the number of photons is so high that the spectral hole burning becomes an important factor, formula equation (A4) is no longer valid, and, as presented in figures 3 and 4, f r is saturated. The other way to increase ∂G ∂N could be increasing the differential material gain ∂g ∂N . Its highest value is near the transparency conditions. A laser cannot operate with a transparent active region, but having as low as possible threshold gain should be beneficial in this respect. But this contradicts the condition of having a low-Q cavity (high ω i ). Moreover, as it can be seen in equation (A2), with higher ω i we get higher resonant peaks which may have an adverse impact on large-signal modulation patterns. The considerations presented above show that a proper choice of the number of QWs can help to increase the modulation bandwidth. Numerical simulations can be extremely useful in analyzing such questions, if a proper model is used.
We performed simulations of the impact of the number of QWs in the active region of the VCSEL under consideration. Apart from the original structure with 5 QWs, we calculated lasers with 3, 4 and 6 QWs. The only part of the laser that has been changed is the QWs, their barriers and the two adjacent layers. The thickness of the adjacent layers has been adjusted in such a way that all four structures emitted at the same wavelength. In each case, the QWs were placed symmetrically on both sides of the anti-node of the standing way, so as to maximize the confinement factor. In figure A6 f 3 dB as a function of the current, for the four considered active regions, are presented. Because the analyzed lasers have the top DBRs of only 14.5 pairs, the losses are high, and a higher number of the QWs improves the laser's performance. The cavity losses do not depend on the number of the QWs, so the total gain necessary to support lasing, provided by all the wells, is the same in all the cases. When there are more QWs, the gain per a single well is lower, and as a result derivative ∂g ∂N is bigger. As discussed above, this increases the resonance frequency and hence f 3 dB . The standing wave electric field intensity in the two outermost wells is still relatively high, around 75% of its highest value in the active region, so it is possible that 7 QWs could perform even better. However, in the analyzed structure design, maximally 6 QWs can be put in the active region without changing further layers in the resonator, which would make the results difficult to compare.