Stochastic Simulation of Emission Spectra and Classical Photon Statistics of Quantum Dot Superluminescent Diodes

We present a stochastic procedure to investigate the correlation spectra of quantum dot superluminescent diodes. The classical electric field of a diode is formed by a polychromatic superposition of many independent stochastic oscillators. Assuming fields with individual carrier frequencies, Lorentzian linewidths and amplitudes we can form any relevant experimental spectrum using a least square fit. This is illustrated for Gaussian and Lorentzian spectra, Voigt profiles and box shapes. Eventually, the procedure is applied to an experimental spectrum of a quantum dot superluminescent diode which determines the first- and second-order temporal correlation functions of the emission. We find good agreement with the experimental data and a quantized treatment. Thus, a stochastic field represents broadband light emitted by quantum dot superluminescent diodes.


Introduction
Modern day optical applications like optical coherence tomography [3,4] and ghost imaging [1,5,6] make use of the unique emission properties of spectrally broadband light-emitting quantum dot superluminescent diodes (QDSLD). By using specialized waveguide geometries and gain materials, QDSLDs are able to combine high output intensities, spatially directed emission and spectral widths in the THz regime. Hence, they fill the gap in the family of semiconductor-based optical emitters between coherent laser diodes and incoherent light emitting diodes. After being proposed in 1973 [7], research on the characteristics of QDSLDs has been intensified in recent years after Boitier et al. [8] enabled the direct measurement of coherence times in the femtosecond regime using two-photon absorption in semiconductors. This research includes experimental studies of emission and photon statistical characteristics of QDSLDs [9,10,11], as well as theoretical investigations based on rate equations [12,13,14], travelling wave approaches [15], finite element methods [16] and a quantized treatment [2,17]. To this day, effort is being put into developing more efficient and high-powered QDSLDs [18,19,20].
Adding a new perspective to the investigation of QDSLDs, we discuss a stochastic model for the emission in this article. Stochastic approaches have long proven to have a wide-ranging field of applications in biology [21,22,23], engineering [24], finance [25], quantum many-body physics [26], soft-matter physics [27,28], optics [29,30] and many more scientific areas. We develop a model to describe experimental emission spectra from a superposition of stochastic fields. Using least square fits, we determine Lorentzian linewidths, carrier frequencies and amplitudes of the individual fields to model experimental spectra. This is illustrated for Gaussian-, Lorentzian-, Voigt-as well as bandpass spectra and applied to the experimental spectrum of a QDSLD [1]. Using numerical simulations, we determine first-and second-order temporal correlation functions of the electric field and calculate the spectral power density of the resulting emission.
The article is organized as follows: the stochastic model of emission spectra is developed in Sec. 2. It consists of individual classical fields, which are described by a distinct stochastic differential equation. After investigating properties of these fields, relevant spectra are modelled as a superposition. This is applied to a specific experimental spectrum produced by a QDSLD [1]. The model is subsequently used to calculate the emission spectrum of the diode in Sec. 3 and the normalized stationary second-order temporal correlation function in Sec. 4. A conclusion is given in Sec. 5. An appendix summarizes the convergence properties of the simulation schemes.

Stochastic model of emission spectra
The classical electric field of a diode results from a superposition of stochastic fields. Hence, the electric field outside of the diode reads with the number of fields N and ε j (t) the j-th complex field amplitude.

Ornstein-Uhlenbeck process
An individual classical field ε(t) ∈ C is modelled as a complex Ornstein-Uhlenbeck process [31,32]. This is described by the Ito stochastic differential equation with the carrier frequency ν 0 , the linewidth γ, the diffusion constant D = √ γI, the mean intensity of the electric field I = lim t→∞ |ε(t)| 2 and the complex Wiener noise increment dW (t) ∈ C, whose properties are given by dW (t) = 0 and |dW (t)| 2 = dt.
The stationary first-order temporal correlation function reads [32,33] The spectral power density is given by the Fourier transform of eq. (3) in accordance to the Wiener-Khintchine theorem [34,35]. This yields Furthermore, the stationary normalized second-order temporal correlation function is given by the Siegert relation [33,36]

