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

Kai Niklas Hansmann, Reinhold Walser

Technische Universität Darmstadt, Institut für Angewandte Physik, Hochschulstraße 4a, Darmstadt, Germany.

DOI: 10.4236/jmp.2021.121003                

Abstract

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 superposition of independent stochastic oscillators represents the first- and second-order correlation properties of broadband light emitted by quantum dot superluminescent diodes.

Keywords

Stochastic Simulation, Quantum Dot, Superluminescent Diode

1. Introduction

Modern-day optical applications like optical coherence tomography [1] [2] and ghost imaging [3] [4] [5] 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 [6], research on the characteristics of QDSLDs has been intensified in recent years after Boitier et al. [7] 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 [8] [9] [10], travelling wave approaches [11], finite element methods [12] and quantized treatments [13] [14]. Experimental studies focus on the determination of first- and second-order temporal correlation properties of QDSLDs [15] [16]. In 2011, Blazek et al. were able to observe a temperature dependent suppression of intensity fluctuations $g^{\left(2\right)}\left(\tau =0\right)<2$ using a broadband emitting QDSLD [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] [27], soft-matter physics [28] [29], optics [30] [31] 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 [3]. 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 [3]. 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

${\epsilon }_{\text{d}}\left(t\right)={\displaystyle \underset{j=1}{\overset{N}{\sum }}}{\epsilon }_j\left(t\right)\mathrm{,}$ (1)

with the number of fields N and ${\epsilon }_j\left(t\right)$ the j-th complex field amplitude.

2.1. Ornstein-Uhlenbeck Process

An individual classical field $\epsilon \left(t\right)\in ℂ$ is modelled as a complex Ornstein-Uhlenbeck process [32]. This is described by the Ito stochastic differential equation [33]

$\text{d}\epsilon \left(t\right)=\left(\text{i}{\nu }_0-\gamma \right)\epsilon \left(t\right)\text{d}t+D\text{d}W\left(t\right)\mathrm{,}$ (2)

with the carrier frequency ${\nu }_0$, the linewidth $\gamma$, the diffusion constant $D=\sqrt{\gamma I}$, the mean intensity of the electric field $I={\lim }_{t\to \infty }\langle {\left|\epsilon \left(t\right)\right|}^2\rangle$ and the complex Wiener noise increment $\text{d}W\left(t\right)\in ℂ$, whose properties are given by $\langle \text{d}W\left(t\right)\rangle =0$ and $\langle {\left|\text{d}W\left(t\right)\right|}^2\rangle =\text{d}t$.

The stationary first-order temporal correlation function reads [33] [34]

$G_{\text{s}}^{\left(1\right)}\left(\tau \right)=\underset{t\to \infty }{lim}\langle {\epsilon }^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\epsilon \left(t+\tau \right)\rangle =I{\text{e}}^{-\gamma \left|\tau \right|-\text{i}{\nu }_0\tau }\mathrm{.}$ (3)

The spectral power density is given by the Fourier transform

$S\left(\nu \right)=\frac{1}{\sqrt{2\pi }}{\displaystyle {\int }_{-\infty }^{\infty }}\text{d}\tau G_{\text{s}}^{\left(1\right)}\left(\tau \right){\text{e}}^{\text{i}\nu \tau }$ (4)

of (3) in accordance to the Wiener-Khintchine theorem [35] [36]. This yields

$S\left(\nu \right)=I\sqrt{\frac{2}{\pi }}\frac{\gamma }{{\left(\nu -{\nu }_0\right)}^2+{\gamma }^2}\mathrm{,}\frac{1}{\sqrt{2\pi }}{\displaystyle {\int }_{-\infty }^{\infty }}\text{d}\nu S\left(\nu \right)=I\mathrm{.}$ (5)

Furthermore, the stationary normalized second-order temporal correlation function is given by the Siegert relation [34] [37]

