Modern-day optical applications like optical coherence tomography   and ghost imaging    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 , research on the characteristics of QDSLDs has been intensified in recent years after Boitier et al.  enabled the direct measurement of coherence times in the femtosecond regime using two-photon absorption in semiconductors. This research is both theoretical and experimental nature. From the theoretical side, the understanding of light generation processes inside the diode is a main focus. For this, a plethora of approaches is being utilized including rate equation models   , travelling wave approaches , finite element methods  and quantized treatments  . Experimental studies focus on the determination of first- and second-order temporal correlation properties of QDSLDs  . In 2011, Blazek et al. were able to observe a temperature dependent suppression of intensity fluctuations using a broadband emitting QDSLD . To this day, effort is being put into developing more efficient and high-powered QDSLDs   .
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   , engineering , finance , quantum many-body physics  , soft-matter physics  , optics   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 . 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 Section 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 . The model is subsequently used to calculate the emission spectrum of the diode in Section 3 and the normalized stationary second-order temporal correlation function in Section 4. A conclusion is given in Section 5. An Appendix summarizes the convergence properties of the simulation schemes.
2. 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 the j-th complex field amplitude.
2.1. Ornstein-Uhlenbeck Process
An individual classical field is modelled as a complex Ornstein-Uhlenbeck process . This is described by the Ito stochastic differential equation 
with the carrier frequency , the linewidth , the diffusion constant , the mean intensity of the electric field and the complex Wiener noise increment , whose properties are given by and .
The stationary first-order temporal correlation function reads  
The spectral power density is given by the Fourier transform
of (3) in accordance to the Wiener-Khintchine theorem  . This yields
Furthermore, the stationary normalized second-order temporal correlation function is given by the Siegert relation  
2.2. Stochastic Simulation
In addition to analytical results, we perform numerical simulations of (2). In order to obtain an efficient simulation procedure, we separate the rapid oscillating carrier frequency by the transformation yielding
As the diffusion constant D is independent of the electric field amplitude , the Euler scheme  can be used to achieve strong convergence of order 1.0 (see Appendix). Therefore the electric field amplitude can be simulated iteratively
with the discrete time step and a complex Gaussian random process with mean and variance .
The first-order temporal correlation function of (3) is calculated from a sample average over M realizations
long after the transient regime . 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 (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 Figure 1.
Figure 1. (a) Stationary first-order temporal correlation function versus time , (b) spectral power density versus frequency with width (FWHM) THz and (c) stationary normalized second-order temporal correlation function versus time for the electric field described by (2). Simulation results (black, solid) and analytical expressions (yellow, dashed) were calculated with and . The correlation functions were determined using and realizations.
The simulations show good agreement with the analytical results for the given parameters ( , , , , ).
2.3. Modelling of Emission Shapes
The emission of a diode (1) is described as the superposition of N independent classical fields with individual linewidths , mean intensities and central frequencies . 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 3 N free parameters , and in (12) by means of a least square fit, minimizing the error functional
for a test spectrum at discrete frequencies . Examples of interest are given by Gaussian spectra 
Lorentzian spectra 
Voigt profiles 
with the complementary error function , and bandwidth limited box shapes
This is illustrated in Figure 2. Please note that all spectra are normalized to
2.4. 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  had an active medium consisting of inhomogeneously broadened InAs/InGaAs quantum dot layers. The optical power spectrum has a Gaussian shape (see Figure 3).
Figure 2. (a) Gaussian spectrum , (b) Lorentzian spectrum , (c) Voigt profile and (d) bandwidth limited box shape versus frequency (black, solid; , , ) modelled according to (12) with elementary oscillators (yellow, dashed). The oscillations at the edge of the box are a typical Gibbs phenomenon .
Figure 3. Experimental optical power spectrum (black, solid) versus frequency  with Gaussian fit (red, dotted), stochastic emission for oscillators according to (12) (yellow, dashed) and spectra of individual oscillators (solid, blue). The results of the Gaussian fit are a central frequency and a width .
Therefore the developed formalism is used to describe the emission of the diode, which is modelled by individual oscillators. Using a least square fit (see (13)) to an experimental spectrum  the linewidths , mean intensities and central frequencies describing the emission are determined. This is illustrated in Figure 3.
3. QDSLD Emission Spectrum
The optical power spectrum emitted by the QDSLD is simulated numerically. For this, the central frequencies , linewidths and mean intensities describing the QDSLD emission determined in Sec. 2.4 are used to calculate the individual electric fields according to (8). The electric field emitted by the diode results as a superposition of the individual field according to (1). Subsequently, the stationary first-order temporal correlation function is calculated according to (9) using realizations of the diode field . The spectral power density of the emission is determined using a Fourier transformation.
The result of the simulation (see Figure 4) shows good agreement with the experimental optical power spectrum . We define the width of the spectral power density  
This yields , implying a coherence time of , which matches the experimental results of and very well.
The method of modelling emission spectra as a superposition of individual oscillators is therefore suitable to describe the first-order temporal correlation properties of QDSLDs.
Figure 4. Experimental optical power spectrum (black, solid) versus frequency  with simulation results (yellow, dashed) for oscillators resulting from a Fourier transformation of the stationary first-order temporal correlation function calculated according to (9) with .
By extracting appropriate simulation parameters from an experimental spectrum, the electric field emitted by the diode can be simulated numerically and can be used to calculate the stationary first-order temporal correlation function and optical power spectrum of the emission.
4. 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, the electric field emitted by the diode already calculated in Sec. 3 can be reused. Instead of calculating first-order temporal correlation properties of the field, M realizations of are used to calculate the stationary normalized second-order temporal correlation function of the emission according to (6, 9, 10).
The result for the central frequencies , linewidths and mean intensities determined in Section 2.4 is illustrated in Figure 5. There is good agreement between the simulation and the experimental data . The central degree of second-order temporal coherence indicates a Gaussian photon distribution, which was also shown experimentally. Therefore, the developed formalism is suited for classical photon statistical investigations of QDSLDs.
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 first- and 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.
Figure 5. Experimental stationary normalized second-order temporal correlation function (black, solid) versus time  with simulation results (yellow, dashed) for oscillators calculated according to (9, 10) with .
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 . 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 central degree of second-order temporal coherence . 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 utilized in the investigation of other properties of QDSLDs. As it explains the statistical properties of the electric field emitted by the diode, it can be used in a classical explanation of temperature dependent intensity fluctuation suppression observed by Blazek et al.
We thank Sébastien Blumenstein for the provision of experimental data and Prof. Wolfgang Elsäßer for stimulating discussions.
Appendix. Convergence of Stochastic Simulations
Consider the Ito stochastic differential equation 
with the drift term and the diffusion term . Identifying , the formal solution of this equation is given by integration:
The goal of time discrete maps of stochastic differential equations is the approximation of a solution up to a order of convergence . Such a scheme is said to converge strongly with order , if for the final time instant T and there is a finite and such that 
for any time discretization . 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 
with and , to nonconstant terms inside the integrals of the formal solution (21). A criterium  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 being of order 1 and Wiener noise increments being of order 1/2. Simulation schemes of half-integer order additionally require the inclusion of the pure time integral of order .
A strong convergence scheme of order 1/2 is the Euler scheme
Discretizing the time steps, discrete map can be developed which yields
where is a Gaussian random process with and . To expand the Euler scheme to order 1.0 of strong convergence, the double stochastic integral appearing in the remainder in (24) has to be included, yielding
This is called the Milstein scheme . The double stochastic integral can be calculated 
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.
 Huang, D., Swanson, E.A., Lin, C.P., Schuman, J.S., Stinson, W.G., Chang, W., Hee, M.R., Flotire, T., Gregory, K., Puliafito, C.A. and Fujimoto, J.G. (1991) Science, 22, 1178-1181.
 Friedrich, F., Elsäßer, W. and Walser, R. (2020) Optics Communications, 458, Article ID: 124449.
 Kiethe, J., Heuer, A. and Jechow, A. (2017) Laser Physics Letters, 14, Article ID: 086201.
 Ozaki, N., Yamauchi, S., Hayashi, Y., Watanabe, E., Ohsato, H., Ikeda, N., Sugimoto, Y., Furuki, K., Oikawa, Y., Miyaji, K., Childs, D.T.D. and Hogg, R.A. (2019) Journal of Physics D: Applied Physics, 52, Article ID: 225105.
 Olver, F.W.J., Olde Daalhuis, A.B., Lozier, D.W., Schneider, B.I., Boisvert, R.F., Clark, C.W., Miller, B.R., Saunders, B.V., Cohl, H.S. and McClain, M.A. (2020) NIST Digital Library of Mathematical Functinos, Release 1.0.28 of 2020-09-15.