Stochastic simulation
In addition to analytical results we perform numerical simulations of eq. (2). In order to obtain an efficient simulation procedure we separate the rapid oscillating carrier frequency by the transformation As the diffusion constant D is independent of the electric field amplitude η(t), the Euler scheme [37] can be used to achieve strong convergence of order 1.0 (see Appendix A). Therefore the electric field amplitude can be simulated iteratively with the discrete time step ∆t = t i+1 − t i and ∆W a complex Gaussian random process with mean ∆W = 0 and variance |∆W | 2 = ∆t. The first-order temporal correlation function of eq. (3) is calculated from a sample average over M realizations long after the transient regime t s 1/γ j . ε (m) (t) is the m-th realization of the electric field. This result can be used to calculate the spectral power density of the emission using a Fourier transformation.
The determination of the normalized second-order temporal correlation function (6) can be split into two separate calculations. The first-order temporal correlation functions in the denominator can be simulated according to eq. (9), while the second-order correlation function in the numerator can be calculated as Simulation results as well as analytical calculations of the first-and second-order temporal correlation properties of an individual field can be seen in Fig. 1. The simulations show good agreement with the analytical results for the given parameters (γ = 0.5 THz, I = 1, ν 0 = 10 THz, ∆t = 0.01 ps, M = 10 4 ).  s (τ ) versus time τ for the electric field described by eq. (2). Simulation results (black, solid) and analytical expressions (yellow, dashed) were calculated with I = 1 and ν 0 = 10 THz. The correlation functions were determined using ∆t = 0.01 ps and M = 10 4 realizations.

Modelling of emission shapes
The emission of a diode (1) is described as the superposition of N independent classical fields with individual linewidths γ j , mean intensities I j and central frequencies ν j . The stationary first-order temporal correlation function reads Thus, the spectral power density is the incoherent sum of the individual spectra This model can be used to approximate a wide range of shapes through the adjustment of the 3N free parameters γ j , I j and ν j in eq. (12) by means of a least square fit, minimizing the error functional for a test spectrum S t (ν) at discrete frequencies ν i . Examples of interest [38] are given by Gaussian spectra Lorentzian spectra Voigt profiles with the complementary error function erfc(z), and bandwidth limited box shapes This is illustrated in Fig. 2. Please note that all spectra are normalized to

Model of quantum dot superluminescent diode emission
Superluminescent diodes are semiconductor-based light sources, which are characterized by spatially directed emission and spectral widths in the THz regime. The experiments with QDSLDs [1] had an active medium consisting of inhomogeneously broadened InAs/InGaAs quantum dot layers. The optical power spectrum has a Gaussian shape (see Fig. 3). The emission of the diode is modelled by N = 30 individual oscillators using a least square fit (see eq. (13)) to an experimental spectrum S e (ν) [1] to determine linewidths γ j , mean intensities I j and central frequencies ν j describing the emission. This is illustrated in Fig. 3.

QDSLD emission spectrum
The description of the QDSLD is used to simulate the output spectrum of such a diode. For this, the central frequencies ν j , linewidths γ j and mean intensities I j determined in Sec. 2.4 are used to simulate the individual electric fields ε j (t) according to eq. (8). The electric field produced by the diode ε d (t) results according to eq. (1). Subsequently, the stationary first-order temporal correlation function G d (τ ) is calculated according to eq. (9) using M = 10 4 realizations of the diode field ε d (t). The spectral power density of the emission S d (ν) is determined using a Fourier transformation.  The result of the simulation (see Fig. 4) shows good agreement with the experimental optical power spectrum [1]. We define the width of the spectral power density [40,41] This yields b d = 4.51 THz, implying a coherence time of τ c,d = 1/b d = 221.9 fs, which matches the experimental results of b e = 4.29 THz and τ c,e = 233 fs very well. The developed formalism of modelling the emission of a quantum dot superluminescent diode as a superposition of individual oscillators is therefore suitable to describe the spectral power density of the diode.