$g_{\text{s}}^{\left(2\right)}\left(\tau \right)=\underset{t\to \infty }{lim}\frac{\langle {\epsilon }^{\mathrm{<mo>\; *\; </mo>}}\left(t\right){\epsilon }^{\mathrm{<mo>\; *\; </mo>}}\left(t+\tau \right)\epsilon \left(t+\tau \right)\epsilon \left(t\right)\rangle }{\langle {\epsilon }^{\mathrm{<mo>\; *\; </mo>}}\left(t\right)\epsilon \left(t\right)\rangle \langle {\epsilon }^{\mathrm{<mo>\; *\; </mo>}}\left(t+\tau \right)\epsilon \left(t+\tau \right)\rangle }=1+{\text{e}}^{-2\gamma \left|\tau \right|}\mathrm{.}$ (6)

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 $\epsilon \left(t\right)=\eta \left(t\right){\text{e}}^{-\text{i}{\nu }_{0}t}$ yielding

$\text{d}\eta \left(t\right)=-\gamma \eta \left(t\right)\text{d}t+D\text{d}W\left(t\right)\mathrm{.}$ (7)

As the diffusion constant D is independent of the electric field amplitude $\eta \left(t\right)$, the Euler scheme [38] can be used to achieve strong convergence of order 1.0 (see Appendix). Therefore the electric field amplitude can be simulated iteratively

$\eta \left(t_{i+1}\right)=\eta \left(t_i\right)-\gamma \eta \left(t_i\right)\Delta t+D\Delta W\mathrm{,}$ (8)

with the discrete time step $\Delta t=t_{i+1}-t_i$ and $\Delta W$ a complex Gaussian random process with mean $\langle \Delta W\rangle =0$ and variance $\langle {\left|\Delta W\right|}^2\rangle =\Delta t$.

The first-order temporal correlation function of (3) is calculated from a sample average over M realizations

$G^{\left(1\right)}\left(\tau \right)=\frac{1}{M}{\displaystyle \underset{m=1}{\overset{M}{\sum }}}{\left({\epsilon }^{\left(m\right)}\left(t_s\right)\right)}^{\mathrm{<mo>\; *\; </mo>}}{\epsilon }^{\left(m\right)}\left(t_s+\tau \right)$ (9)

long after the transient regime $t_s\gg 1/{\gamma }_j$. ${\epsilon }^{\left(m\right)}\left(t\right)$ 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

$G^{\left(2\right)}\left(\tau \right)=\frac{1}{M}{\displaystyle \underset{m=1}{\overset{M}{\sum }}}{\left|{\epsilon }^{\left(m\right)}\left(t_s+\tau \right){\epsilon }^{\left(m\right)}\left(t_s\right)\right|}^2\mathrm{.}$ (10)

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 $G_{\text{s}}^{\left(1\right)}\left(\tau \right)$ versus time $\tau$, (b) spectral power density $S\left(\nu \right)$ versus frequency $\nu$ with width (FWHM) $2\gamma =1$ THz and (c) stationary normalized second-order temporal correlation function $g_{\text{s}}^{\left(2\right)}\left(\tau \right)$ versus time $\tau$ for the electric field described by (2). Simulation results (black, solid) and analytical expressions (yellow, dashed) were calculated with $I=1$ and ${\nu }_0=10\text{THz}$. The correlation functions were determined using $\Delta t=0.01\text{ps}$ and $M={10}^4$ realizations.

The simulations show good agreement with the analytical results for the given parameters ( $\gamma =0.5\text{THz}$, $I=1$, ${\nu }_0=10\text{THz}$, $\Delta t=0.01\text{ps}$, $M={10}^4$ ).

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 ${\gamma }_j$, mean intensities $I_j$ and central frequencies ${\nu }_j$. The stationary first-order temporal correlation function reads

$G_{\text{d}}^{\left(1\right)}\left(\tau \right)=\underset{t\to \infty }{lim}\langle {\epsilon }_{\text{d}}^{\mathrm{<mo>\; *\; </mo>}}\left(t\right){\epsilon }_{\text{d}}\left(t+\tau \right)\rangle =\underset{t\to \infty }{lim}{\displaystyle \underset{j=1}{\overset{N}{\sum }}}\langle {\epsilon }_j^{\mathrm{<mo>\; *\; </mo>}}\left(t\right){\epsilon }_j\left(t+\tau \right)\rangle \mathrm{.}$ (11)

