Document [original]

Journal of Physics D: Applied Physics

J. Phys. D: Appl. Phys. 53 (2020) 345101 (13pp) https://doi.org/10.1088/1361-6463/ab8b94

Numerical model for small-signal

modulation response in vertical-cavity

surface-emitting lasers

Michał Wasiak1, Patrycja

Spiewak1, Nasibeh Haghighi2, Marcin Ge¸bski1,2, Emilia

Pruszy´

nska-Karbownik1, Paulina Komar1, James A Lott2and Robert P Sarzała1

1Institute of Physics, Lodz University of Technology, Ł´

od´

z, Poland

2Institute of Solid State Physics and Center of Nanophotonics, Technical University of Berlin, Berlin,

Germany

E-mail: [email protected]

Received 20 January 2020, revised 4 April 2020

Accepted for publication 21 April 2020

Published 8 June 2020

Abstract

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.

Keywords: small-signal modulation response, semiconductor lasers, numerical modeling, optical

data transfer

(Some figures may appear in colour only in the online journal)

1. 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–6], approaching even

100Gbps [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 val-

ues can be chosen arbitrarily to a certain extent, or have to be

fitted to experiment [10–16]. These models have shown that

the approach based on rate equations can give good agree-

ment with the small-signal modulation (SSM) response exper-

iment. This approach allows for determination where the lim-

itations of the obtained modulation bandwidth come from, but

not always how to overcome the problem. The greatest chal-

lenge is to calculate, not fit, the parameters in these equations

in such a way that the results of SSM experiments can be pre-

dicted 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 modu-

lation, a reliable model of above-threshold continuous-wave

This is an open access article distributed under the terms of the Creative Commons Attribution 4.0 License.

Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI.

J. Phys. D: Appl. Phys. 53 (2020) 345101 M Wasiak et al

(CW) operation of the VCSEL is needed. In this paper we

show that based on such a model developed by the Photon-

ics 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–

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 character-

istics of a GaAs-based high-speed VCSEL.

2. 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 mod-

els 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. Sec-

tion 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 usu-

ally 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.

2.1. CW operation model

The numerical model of CW above-threshold operation is

based on models of the following physical phenomena:

Figure 1. Flowchart of the calculations in the CW VCSEL

operation model.

•a self-consistent model of the electrical and thermal phe-

nomena in the laser

•a model of carrier diffusion and recombination in the active

region

•an optical model for the electromagnetic radiation in the

laser cavity

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 memory-

consuming 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 temper-

ature 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:

ω=ωr+iωi.(1)

In the case of a lossy cavity, the imaginary part ωi>0. This

means, that the mode’s energy in the cavity EMwill decay

according to the following formula:

EM(t) = EM(0)exp(−2ωit)(2)

so the number ωidescribes all the mode losses. In our case,

where the gain in the active-region is 0, it includes optical

J. Phys. D: Appl. Phys. 53 (2020) 345101 M Wasiak et al

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 distri-

bution 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:

G+dEM

dt =0.(3)

The number of photons generated in the active region is given

by the following equation:

G=ˆAct

cε0nrg(r)

2ℏωr∥E(r)∥2dxdydz (4)

where Act denotes the active region, r= (x,y,z),cis the speed

of light, ε0is the dielectric constant, nris the refractive index

(it also can depend on r), gis the material gain, ℏωris the

photon energy, and Eis the mode’s electric field in the stand-

ard complex-number notation. The impact of the stimulated

emission on the carrier (and hence gain) distribution, mani-

festing as the spatial hole burning, is taken into account by

the presence of the stimulation-emission losses in the carrier

diffusion equation:

D∆u(r)−(Au(r) + Bu2(r) + Cu3(r))

−cε0nrg(r)

2ℏωr∥E(r)∥2+j⊥(r)

ed =0 (5)

whereuisthecarrier densityin theactive region,Disthe ambi-

polar diffusion coefficient, A,B,Care, in general temperature-

dependent, mono-molecular, bi-molecular and Auger recom-

bination coefficients [27]; j⊥is the perpendicular component

of the current density injected into the active region, dis the

total thickness of the active-region QWs, and eis the element-

ary charge. In our implementation, this equation is solved in

2D, and the mode’s distribution is averaged in the vertical dir-

ection in the QWs. The optical gain distribution is calculated

based on the carrier, temperature and mode distributions. First,

optical gain ˜

gis calculated assuming thermal (Fermi–Dirac)

electronic states occupancy. Our simulations showed that such

a gain model is not sufficient to describe the observed rolling-

over 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]:

g=˜

1+αε0n2

2ℏωr∥E∥2.(6)

The factor ε0n2

2ℏωr∥E∥2is the density of photons that interact with

the carriers. Because only the electric field of the electromag-

netic 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 m3and describes

the strength of SHB in this system. In our calculations, we fit-

ted αto measured CW LI curves.

The total energy of the mode is given by the following for-

mula:

EM=ε0

2×ˆRes (n2

r(r)∥E(r)∥2

2+η2

0∥H(r)∥2

2)dxdydz (7)

where Res denotes the entire resonator, η0is the free space

impedance and His the mode’s magnetic field. Although loc-

ally the energy densities of the electric and the magnetic fields

are usually different (in a standing wave the nodes of the mag-

netic 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:

EM=ε0

2ˆRes

r(r)∥E(r)∥2dxdydz

=µ0

2ˆRes

∥H(r)∥2dxdydz.(8)

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 dens-

ity of the electric field. Because the optical model (Maxwell’s

equations) makes no distinctions between solutions that dif-

fer in the amplitude only, we can use in the calculations elec-

tromagnetic 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 MEand MH. Then:

ℏωr=ε0

2ˆRes (n2

r(r)ME(r)

2+η2

MH(r)

2)dxdydz.(9)

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:

E2(r) = PME(r)H2(r) = PMH(r)(10)

where F2stands for ∥F(r)∥2for Fbeing either Eor H. With

this notation, we can write equation (3), using equation (2), in

the following form:

PˆAct

cε0nrg(r)

2ℏωrME(r)dxdydz =2ωiP.(11)

Although the amplitude Pcancels out on both sides, the left-

hand side does depend on P, because the carrier density and

hence material gain gdo. The above equation can be written

as:

ˆAct

cε0nrg(r,P)

2ℏωrME(r)dxdydz =2ωi.(12)

There is a fundamental question, whether this equation has

solutions. The optical gain gis a decreasing function of P

J. Phys. D: Appl. Phys. 53 (2020) 345101 M Wasiak et al

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ε0nrg(r,0)

2ℏωrME(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 presen-

ted in this section, we can write an equation for the average

threshold gain:

⟨gth⟩c´Act nrE2

´Res n2

rE2=2ωi.(14)

Or, using the light speed in the active-region medium v=c/nr:

⟨gth⟩v´Act n2

rE2

´Res n2

rE2=2ωi.(15)

We can write this equation in a form widely used in the liter-

ature [10,16]:

⟨gth⟩vΓ=1

τ(16)

Γ=´Act n2

rE2

´Res n2

rE2(17)

where τ=1/(2ωi)is the photon cavity lifetime and Γis called

the confinement factor. However, in the literature the confine-

ment factor is often defined as [16,28]:

Γ=´Act E2

´Res E2.(18)

The difference, i.e. the lack of the medium permittivity ε=n2

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.

2.2. SSM response model

Our model of SSM response takes into account three phenom-

ena:

(a) electric capacitance-related phenomena;

(b) interaction between the carriers and photons in the active

region; and

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 mod-

ulation of the temperature when the current modulation fre-

quencies 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 oscilla-

tion or modulation response in semiconductor lasers since the

1970s [29–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 num-

ber of photons in the active region. To the best of our know-

ledge, 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]:

N(t) = I(t)

q−R(t)−G(t)(19a)

P(t) = G(t)−L(t)(19b)

where Nis the number of carriers in the active region, Iis

the current injected into the active region, Rare the carrier

losses in all the processes other than stimulated emission. P

is the number of photons in the whole cavity, and Lstands for

the photon losses. As defined in the previous section, Gis 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 I0with a sinusoidal modulation

δI(t)=BIsin(Ωt) with a small amplitude. Then, both N(t) and

P(t) can be presented as an average value modified by a small

modulation:

I(t) = I0+δI(t)(20)

N(t) = N0+δN(t)(21)

P(t) = P0+δP(t).(22)

J. Phys. D: Appl. Phys. 53 (2020) 345101 M Wasiak et al

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 Nand P. Then, equations (19)

lead to a system of equations where there are only the amp-

litudes of the oscillating functions:

δN=δI

q−(∂R

∂N(N0,P0)δN+∂R

∂P(N0,P0)δP)

−(∂G

∂N(N0,P0)δN+∂G

∂P(N0,P0)δP)(23a)

δP=(∂G

∂N(N0,P0)δN+∂G

∂P(N0,P0)δP)

−(∂L

∂N(N0,P0)δN+∂L

∂P(N0,P0)δP).

(23b)

This is a non-homogeneous system of two linear differential

equations that can be written in a matrix form:

dt [δN(t)

δP(t)]=[γ11 γ12

γ21 γ22][δN(t)

δP(t)]+[δI(t)

0.].(24)

Because the current oscillation δIis sinusoidal, so are oscilla-

tions δNand δP. They have the same frequency, but different

amplitudes and phases. We are mostly interested in the amp-

litude 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:

BP(Ω) = 1

√(Ω2−detγ)2+ Ω2tr2γ

(25)

where Ωis the angular frequency of the modulation and γ

is the square matrix from equation (24). So far, our calcula-

tions 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 mod-

ulation 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 tis described

by the following formula:

pem(t) = {0 for t<0

2ωiexp(−2ωit)for t>0.(26)

The number of photons emitted Pem (or simply the instantan-

eous output power) is then given by the following convolution:

Pem(t) = P(t)∗pem(t).(27)

Because we assume that P(t) is a sinusoidal function, this con-

volution can be calculated analytically. As a result, we get a

simple formula for the amplitude of the output power oscilla-

tion, assuming a constant amplitude of the current modulation:

Bem(Ω) = BP(Ω) 1

√1+(Ω

2ωi)2=BP(Ω)BC(Ω) (28)

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 mod-

ulation 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

BI(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:

Bres(Ω) = BI(Ω)BP(Ω) 1

√1+(Ω

2ωi)2(29)

S21(f) = B2

res(2πf)

res(0)(30)

where fis 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 capacitance-

related phenomena. Our main goal, however, was not to rely

on fitting the parameters of this model to dynamic measure-

ments. In order to achieve this goal, we have to calculate the

elements of matrix γfrom equation (24). It is worth emphas-

izing that thanks to this fundamental difference it is possible

to predict the performance of a laser already at the designing

stage.

2.3. Determination of matrix γ

The physical sense of the elements in matrix γis given by

equations (23) and (24), namely:

γ11 =−∂R

∂N−∂G

∂Nγ12 =−∂R

∂P−∂G

∂P(31)

γ21 =∂G

∂Nγ22 =∂G

∂P−∂L

∂P.(32)

All the derivatives are calculated at N0,P0, which are the num-

ber of carriers in the active regionand the number of photons in

the cavity in a laser driven by a constant current I0. These val-

ues are found using our CW above-threshold model described

before.

J. Phys. D: Appl. Phys. 53 (2020) 345101 M Wasiak et al

Here we present how the derivatives are calculated. Deriv-

ative ∂L

∂Pis pretty straightforward. The losses of the mode are

proportional to the number of photons it contains, so this deriv-

ative is the proportionality constant, in this case (see equa-

tion (2)):

∂L

∂P(N0,P0) = 2ωi.(33)

In order to calculate the derivatives of the recombination loses

Rwe need to define this function in our model. It describes all

the carrier losses other than those caused by stimulated emis-

sion and is calculated using the following formula:

R=ˆAct (Au(r) + Bu2(r) + Cu3(r))dxdydz.(34)

And the number of carriers Nis related to the carrier concen-

tration uin the following way:

N=ˆAct

u(r)dxdydz.(35)

Because none of the terms in equation (34) depends explicitly

on P, we see that:

∂R

∂P(N0,P0) = 0.(36)

The derivative ∂R

∂Nis not properly defined in this picture,

because Ras presented in equation (34) is not a function of

Nbecause Rdepends on the distribution of carriers, not only

on their number. In our case, however, when the possible car-

rier distributions under CW operation are restricted to those

given by the solution of equation (5), we know how to cal-

culate R(N0). We approximate the derivative by a difference

quotient:

∂R

∂N(N0,P0) = R(N0+δN,P0)−R(N0,P0)

δN.(37)

In order to calculate R(N0+δN,P0)we need to know the dis-

tribution of δN, i.e. a function δu(r)such that

δN=ˆAct

δu(r)dxdydz.(38)

Then

R(N0+δN,P0) = ˆAct [A(u+δu)(r) + B(u+δu)2(r)

+C(u+δu)3(r)]dxdydz.(39)

Distribution δuhas to be associated with the shape of the cur-

rent 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 cur-

rent 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:

FWHM(t) = 2√2log2Dt (40)

where Dis the diffusion constant. In our simulations we solve

diffusion equation (5) for D=10 cm2

s. For t=0.05ns, being

the current rise time for modulation frequency of 10GHz, the

spreading radius is equal to only around 0.26μm. In our sim-

ulations we consider frequencies up to 40GHz, so the second

assumption presented above is much better than the first one.

We also verified this conclusion by performing S21 simulations

in both cases and we got a much better comparison with exper-

iment in the second case.

We use the same δudistribution to calculate ∂G

∂Nas an

appropriate difference quotient, where G(N0,P0)and G(N0+

δN,P0)are calculated using equation (4). In this equation, the

only function that depends on uis the gain g.

The last remaining derivative is ∂G

∂P. As shown in the

description of the CW model, function Gcan be expressed as:

G=PˆAct

cε0nrg(r,P)

2ℏωrME(r)dxdydz.(41)

The dependence of the gain on Pis a result of spectral hole

burning. This phenomenon is decribed in the CW model by

the following formula (see equation (6) and (10))

g=˜

1+αε0n2

2ℏωrPME

.(42)

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 than

~ 10ps (half of the period for frequency of 50GHz), 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 presen-

ted dynamic model we wanted to avoid any parameters that

would have to be fitted. Like in the case of the shape of function

δudiscussed 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 theor-

etically, 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 S21(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 val-

ues we used in modelling of CW LI curves) did not resemble

measured S21 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:

J. Phys. D: Appl. Phys. 53 (2020) 345101 M Wasiak et al

∂G

∂P(N0,P0) = ˆAct

cε0nrg(r)

2ℏωrME(r)dxdydz

equation(11)

=2ωiequation(33)

=∂L

∂P(N0,P0).(43)

Collecting the results presented above, we can write the ele-

ments of matrix γin a form reduced when compared to equa-

tions (31) and (32):

γ11 =−∂R

∂N−∂G

∂Nγ12 =−2ωi(44a)

γ21 =∂G

∂Nγ22 =0 (44b)

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 com-

pared 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 soph-

isticated, much faster, models can be used). Next were the iter-

ations 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.

2.4. 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, ana-

lytically 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 sub-

ject, 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-optical-

mode 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 semi-

conductor laser, the EEL, this is not the case, except for soph-

isticated 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 effect-

ive parameters describing the whole group of the modes the

laser switches between. This subject, however is not within

the scope of this paper.

3. Results and comparison with experiment

In order to verify our model, we compared the results of sim-

ulations with measured characteristics of gallium-arsenide-

based VCSELs emitting at ~ 980nm. The processing of the

VCSELs was performed at the Technical University Berlin

using planar procedures developed there. First, top metal p-

contact 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 photolito-

graphy was used to define mesas. The mesas were then dry-

etched in an inductively-coupled plasma reactive-ion etching

(ICP-RIE) reactor using Cl/BCl3plasma. 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 atmo-

sphere in order to form the oxide apertures. The photolitho-

graphy 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 n-

contacts were patterned using photolitography and deposited

onto the 1.5 µm-thick (n+) GaAs layer using thermal evap-

oration. In order to lower the contact resistance of the con-

tacts, the sample was annealed at 420◦C for 1 minute using a

rapid thermal annealing furnace. The sample was then planar-

ized 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 amp-

litude sinusoidal AC current. The AC part was modulated with

frequencies from 0.05 to 40GHz. 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 mul-

timode 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 per-

formed 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

J. Phys. D: Appl. Phys. 53 (2020) 345101 M Wasiak et al

Figure 2. Experimental and simulated LI curves. The dashed line

represents the measured values of the power multiplied by 1.8.

photoluminescence (PL) spectra from a test wafer, where the

spectrum is not modified by the presence of an optical cav-

ity. 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 res-

ults 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 sur-

face 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 cal-

ibration structures and the fact that the simulated lasers oper-

ated 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 recom-

bination coefficient Cin 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 simu-

lations we assumed C=3·10−29 cm6/s at 300 K with a lin-

ear dependence on temperature dC/dT =3·10−31 cm6/(sK).

The value assumed lies within the range of the values for act-

ive regions emitting at 980 nm reported in the literature [33,

34]. The value of the output power depends also on mater-

ial absorption in the layers of the laser. For some layers, for

instance the gradient layers in the DBRs, there is little inform-

ation 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 2is almost twice

as high as the measured one.

Figure 3. Simulated (smoothed) and measured SSM response

curves for driving currents from 1 mA to 5mA. The dashed

horizontal line denotes the −3dB level.

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 cm3we 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 para-

meter αis difficult to determine theoretically. The strength of

SHB effect depends on factors such as (among others) carrier-

phonon 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 fit-

ting 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 3the simulated and meas-

ured SMM response (S21) for the modulation frequency up

to 40GHz 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 5mA.

•The height of the peaks decreases with current, while their

width increases. Except for the lowest currents, the simu-

lated widths and heights are in good agreement with exper-

iment.

•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

frand the −3 dB frequency f3 dB are not the same. Figure 4

shows more clearly where this difference comes from. In

J. Phys. D: Appl. Phys. 53 (2020) 345101 M Wasiak et al

Figure 4. Measured and simulated f3 dB as a function of the driving

current.

the experiment, f3 dB is higher at low currents, but the slope

df3 dB/dI is very similar to the modelled one. Near the current

of 2mA, the experimental f3 dB(I)curve bends, unlike in the

model. As a result, both curves intersect at a current between

2 mA and 3mA. Both curves reach their maximal value of sim-

ilar currents between 5mA and 6 mA.

4. 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 simu-

lations 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 qual-

itative agreement, although the model predicts slightly higher

values of the maximal frequency f3 dB. We have shown the use-

fulness 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 non-

radiative 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

Figure A1. The number of carriers in the active region (dashed line,

right axis) and recombination rates R(losses) and G(stimulated

emission) as functions of the driving current.

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 emis-

sion rate is given by equation (4). The rates sum up to a linear

function of the current, namely the number of carries injec-

ted 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 num-

ber of photons, becomes significant. It reduces the gain and

in order to overcome this decrease, more carriers are neces-

sary. At higher carrier concentrations and temperatures, Auger

recombination (with the rate of Cu3) becomes dominant, what

explains the superlinear increase in the losses. Figure A2

shows also the effect of spatial hole burning, namely the

J. Phys. D: Appl. Phys. 53 (2020) 345101 M Wasiak et al

Figure A2. Profiles of the carrier concentration in the active region

(solid lines) and corresponding mode profiles (dashed line of the

same color) for different values of the driving current. The mode

profiles are normalized so that their integral (in 2D, assuming

cylindrical symmetry) is equal to 1.

Figure A3. An S21 curve (see equation (30) for the relation between

S21 and the Bfunctions) for the current of 3mA (solid line) and the

three factors that constitute it (dashed lines).

decrease in the carrier density at the center, where the mode

intensity is strongest, relative to the density near the aper-

ture edge, where the current injected is highest and the mode

intensity is smaller than at the center.

The simulated SSM response spectrum S21(f) is a product

of three functions as shown in equation (29). In figure A3 the

S21 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 BP.

The parasitic effects, i.e. function BIdo have an important

impact on the height of the resonance peak and the frequency

f3 dB. In the case of this laser, where the top DBR has only

14.5 pairs, the impact of the cavity-lifetime (function BC) is

very weak in this range of frequencies. This, however, may not

be true in the case of VCSELs having high-Q-factor cavities.

Figure A4. Values of the derivatives from equation (44) in

logarithmic scale as function of the driving current.

For example, calculations for VCSELs with 22.5 and 26.5 top

(coupling) DBR pairs, otherwise identical to the VCSEL sim-

ulated here, show that the reduction in the frequency response

caused by the cavity-lifetime is 5.5% and 8.5% respectively,

at 40GHz.

In order to analyze the behaviour of function BPwe plot-

ted, in figure A4, the derivatives that define this function as

described in equations (25) and (44).

The term ωiis 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

∂Nis 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

∂Nis

much smaller than ∂G

∂Nfor currents other than those near the

threshold and the roll-over, because in this region losses Rare

much smaller than the photon generation Gas presented in

figure A1.

The position of the resonance peak, i.e. the frequency fr

is approximately equal to the position of the maximum of

function BP, 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 for-

mula:

4π2f2

r=detγ−tr2γ

2=2ωi∂G

∂N−(∂R

∂N+∂G

∂N)2

.(A1)

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

∂Nis close to 0; and

(b) the sum ∂G

∂N+∂G

∂Nis big enough compared with ωi. This

can be approximated by the condition that ∂G

∂Nis at least

twice as big as ωiwhen ∂G

∂N≫∂R

∂Nwhich should be fulfilled

whenever a VCSEL operates far enough from its threshold.

J. Phys. D: Appl. Phys. 53 (2020) 345101 M Wasiak et al

Figure A5. The height of the resonant peak (dashed curve, left axis)

and its position as functions of the driving current.

The first condition must be satisfied near the threshold, but

when, like in the analyzed VCSEL, ωi≫∂R

∂Nthe 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 ωilow enough)

mayhave response curveswithout the resonance peak. A trans-

ition from S21 curves with a peak to curves without a peak

caused by changes in the Q-factor has been shown experiment-

ally 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 S21 have a resonant peak.

The height of the resonance peak can be roughly approx-

imated (neglecting the impact of functions BIand BC) by

the maximum (if present, as discussed above) of the function

P(Ω)/B2

P(0). This is described by the following formula:

H=(4

ξ(4−ξ))2

(A2a)

ξ=tr2γ

detγ=(∂R

∂N+∂G

∂N)2

2ωi∂G

∂N

.(A2b)

Function Hhas a maximum when the following condition:

∂G

∂N=∂R

∂N(A3)

is satisfied. Derivative ∂G

∂Nstarts from 0 at the threshold and

rises with the current. Derivative ∂R

∂Nis initially approximately

constant (see figure A4) and its value is related to the recom-

bination itself. Because the recombination rate is a superlinear

function of carrier concentration, a higher value of the deriv-

ative means a higher value of the recombination rate. As a res-

ult, the current IHfor which the S21 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 con-

tain ωi, the value of IH(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 3one can

notice that our simulations probably underestimated the rate

of the non-radiative recombination, because in the measure-

ments the highest resonant peak appears at a current some-

where in the vicinity of the current of 2mA, while in the sim-

ulations it is slightly above 1 mA. The formula (A1) for the

resonant frequency can be radically simplified under the con-

dition ωi≫∂G

∂N≫∂R

∂N:

fr(I) = √2ωi(I)∂G

∂N(I)

2π.(A4)

Although this formula is an approximation and its validity is

limited by the conditions discussed above, it gives an idea

how frand hence f3 dB can be increased. Frequency f3 dB itself,

unlike fr, 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 fr, one should increase

ωiby reducing the number of DBR pairs (reducing the Q-

factor) and increase ∂G

∂Nby, 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 import-

ant factor, formula equation (A4) is no longer valid, and, as

presented in figures 3and 4,fris saturated. The other way

to increase ∂G

∂Ncould 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 equa-

tion (A2), with higher ωiwe get higher resonant peaks which

may have an adverse impact on large-signal modulation pat-

terns. The considerations presented above show that a proper

choice of the number of QWs can help to increase the mod-

ulation 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 adjus-

ted in such a way that all four structures emitted at the same

wavelength. In each case, the QWs were placed symmetric-

ally on both sides of the anti-node of the standing way, so as

to maximize the confinement factor. In figure A6 f3 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

J. Phys. D: Appl. Phys. 53 (2020) 345101 M Wasiak et al

Figure A6. f3 dB as a function of the driving current for VCSELs

with different number of QWs in the active region. The curve for 5

QW is the same with the one in figure 4.

same in all the cases. When there are more QWs, the gain per

a single well is lower, and as a result derivative ∂g

∂Nis bigger.

As discussed above, this increases the resonance frequency

and hence f3 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.

ORCID iDs

Michał Wasiak https://orcid.org/0000-0002-9569-4265

Patrycja ´

Spiewak https://orcid.org/0000-0002-2820-7420

Marcin Ge¸bski https://orcid.org/0000-0002-7307-1761

Emilia Pruszy´

nska-Karbownik https://orcid.org/0000-

0002-5973-9825

Paulina Komar https://orcid.org/0000-0003-1968-5591

James A Lott https://orcid.org/0000-0003-4094-499X

Robert P Sarzała https://orcid.org/0000-0002-2929-0844

References

[1] Haglund E, Westbergh P, Gustavsson J S, Haglund E P,

Larsson A, Geen M and Joel A 2015 30 GHz bandwidth

850 nm VCSEL with sub-100 fJ/bit energy dissipation at

25–50 Gbit/s Electron. Lett. 51 1096–8

[2] Larsson A et al 2018 1060 nm VCSELs for long-reach

optical interconnects Optical Fiber Technology

44 36–42

[3] Haghighi N, Moser P and Lott J A 2019 Bandwidth and

optical output power of VCSELs and VCSEL arrays Proc.

SPIE vol 10938

[4] Dave H, Liao P, Fryslie S T M, Gao Z, Thompson B J, Willner

A E and Choquette K D 2019 Digital modulation of

coherently-coupled 2×1 vertical-cavity surface-emitting

laser arrays IEEE Photon. Technol. Lett. 31 173–6

[5] Cheng C, Ledentsov N, Khan Z, Yen J, Ledentsov N N and Shi

J 2019 Ultrafast Zn-diffusion and oxide-relief 940 nm

vertical-cavity surface-emitting lasers under

high-temperature operation IEEE J. Sel. Top. Quantum

Electron. 25 1–7

[6] Kuchta D M, Rylyakov A V, Schow C L, Proesel J E, Baks C

W, Westbergh P, Gustavsson J S and Larsson A 2014 A 50

Gb/s NRZ modulated 850 nm VCSEL transmitter operating

error free to 90 ◦CJ. Lightwave Technol. 33 802–10

[7] Szczerba K, Lengyel T, Karlsson M, Andrekson P A and

Larsson A 2016 94-Gb/s 4-PAM using an 850-nm VCSEL,

pre-emphasis and receiver equalization IEEE Photon.

Technol. Lett. 28 2519–21

[8] Liu A, Wolf P, Lott J A and Bimberg D 2019 Vertical-cavity

surface-emitting lasers for data communication and sensing

Photon. Res. 7121–36

[9] Kuchta D M 2017 High capacity VCSEL-based links 2017

Optical Fiber Conf. and Exhibition (OFC) pp 1–94

[10] Coldren L A and Corzine S W 1995 Dynamics Effects Diode

Lasers and Photonic Integrated Circuits (New York: Wiley)

[11] Siu Y, Wong W N, Shum P and Li E H 1997 Theoretical

analysis of modulation response and second-order harmonic

distortion in vertical-cavity surface-emitting lasers J.

Quantum Electron. 32 2139–47

[12] Chow W W, Schneider H, Koch S W, Chang C-H,

Chrostowski L and Chang-Hasnain C 2002 Nonequilibrium

model for semiconductor laser modulation response J.

Quantum Electron. 38 402–9

[13] Tucker R S 1981 Large-signal circuit model for simulation of

injection-laser modulation dynamics IEE Proc. I

(Solid-State and Electron Devices) 128 180–4

[14] Arnold G, Russer P and Petermann K 1982 Modulation of

Laser Diodes (Berlin: Springer) pp 213–42

[15] Petermann K 2012 Intensity modulation characteristics of

laser diodes (Springer Science & Business Media) vol 3

pp 78–118

[16] Michalzik R 2013 Vcsel Fundamentals (Berlin: Springer)

pp 19–75

[17] Xu D et al 2009 Room-temperature continuous-wave

operation of the In(Ga)As/GaAs quantum-dot VCSELs for

the 1.3 µm optical-fibre communication Semicond. Sci.

Technol. 24 055003

[18] Sarzał A R, Marciniak M and Czyszanowski T 2018 Thermal

properties of GaN-based semiconductor-metal

subwavelength grating VCSELs and novel current injection

scheme J. Phys. D: Appl. Phys. 51 285102

[19] Sarzał a R P, ´

Spiewak P and Wasiak M 2019 Influence of

resonator length on performance of nitride TJ VCSEL IEEE

J. Quantum Electron. 55 1–9

[20] Wasiak M, ´

Spiewak P, Moser P, Walczak J, Sarzał a R P,

Czyszanowski T and Lott J A 2016 Numerical model of

capacitance in vertical-cavity surface-emitting lasers J.

Phys. D: Appl. Phys. 49 175104

[21] Tibaldi A, Bertazzi F, Goano M, Michalzik R and Debernardi

P 2019 VENUS: A vertical-cavity surface-emitting laser

electro-opto-thermal numerical simulator IEEE J. Sel.

Topics Quantum Electron. 25 1–12

[22] Kalosha V P, Shchukin V A, Ledentsov N and Ledentsov N N

2019 Comprehensive analysis of electric properties of

oxide-confined vertical-cavity surface-emitting lasers IEEE

J. Sel. Topics Quantum Electron. 25 1–9

[23] Wasiak M 2011 Mathematical rigorous approach to simulate

an over-threshold VCSEL operation Phys. E Low Dimens.

Syst. Nanostruct. 43 1439–44

[24] Dems M, Kotynski R and Panajotov K 2005 PlaneWave

Admittance Method—a novel approach for determining the

electromagnetic modes in photonic structures Opt. Express

13 3196–207

J. Phys. D: Appl. Phys. 53 (2020) 345101 M Wasiak et al

[25] Sarzał a R P, Mendla P, Wasiak M, Ma´

ckowiak P, Bugajski M

and Nakwaski W 2004 Comprehensive self-consistent

three-dimensional simulation of an operation of the

gaas-based oxide-confined 1.3-µm quantum-dot

(InGa)As/GaAs vertical-cavity surface-emitting lasers Opt.

Quantum Electron. 36 331–47

[26] Sarzał a R P 2005 Designing strategy to enhance mode

selectivity of higher-output oxide-confined

vertical-cavity surface-emitting lasers Appl. Phys. A

81 275–83

[27] Agrawal G P and Dutta N K 2013 Basic Concepts

Semiconductor Lasers (Dordrecht: Kluwer)

https://doi.org/10.1007/978-1-4613-0481-4

[28] Botez D 1978 Analytical approximation of the radiation

confinement factor for the TE0mode of a double

heterojunction laser IEEE J. Quantum Electron.

14 230–2

[29] Paoli T L and Ripper J E 1970 Direct modulation of

semiconductor lasers Proc. IEEE 58 1457–65

[30] Adams M J 1973 Rate equations and transient phenomena in

semiconductor lasers Opto-Electron. 5201–15

[31] Boers P M, Vlaardingerbroek M T and Danielsen M 1975

Dynamic behaviour of semiconductor lasers Electron. Lett.

11 206–8

[32] Haglund E 2015 VCSELs for High-Speed, Long-Reach, and

Wavelength-Multiplexed Optical Interconnects src PhD

thesis Chalmers University of Technology

[33] Liu Y, Ng W-C, Choquette K D and Hess K 2005 Numerical

investigation of self-heating effects of oxide-confined

vertical-cavity surface-emitting lasers IEEE J. Quantum

Electron. 41 15–25

[34] Giudice G E, Kuksenkov D V, Temkin H and Lear K L 1999

Differential carrier lifetime in oxide-confined vertical cavity

lasers obtained from electrical impedance measurements

Appl. Phys. Lett. 74 899–901

[35] Tsai C-Y, Tsai C-Y, Yu-Hwa L and Eastman L F 1996

Nonlinear gain coefficients in semiconductor lasers: effects

of carrier heating IEEE J. Quantum Electron. 32 201–12

[36] Westbergh P, Gustavsson J S, Kögel B, Haglund Å and

Larsson A 2011 Impact of photon lifetime on high-speed

VCSEL performance IEEE J. Sel. Topics Quantum

Electron. 17 1603–13