Second-order temporal correlation function
In addition to the investigation of the optical power spectrum, the developed formalism can be used to investigate the classical photon statistics of the QDSLD emission. For this, M realizations of the electric field ε d (t) determined in sec. 3 are used to calculate the stationary normalized second-order temporal correlation function g (2) d (τ ) of the emission (6,9,10). The result for the central frequencies ν j , linewidths γ j and mean intensities I j determined in section 2.4 is illustrated in Fig. 5. There is good agreement between the simulation and the experimental data g (2) e (τ ). The central degree of second-order temporal coherence g (2) d (τ = 0) 2 indicates a Gaussian photon distribution, which was also shown experimentally. Therefore, the developed formalism is suited for classical photon statistical investigations of QDSLDs.

Conclusion
In this article we study a stochastic model to describe experimental emission spectra. These are considered to result as a superposition of individual complex Ornstein-Uhlenbeck processes. The firstand second-order temporal correlation properties of these oscillators are investigated analytically and numerically. We can approximate Gaussian-, Lorentzian-, Voigt-and bandwidth limited spectra by determination of Lorentzian linewidths, carrier frequencies and amplitudes of the individual oscillators using least square fits.
The developed procedure is applied to the emission properties of quantum dot superluminescent diodes. Simulation parameters are extracted from a least square fit to an experimental spectrum [1]. These are used to simulate the QDSLD emission and calculate first-and second-order temporal correlation properties. The determined spectral power density of the emission, resulting from a Fourier transformation of the stationary first-order temporal correlation function, shows good agreement with the experimental results regarding the shape of the spectral line, as well as spectral width and coherence time. Additionally, calculating the stationary normalized second-order temporal correlation function results in a a central degree of second-order temporal coherence g (2) (τ = 0) 2. This indicates a Gaussian photon distribution, which is in agreement with experiments and former theoretical investigations.
The stochastic description of QDSLD emission offers a straightforward perspective on the process of light generation inside QDSLDs, describing it as a superposition of individual classical oscillators. More data on the emission characteristics of the constituents of QDSLDs can lead to a better understanding and contribute to the design of new diodes. Furthermore, this approach can be used in the investigation of other phenomena shown by QDLSDs such as temperature dependent intensity noise reduction [11].

Acknowledgement
We thank Sébastien Blumenstein for the provision of experimental data and Prof. Wolfgang Elsäßer for stimulating discussions.

A Convergence of stochastic simulations
Consider the Ito stochastic differential equation [32] with the drift term a(x) and the diffusion term b(x). Identifying x(t 0 ) = x 0 , the formal solution of this equation is given by integration: The goal of time discrete maps x(t i+1 ) = F (x i ) of stochastic differential equations is the approximation of a solution x(t) up to a order of convergence γ. Such a scheme is said to converge strongly with order γ > 0, if for the final time instant T and N = T /∆ there is a finite and ∆ 0 > 0 such that [37] for any time discretization 0 < ∆ < ∆ 0 . A strong Taylor scheme of order γ can be constructed by considering the Ito-Taylor expansion, which is obtained by continuously applying the integral form of Ito's formula [32] f ( with L 0 = a(x(t ))∂ x + (1/2)b 2 (x(t ))∂ 2 x and L 1 = b(x(t ))∂ x , to nonconstant terms inside the integrals of the formal solution (21). A criterium [37] for the terms of the Ito-Taylor expansion required for the associated strong Taylor scheme to achieve a desired order of strong convergence γ states, that a simulation scheme strongly converges to the order of an integer γ if it includes all combinations of integrals up to this order, with time differentials dt being of order 1 and Wiener noise increments dW (t) being of order 1/2. Simulation schemes of half-integer order γ additionally require the inclusion of the pure time integral of order γ + 1/2.
A strong convergence scheme of order 1/2 is the Euler scheme Discretizing the time steps, discrete map can be developed which yields x(t i+1 ) = x(t i ) + a(x(t i )) ∆t + b(x(t i )) ∆W, where ∆W is a Gaussian random process with ∆W = 0 and ∆W 2 = ∆t. To expand the Euler scheme to order 1.0 of strong convergence, the double stochastic integral appearing in the remainder R in eq. (24) has to be included, yielding This is called the Milstein scheme [42]. The double stochastic integral can be calculated [32] t t0 t t0 dW (t )dW (t ) = 1 2 which leads to the iteration rule for the Milstein method With increasing order of convergence γ, the simulation schemes become more complex and include an increasing number of stochastic increments.