Thus, the spectral power density is the incoherent sum of the individual spectra

$S_{\text{d}}\left(\nu \right)={\displaystyle \underset{j=1}{\overset{N}{\sum }}}{S}_j\left(\nu \right)\mathrm{.}$ (12)

This model can be used to approximate a wide range of shapes through the adjustment of the 3 N free parameters ${\gamma }_j$, $I_j$ and ${\nu }_j$ in (12) by means of a least square fit, minimizing the error functional

$e={\displaystyle \underset{i}{\sum }}{\left(S_{\text{t}}\left({\nu }_i\right)-S_{\text{d}}\left({\nu }_i\right)\right)}^2$ (13)

for a test spectrum $S_{\text{t}}\left(\nu \right)$ at discrete frequencies ${\nu }_i$. Examples of interest are given by Gaussian spectra [39]

$S_{\text{g}}\left(\nu \right)=\frac{1}{\sqrt{{\sigma }^2}}{\text{e}}^{-\frac{{\left(\nu -{\nu }_0\right)}^2}{2{\sigma }^2}}\mathrm{,}$ (14)

Lorentzian spectra [39]

$S_{\text{l}}\left(\nu \right)=\sqrt{\frac{2}{\pi }}\frac{\gamma }{{\left(\nu -{\nu }_0\right)}^2+{\gamma }^2}\mathrm{,}$ (15)

Voigt profiles [39]

$S_{\text{v}}\left(\nu \right)=\frac{1}{\sqrt{{\sigma }^2}}\text{Re}\left\{{\text{e}}^{-z^2}\text{erfc}\left(-\text{i}z\right)\right\}\mathrm{,}z=\frac{\nu +\text{i}\gamma }{\sigma \sqrt{2}}\mathrm{,}$ (16)

with the complementary error function $\text{erfc}\left(z\right)$, and bandwidth limited box shapes

$S_{\text{b}}\left(\nu \right)=\{\begin{array}{ll}\sqrt{2\pi }/\gamma \mathrm{,}\hfill & \text{for}\left|\nu \right|\le \gamma /2\mathrm{,}\hfill \\ \mathrm{0,}\hfill & \text{else}\text{.}\hfill \end{array}$ (17)

This is illustrated in Figure 2. Please note that all spectra are normalized to

$\frac{1}{\sqrt{2\pi }}{\displaystyle {\int }_{-\infty }^{\infty }}\text{d}\nu S\left(\nu \right)=1.$ (18)

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 [3] 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 $S_{\text{g}}\left(\nu \right)$, (b) Lorentzian spectrum $S_{\text{l}}\left(\nu \right)$, (c) Voigt profile $S_{\text{v}}\left(\nu \right)$ and (d) bandwidth limited box shape $S_{\text{b}}\left(\nu \right)$ versus frequency $\nu$ (black, solid; ${\nu }_0=0$, $\sigma =1$, $\gamma =1$ ) modelled according to (12) with $N=30$ elementary oscillators (yellow, dashed). The oscillations at the edge of the box are a typical Gibbs phenomenon [40].

Figure 3. Experimental optical power spectrum $S_{\text{e}}\left(\nu \right)$ (black, solid) versus frequency $\nu$ [3] with Gaussian fit $S_{\text{g}}\left(\nu \right)$ (red, dotted), stochastic emission $S_{\text{d}}\left(\nu \right)$ for $N=30$ oscillators according to (12) (yellow, dashed) and spectra of individual oscillators $S_j\left(\nu \right)$ (solid, blue). The results of the Gaussian fit are a central frequency ${\nu }_0=239.9\text{THz}$ and a width $\sigma =1.17\text{THz}$.

Therefore the developed formalism is used to describe the emission of the diode, which is modelled by $N=30$ individual oscillators. Using a least square fit (see (13)) to an experimental spectrum $S_{\text{e}}\left(\nu \right)$ [3] the linewidths ${\gamma }_j$, mean intensities $I_j$ and central frequencies ${\nu }_j$ 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 ${\nu }_j$, linewidths ${\gamma }_j$ and mean intensities $I_j$ describing the QDSLD emission determined in Sec. 2.4 are used to calculate the individual electric fields ${\epsilon }_j\left(t\right)$ according to (8). The electric field emitted by the diode ${\epsilon }_{\text{d}}\left(t\right)$ results as a superposition of the individual field ${\epsilon }_j\left(t\right)$ according to (1). Subsequently, the stationary first-order temporal correlation function $G_{\text{d}}^{\left(1\right)}\left(\tau \right)$ is calculated according to (9) using $M={10}^4$ realizations of the diode field ${\epsilon }_{\text{d}}\left(t\right)$. The spectral power density of the emission $S_{\text{d}}\left(\nu \right)$ is determined using a Fourier transformation.

The result of the simulation (see Figure 4) shows good agreement with the experimental optical power spectrum [3]. We define the width of the spectral power density [41] [42]

$b=\frac{1}{{\displaystyle {\int }_{-\infty }^{\infty }}\text{d}\nu S^2\left(\nu \right)}\mathrm{.}$ (19)

This yields $b_{\text{d}}=4.51\text{THz}$, implying a coherence time of ${\tau }_{\text{c,d}}=1/b_{\text{d}}=221.9\text{fs}$, which matches the experimental results of $b_{\text{e}}=4.29\text{THz}$ and ${\tau }_{\text{c,e}}=233\text{fs}$ 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 $S_{\text{e}}\left(\nu \right)$ (black, solid) versus frequency $\nu$ [3] with simulation results $S_{\text{d}}\left(\nu \right)$ (yellow, dashed) for $N=30$ oscillators resulting from a Fourier transformation of the stationary first-order temporal correlation function $G_{\text{d}}^{\left(1\right)}\left(\tau \right)$ calculated according to (9) with $M={10}^4$.

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 ${\epsilon }_{\text{d}}\left(t\right)$ 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 ${\epsilon }_{\text{d}}\left(t\right)$ are used to calculate the stationary normalized second-order temporal correlation function $g_{\text{d}}^{\left(2\right)}\left(\tau \right)$ of the emission according to (6, 9, 10).

The result for the central frequencies ${\nu }_j$, linewidths ${\gamma }_j$ and mean intensities $I_j$ determined in Section 2.4 is illustrated in Figure 5. There is good agreement between the simulation and the experimental data $g_{\text{e}}^{\left(2\right)}\left(\tau \right)$. The central degree of second-order temporal coherence $g_{\text{d}}^{\left(2\right)}\left(\tau =0\right)\simeq 2$ indicates a Gaussian photon distribution, which was also shown experimentally. Therefore, the developed formalism is suited for classical photon statistical investigations of QDSLDs.

5. Conclusions

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 $g_{\text{e}}^{\left(2\right)}\left(\tau \right)$ (black, solid) versus time $\tau$ [3] with simulation results $g_{\text{d}}^{\left(2\right)}\left(\tau \right)$ (yellow, dashed) for $N=30$ oscillators calculated according to (9, 10) with $M={10}^4$.

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 [3]. 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 $g^{\left(2\right)}\left(\tau =0\right)\simeq 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 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.

Acknowledgements

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 [33]

$\text{d}x\left(t\right)=a\left(x\left(t\right)\right)\text{d}t+b\left(x\left(t\right)\right)\text{d}W\left(t\right)\mathrm{,}$ (20)

with the drift term $a\left(x\right)$ and the diffusion term $b\left(x\right)$. Identifying $x\left(t_0\right)=x_0$, the formal solution of this equation is given by integration:

$x\left(t\right)=x_0+{\displaystyle {\int }_{t_0}^t}\text{d}{t}^{\prime }a\left(x\left(t^{\prime }\right)\right)+{\displaystyle {\int }_{t_0}^t}\text{d}W\left(t^{\prime }\right)b\left(x\left(t^{\prime }\right)\right)$ (21)

The goal of time discrete maps $x\left(t_{i+1}\right)=F\left(x_i\right)$ of stochastic differential equations is the approximation of a solution $x\left(t\right)$ up to a order of convergence $\gamma$. Such a scheme is said to converge strongly with order $\gamma >0$, if for the final time instant T and $N=T/\Delta$ there is a finite $\epsilon$ and ${\Delta }_0>0$ such that [38]

$\langle \left|x\left(T\right)-x\left(t_N\right)\right|\rangle \le \epsilon {\Delta }^{\gamma }$ (22)

for any time discretization $0<\Delta <{\Delta }_0$. A strong Taylor scheme of order $\gamma$ can be constructed by considering the Ito-Taylor expansion, which is obtained by continuously applying the integral form of Ito’s formula [33]

$f\left(x\left(t\right)\right)=f\left(x_0\right)+{\displaystyle {\int }_{t_0}^t}\text{d}{t}^{\prime }{L}^{0}f\left(x\left(t^{\prime }\right)\right)+{\displaystyle {\int }_{t_0}^t}\text{d}W\left(t^{\prime }\right)L^{1}f\left(x\left(t^{\prime }\right)\right)\mathrm{,}$ (23)

with $L^0=a\left(x\left(t^{\prime }\right)\right){\partial }_x+\left(1/2\right)b^2\left(x\left(t^{\prime }\right)\right){\partial }_x^2$ and $L^1=b\left(x\left(t^{\prime }\right)\right){\partial }_x$, to nonconstant terms inside the integrals of the formal solution (21). A criterium [38] for the terms of the Ito-Taylor expansion required for the associated strong Taylor scheme to achieve a desired order of strong convergence $\gamma$ states, that a simulation scheme strongly converges to the order of an integer $\gamma$ if it includes all combinations of integrals up to this order, with time differentials $\text{d}t$ being of order 1 and Wiener noise increments $\text{d}W\left(t\right)$ being of order 1/2. Simulation schemes of half-integer order $\gamma$ additionally require the inclusion of the pure time integral of order $\gamma +1/2$.

A strong convergence scheme of order 1/2 is the Euler scheme

$x\left(t\right)=x_0+a\left(x_0\right){\displaystyle {\int }_{t_0}^t}\text{d}{t}^{\prime }+b\left(x_0\right){\displaystyle {\int }_{t_0}^t}\text{d}W\left(t^{\prime }\right)+\mathcal{R}\mathrm{.}$ (24)

Discretizing the time steps, discrete map can be developed which yields

$x\left(t_{i+1}\right)=x\left(t_i\right)+a\left(x\left(t_i\right)\right)\Delta t+b\left(x\left(t_i\right)\right)\Delta W\mathrm{,}$ (25)

where $\Delta W$ is a Gaussian random process with $\langle \Delta W\rangle =0$ and $\langle \Delta W^2\rangle =\Delta t$. To expand the Euler scheme to order 1.0 of strong convergence, the double stochastic integral appearing in the remainder $\mathcal{R}$ in (24) has to be included, yielding

$x\left(t\right)=x_0+a\left(x_0\right){\displaystyle {\int }_{t_0}^t}\text{d}{t}^{\prime }+b\left(x_0\right){\displaystyle {\int }_{t_0}^t}\text{d}W\left(t^{\prime }\right)+L^{1}b\left(x_0\right){\displaystyle {\int }_{t_0}^t}{\displaystyle {\int }_{t_0}^{t^{\prime }}}\text{d}W\left(t^{\prime }\right)\text{d}W\left(t^{″}\right)+\mathcal{R}\mathrm{.}$ (26)

This is called the Milstein scheme [43]. The double stochastic integral can be calculated [33]

${\displaystyle {\int }_{t_0}^t}{\displaystyle {\int }_{t_0}^{t^{\prime }}}\text{d}W\left(t^{\prime }\right)\text{d}W\left(t^{″}\right)=\frac{1}{2}\left[{\left(W\left(t\right)-W\left(t_0\right)\right)}^2-\left(t-t_0\right)\right]\mathrm{,}$ (27)

which leads to the iteration rule for the Milstein method

$\begin{array}{c}x\left(t_{i+1}\right)=x\left(t_i\right)+\left[a\left(x\left(t_i\right)\right)-\frac{1}{2}b\left(x\left(t_i\right)\right){\partial }_{x}b\left(x\left(t_i\right)\right)\right]\Delta t\\ +b\left(x\left(t_i\right)\right)\Delta W+\frac{1}{2}b\left(x\left(t_i\right)\right){\partial }_{x}b\left(x\left(t_i\right)\right)\Delta W^2\mathrm{.}\end{array}$ (28)

With increasing order of convergence $\gamma$, the simulation schemes become more complex and include an increasing number of stochastic increments.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	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. https://doi.org/10.1126/science.1957169
[2]	Judson, P.D.L., Groom, K.M., Childs, D.T.D., Hopkinson, M., Krstajic, N. and Hogg, R.A. (2009) Microelectronics Journal, 40, 588-591. https://doi.org/10.1016/j.mejo.2008.06.037
[3]	Hartmann, S. and Elsäßer, W. (2017) Journal of Modern Optics, 7, Article No. 41866. https://doi.org/10.1038/srep41866
[4]	Pittman, T.B., Shih, Y.H., Strekalov, D.V. and Sergienko, A.V. (1995) Physical Review A, 52, R3429-R3432. https://doi.org/10.1103/PhysRevA.52.R3429
[5]	Gatti, A., Bache, M., Magatti, D., Brambilla, E., Ferri, F. and Lugiato, L.A. (2006) Journal of Modern Optics, 53, 739-760. https://doi.org/10.1080/09500340500147240
[6]	Lee, T.-P., Burrus, C.A. and Miller, B.I. (1973) IEEE Journal of Quantum Electronics, QE-9, 820-828. https://doi.org/10.1109/JQE.1973.1077738
[7]	Boitier, F., Godard, A., Rosencher, E. and Fabre, C. (2009) Nature Physics, 5, 267-270. https://doi.org/10.1038/nphys1218
[8]	Uskov, A.V., Berg, T.W. and Mørk, J. (2004) IEEE Journal of Quantum Electronics, 40, 306-320. https://doi.org/10.1109/JQE.2003.823032
[9]	Bardella, P., Rossetti, M. and Montrosset, I. (2009) IEEE Journal of Selected Topics in Quantum Electronics, 15, 785-791. https://doi.org/10.1109/JSTQE.2009.2013128
[10]	Forrest, A.F., Krakowski, M., Bardella, P. and Cataluna, M.A. (2020) Optics Express, 28, 846-859. https://doi.org/10.1364/OE.377768
[11]	Rossetti, M., Bardella, P. and Montrosset, I. (2011) IEEE Journal of Quantum Electronics, 47, 139-150. https://doi.org/10.1109/JQE.2010.2055550
[12]	Li, Z.Q. and Li, Z.M.S. (2010) IEEE Journal of Quantum Electronics, 46, 454-461. https://doi.org/10.1109/JQE.2009.2032426
[13]	Friedrich, F., Elsäßer, W. and Walser, R. (2020) Optics Communications, 458, Article ID: 124449. https://doi.org/10.1016/j.optcom.2019.124449 http://www.sciencedirect.com/science/article/pii/S0030401819307515
[14]	Hartmann, S., Friedrich, F., Molitor, A., Reichert, M., Elsäßer, W. and Walser, R. (2015) New Journal of Physics, 17, Article ID: 043039. https://doi.org/10.1088/1367-2630/17/4/043039
[15]	Jechow, A., Seefeldt, M., Kurzke, H., Heuer, A. and Menzel, R. (2013) Nature Photonics, 7, 973-976. https://doi.org/10.1038/nphoton.2013.271
[16]	Kiethe, J., Heuer, A. and Jechow, A. (2017) Laser Physics Letters, 14, Article ID: 086201. http://stacks.iop.org/1612-202X/14/i=8/a=086201 https://doi.org/10.1088/1612-202X/aa772c
[17]	Blazek, M. and Elsäßer, W. (2011) Physical Review A, 84, Article ID: 063840. https://doi.org/10.1103/PhysRevA.84.063840
[18]	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. https://doi.org/10.1088/1361-6463/ab0ea5
[19]	Forrest, A.F., Krakowski, M., Bardella, P. and Cataluna, M.A. (2019) Optics Express, 27, 10981-10990. https://doi.org/10.1364/OE.27.010981
[20]	Aho, A.T., Viheriälä, J., Virtanen, H., Zia, N., Isoaho, R. and Guina, M. (2009) Applied Physics Letters, 115, Article ID: 081104. https://doi.org/10.1063/1.5111012
[21]	Brückner, D.B., Fink, A., Schreiber, C., Röttgermann, P.J.F., Rödler, J.O. and Broedersz, C.P. (2019) Nature Physics, 15, 595-601. https://doi.org/10.1038/s41567-019-0445-4
[22]	Fang, X., Kruse, K., Lu, T. and Wang, J. (2019) Reviews of Modern Physics, 91, Article ID: 045004. https://doi.org/10.1103/RevModPhys.91.045004
[23]	Vagne, Q. and Sens, P. (2018) Physical Review Letters, 120, Article ID: 058102. https://doi.org/10.1103/PhysRevLett.120.058102
[24]	Nesti, T., Zocca, A. and Zwart, B. (2018) Physical Review Letters, 120, Article ID: 258301. https://doi.org/10.1103/PhysRevLett.120.258301
[25]	Kanazawa, K., Sueshige, T., Takayasu, H. and Takayasu, M. (2018) Physical Review Letters, 120, Article ID: 138301. https://doi.org/10.1103/PhysRevLett.120.138301
[26]	Hartmann, M.J. and Carleo, G. (2019) Physical Review Letters, 122, Article ID: 250502. https://doi.org/10.1103/PhysRevLett.122.250502
[27]	Yoshioka, N. and Hamazaki, R. (2019) Physical Review B, 99, Article ID: 214306. https://doi.org/10.1103/PhysRevB.99.214306
[28]	Qi, K., Westphal, E., Gompper, G. and Winkler, R.G. (2020) Physical Review Letters, 124, Article ID: 068001. https://doi.org/10.1103/PhysRevLett.124.068001
[29]	Liebchen, B. and Löwen, H. (2019) The Journal of Chemical Physics, 150, Article ID: 061102. https://doi.org/10.1063/1.5082284
[30]	Walser, R., Cooper, J. and Zoller, P. (1994) Physical Review A, 50, 4304. https://doi.org/10.1103/PhysRevA.50.4303
[31]	McIntyre, D., Fairchild, C., Cooper, J. and Walser, R. (1993) Optics Letters, 18, 1816-1818. https://doi.org/10.1364/OL.18.001816
[32]	Uhlenbeck, G.E. and Ornstein, L.S. (1930) Physical Review, 36, 823-841. https://doi.org/10.1103/PhysRev.36.823
[33]	Gardiner, C. (2009) Stochastic Methods: A Handbook for the Natural and Social Sciences. Springer, Berlin.
[34]	Glauber, R.J. (1963) Physical Review, 130, 2529. https://doi.org/10.1103/PhysRev.130.2529
[35]	Wiener, N. (1930) Acta Mathematica, 55, 117-258. https://doi.org/10.1007/BF02546511
[36]	Khintchine, A. (1934) Mathematische Annalen, 109, 604-615. https://doi.org/10.1007/BF01449156
[37]	Mandel, L. and Wolf, E. (1995) Optical Coherence and Quantum Optics. Cambridge Univeristy Press, Cambridge. https://doi.org/10.1017/CBO9781139644105
[38]	Kloeden, P.E. and Platen, E. (1999) Numerical Solution of Stochastic Differential Equations. Springer, Berlin.
[39]	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. http://dlmf.nist.gov
[40]	Jerri, A.J. (1998) The Gibbs Phenomenon in Fourier Analysis, Splines and Wavelet Approximations. Springer US, New York. https://doi.org/10.1007/978-1-4757-2847-7
[41]	Süßmann, G. (1997) Zeitschrift für Naturforschung, 52, 49-52. https://doi.org/10.1515/zna-1997-1-214
[42]	Schleich, W.P. (2011) Quantum Optics in Phase Space. Wiley, Berlin.
[43]	Milstein, G.N. (1975) Theory of Probability and Its Applications, 19, 557-562. https://doi.org/10.1137/1119062