Statistics of the Optical Phase of a Gain-Switched Semiconductor Laser for Fast Quantum Randomness Generation

Abstract

The statistics of the optical phase of the light emitted by a semiconductor laser diode when subject to periodic modulation of the applied bias current are theoretically analyzed. Numerical simulations of the stochastic rate equations describing the previous system are performed to describe the temporal dependence of the phase statistics. These simulations are performed by considering two cases corresponding to random and deterministic initial conditions. In contrast to the Gaussian character of the phase that has been assumed in previous works, we show that the phase is not distributed as a Gaussian during the initial stages of evolution. We characterize the time it takes the phase to become Gaussian by calculating the dynamical evolution of the kurtosis coefficient of the phase. We show that, under the typical gain-switching with square-wave modulation used for quantum random number generation, quantity is in the ns time scale; that corresponds to the time it takes the system to lose the memory of the distribution of the initial conditions. We compare the standard deviation of the phase obtained with random and deterministic initial conditions to show that their differences become more important as the modulation speed is increased.

1. Introduction

Experimental and theoretical understanding of the fluctuations of laser light began shortly after the invention of the laser [1,2,3,4,5]. Special attention has been devoted to fluctuations of the light emitted by semiconductor lasers [6,7,8,9,10] due to their vast variety of applications. The best available theoretical description of these fluctuations is based on the Fokker–Planck equation, or alternativelly on Langevin’s stochastic rate equations [3,6,7,8,11]. The phase of the laser electrical field is a random quantity, mainly due to the effect of the spontaneous emission noise. The random character of this phase is precisely the basis of some of the available methods for quantum random number generation (QRNG).

Random numbers are a vital resource for numerous applications including cryptography, statistical analysis, stochastic simulations, decission making in engineering processes, quantitative finance, gambling, massive data processing, etc. [12,13]. Random number generators (RNG) use software algorithms (pseudorandom number generators) or hardware physical devices. Typical physical processes used to generate random numbers are radioactive decay, Johnson or Zener’s noise, chaos noise [13,14] and quantum phenomena [12,13]. QRNGs are a particular case of physical RNGs that can generate truly random numbers from the fundamentally probabilistic nature of quantum events [13]. The advantage of using QRNGs relies on its unpredictibility, which can be proven to be based on quantum physics laws. Typical QRNGs are based on quantum optics [13]. These generators can be divided in (i) generators that use single-photon sources, and (ii) generators that use multi-photon sources, typically semiconductor lasers or LEDs. QRNGs based on single-photon detection methods include: Branching path generators [15], generators measuring the time of arrival of photons [16], photon counting generators [17], and attenuated pulse generators [18]. These methods have been experimentally compared in [19]. Multi-photon QRNGs include: Generators based on quantum vacuum fluctuations [20], on amplified spontaneous emission (ASE) signals [21,22], and on phase noise in continuous wave [23,24,25] and in pulsed laser diodes [26,27,28,29,30,31].

Spontaneous emission is a useful mechanism to generate quantum fluctuations, as it can be ascribed to the vacuum fluctuations of the optical field [26]. Randomness due to spontaneous emission is the basis of QRNGs based on pulsed single-mode laser diodes [26,27,28,30,31]. These generators have several advantages. They are made of commercially available components: For instance, standard photodetectors can be used due to the high signal levels. They are simple, low-cost, robust, and fast: Generation rates up to 43 Gbps quantum random bit generation have been experimentally demonstrated [27]. In these QRNGs the laser diode is periodically modulated from below to above threshold in such a way that gain-switching operation is obtained, typically at Gbps rates. While the laser is below threshold the optical phase becomes random due to the spontaneous emission noise. Gain-switching operation produces a periodic train of laser pulses with random phases. Phase fluctuations are then converted into amplitude fluctuations by using interferometric setups [27,31]. Detection and filtering of the amplitude fluctuations provides the generation of random values with an almost uniform distribution.

The applications of QRNGs, for instance in cryptography [31,32], require that the physical processes underlying their operation must be properly understood and described. For QRNGs based on pulsed semiconductor lasers, it is essential an accurate description of the phase diffusion process, that is, laser phase fluctuations must be qualitatively and quantitatively characterized. Modelling of these fluctuations has been performed by numerical integration of the laser stochastic rate equations [27,30,31,33,34,35,36]. Good quantitative agreement between experiments and theory is achieved when the complete set of parameters of the rate equations is known for the specific laser diode. Good agreement between experimental and theoretical phase fluctuations has been recently reported for a discrete mode laser (DML) [36] for which a complete extraction of the intrinsic parameters was performed [35,37]. This permits a quantitative description of the dependence of phase diffusion on the laser and modulation parameters. On the qualitative side, statistics of optical phase has been described as Gaussian in numerical simulations [27,33,34] since spontaneous emission noise has also Gaussian distribution. However, in these generators the bias current is periodically modulated in such a way that the evolution is mainly in a transient regime, specially when operating at fast bit rates. It is then expected that the choice of initial conditions in the simulations must have impact on the statistics of the optical phase and on the time it takes the phase to be distributed as a Gaussian. This is in fact the main objective of this work: The investigation of the conditions for which the phase becomes distributed as a Gaussian.

In this paper we report a theoretical study of the phase diffusion in gain-switched semiconductor lasers. This is done by performing numerical simulations of the stochastic rate equations for the complex electrical field and carrier number. In our modelling we use the set of parameters recently extracted for a DML device. With these parameters we first analyze the impact of the carrier noise on the phase statistics. In the rest of the paper we focus on the calculation of the temporal dependence of the statistical moments and distribution of the phase. We first consider random initial conditions that contrast to previous analysis in which deterministic fixed initial conditions were chosen [34]. We compare the phase statistics obtained for both types of initial conditions. For both cases we show that the phase is not distributed as a Gaussian because of the non-Gaussianity of the initial conditions. This contrasts to the Gaussian character of the phase that has been assumed in previous works [27,33,34]. We characterize the time it takes the phase becomes approximately Gaussian by calculating the temporal evolution of the kurtosis coefficient of the phase. Our calculations indicate that under the typical gain-switching with square-wave modulation used for QRNG, the time it takes to the phase to become Gaussian is in the ns scale. These are the typical times for which the memory of the distribution of the initial conditions is lost. The comparison between the variance of the phase obtained with random and fixed initial conditions show that their differences become more important as the modulation speed is increased.

Our paper is organized as follows. In Section 2, we present our theoretical model. Section 3 is devoted to analyze the dynamical evolution of the relevant variables. In Section 4, we study the temporal evolution of the phase statistics. Finally, in Section 5 we discuss and summarize our results.

2. Theoretical Model

Gain-switched single-mode semiconductor laser dynamics can be modelled by using a set of stochastic rate-equations that read (in Ito’s sense) [6,37,38]

(1)$\begin{array}{ccc}\hfill \frac{dP}{dt}& =& \left[\frac{G_N(N-N_t)}{1+ϵP}-\frac{1}{{\tau }_p}\right]P+R_{sp}\left(N\right)+\sqrt{2R_{sp}\left(N\right)P}{F}_p\left(t\right)\hfill \end{array}$

(2)$\begin{array}{ccc}\hfill \frac{d\varphi }{dt}& =& \frac{\alpha }{2}\left[G_N(N-N_t)-\frac{1}{{\tau }_p}\right]+\sqrt{\frac{R_{sp}\left(N\right)}{2P}}{F}_{\varphi }\left(t\right)\hfill \end{array}$

(3)$\begin{array}{ccc}\hfill \frac{dN}{dt}& =& \frac{I\left(t\right)}{e}-(AN+B{N}^2+C{N}^3)-\frac{G_N(N-N_t)P}{1+ϵP}\hfill \end{array}$

where $P\left(t\right)$ is the number of photons inside the laser, $\varphi \left(t\right)$ is the optical phase, and $N\left(t\right)$ is the number of carriers in the active region. The parameters appearing in these equations are the following: $G_N$ is the differential gain, $N_t$ is the number of carriers at transparency, $ϵ$ is the non-linear gain coefficient, ${\tau }_p$ is the photon lifetime, $\alpha$ is the linewidth enhancement factor, e is the electron charge, and $A,B$ and C are the non-radiative, spontaneous, and Auger recombination coefficients, respectively. In these equations we consider a temporal dependence of the injected current, $I\left(t\right)$, and a rate of the spontaneous emission given by $R_{sp}\left(N\right)=\beta B{N}^2$ where $\beta$ is the fraction of spontaneous emission coupled into the lasing mode. The Langevin terms $F_P\left(t\right)$ and $F_{\varphi }\left(t\right)$ in Equations (1) and (2), represent fluctuations due to spontaneous emission, with the following correlation properties, $<F_i\left(t\right)F_j\left(t^{\prime }\right)>={\delta }_{ij}\delta (t-t^{\prime })$, where $\delta \left(t\right)$ is the Dirac delta function and ${\delta }_{ij}$ the Kronecker delta function with the subindexes i and j referring to the variables P and $\varphi$.

QRNG systems based on gain-switching of single-mode laser diodes are such that a large signal modulation of the bias current is applied to the device. We consider an injected current following a square-wave modulation of period T with $I\left(t\right)=I_{on}$ during $T/2$, and $I\left(t\right)=I_{off}$ during the rest of the period. This modulation is such that $I_{off}<I_{th}$, where $I_{th}$ is the threshold current of the laser, for obtaining a random evolution of the optical phase induced by the spontaneous emission noise. Numerical integration of the previous stochastic rate equations by usual Euler–Maruyama [3,39] or Heun’s predictor-corrector algorithms [37] present instabilities when the photon number is very small, a situation always present in this type of QRNGs: some spontaneous noise events cause negative values of P that lead to numerical instabilities due to the square root factor multipliying the noise term in Equations (1) and (2). Very recently a set of rate equations for the complex electrical field, $E\left(t\right)$, instead of equations for P and $\varphi$ has been proposed to solve this problem [35]. These equations are the following:

(4)$\begin{array}{ccc}\hfill \frac{dE}{dt}& =& \left[\left(\frac{1}{1+ϵ\mid E{\mid }^2}+i\alpha \right)G_N(N-N_t)-\frac{1+i\alpha }{{\tau }_p}\right]\frac{E}{2}+\sqrt{\beta B}N\xi \left(t\right)\hfill \end{array}$

(5)$\begin{array}{ccc}\hfill \frac{dN}{dt}& =& \frac{I\left(t\right)}{e}-(AN+B{N}^2+C{N}^3)-\frac{G_N(N-N_t)\mid E{\mid }^2}{1+ϵ\mid E{\mid }^2}\hfill \end{array}$

where $E\left(t\right)=E_1\left(t\right)+i{E}_2\left(t\right)$ is the complex electrical field and $\xi \left(t\right)={\xi }_1\left(t\right)+i{\xi }_2\left(t\right)$ is the complex Gaussian white noise with zero average and correlation given by $<\xi \left(t\right){\xi }^{*}\left(t^{\prime }\right)>=\delta (t-t^{\prime })$ that represents the spontaneous emission noise, and where we have considered that $R_{sp}\left(N\right)=\beta B{N}^2$. These equations exactly correspond to our initial model because the application of the rules for the change of variables in the Ito’s calculus [11] to $P= \mid E{\mid }^2=E_1^2+E_2^2$ and $\varphi =arctan(E_2/E_1)$ in Equations (4) and (5) gives Equations (1)–(3). Instabilities do not appear because P is not inside the square root factor that multiplies the noise term.

Up to now we have considered an equation for N that has not any noise term. Carrier noise can also be important for describing statistics in semiconductor laser dynamics [6]. These fluctuations can be taken into account if we substitute Equation (5) by

(6)$\begin{array}{c}\frac{dN}{dt}=\frac{I\left(t\right)}{e}-(AN+B{N}^2+C{N}^3)-\frac{G_N(N-N_t)\mid E{\mid }^2}{1+ϵ\mid E{\mid }^2}\hfill \\ \hfill +\sqrt{2\left(AN+B{N}^2+C{N}^3+\frac{I\left(t\right)}{e}\right)}{\xi }_N-2\sqrt{\beta B}N\left(E_1{\xi }_1+E_2{\xi }_2\right)\end{array}$

where ${\xi }_N$ is a real Gaussian white noise of zero average and correlation given by $<{\xi }_N\left(t\right){\xi }_N\left(t^{\prime }\right)>=\delta (t-t^{\prime })$ and statistically independent of $\xi \left(t\right)$ [6,10,37,40].

In this work we will numerically solve Equations (4) and (6) by using the Euler–Maruyama algorithm [3,39] with an integration time step of 0.001 ps. We will use the numerical values of the parameters that have been extracted for a discrete mode laser (DML) [35,37]. This device is a single longitudinal mode semiconductor laser emitting close to 1550 nm wavelength and $I_{th}=14.14$ mA at a temperature of 25 ${}^{\circ }$C. The values of the parameters are $G_N=1.48\times {10}^4$ s${}^{-1}$, $N_t=1.93\times {10}^7$, $ϵ=7.73\times {10}^{-8}$, ${\tau }_p=2.17$ ps, $\alpha =3$, $\beta =5.3\times {10}^{-6}$, $A=2.8\times {10}^8$ s${}^{-1}$, $B=9.8$ s${}^{-1}$, and $C=3.84\times {10}^{-7}$ s${}^{-1}$ [35,37]. Simulation and experimental results have shown not only qualitative but also a remarkable quantitative agreement for a very wide range of gain-switching conditions [35,37,41].

3. Analysis of the Dynamics

We first analyze the dynamical evolution of relevant variables when the laser is modulated with $I_{on}=30$ mA, $I_{off}$ = 7 mA, and $T=1$ ns. The laser is switched-off with a current close to $I_{th}/2$, for obtaining a significant effect of the spontaneous emission noise on the randomness of the phase. Figure 1a–c show the photon number, carrier number, and optical phase, respectively, as a function of time. We integrate the equations for consecutive bias current pulses in such a way that the initial conditions for one period correspond to the final values of the variables at the end of the previous period. Figure 1a shows P for three consecutive pulses. The laser is switched-on with $I_{on}$ at $t=1$ ns. After this time P begins to build-up from very small random values determined by the spontaneous emission noise events. After the emission of the pulse with the corresponding relaxation oscillations, P begins to decrease at $t=1.5$ ns (when $I_{off}$ is applied), reaching the small random values at which spontaneous emission noise dominates the device dynamics. N begins at $t=1$ ns from a value well below the threshold carrier number, $N_{th}=N_t+1/\left(G_N{\tau }_p\right)=5.045\times {10}^7$, as it can be seen in Figure 1b. The characteristics relaxation oscillations of N associated to the pulse emission are followed by a monotonous decrease from $t=1.5$ ns to 2 ns due to the value below threshold of $I_{off}$.

The optical phase is calculated at each integration step from $E_1$ and $E_2$ in such a way that it is a continuous function of t. The dynamical evolution of $\varphi$ is shown in Figure 1c. When P is large (small) the noise term in Equation (2) is much smaller (larger) than the other term in that equation and $\varphi$ mainly evolves in a deterministic (random) way. The deterministic decrease of $\varphi$ is due to the value below threshold of the current when switching-off the laser: $G_N(N-N_t)-\frac{1}{{\tau }_p}<0$ because $N<N_{th}$, and therefore $\varphi$ decreases (see Equation (2)).

Visualization of different random trajectories and calculation of statistical moments of the phase, specially its standard deviation, ${\sigma }_{\varphi }\left(t\right)$, have been usually done by overlaying them in a temporal window with a duration of a few periods [33,34,35]. For instance just one period is considered in references [34,35] to calculate the value of ${\sigma }_{\varphi }\left(t\right)=\sqrt{<{\varphi }^2>\left(t\right)-<\varphi {>}^2\left(t\right)}$ with $0\le t\le T$. To obtain well defined averages, $<\varphi >\left(t\right)$ and $<{\varphi }^2>\left(t\right)$, it is necessary to make a choice of the initial conditions at the beginning of each period because $\varphi$ is an unbounded quantity, as shown in Figure 1. One choice is to take $P\left(0\right)=<P\left(0\right)>$, $N\left(0\right)=<N\left(0\right)>$, and $\varphi \left(0\right)=<\varphi \left(0\right)>$ [34], that is fixed initial conditions. A second choice is to take random initial conditions [35]. Photon and carrier numbers at $t=0$ are those obtained at the end of the previous period, like in Figure 1. The change with respect to Figure 1 is related to the phase and it is based on the cyclic nature of angles: We consider that $\varphi$ at the beginning of the period, $\varphi \left(0\right)$, is that corresponding to $\varphi$ at the end of the previous period, $\varphi \left(T\right)$, but converted into the $[0,2\pi )$ range, that is we consider that $\varphi \left(0\right)$ is given by $\varphi \left(T\right)-\mathrm{int}\left(\frac{\varphi \left(T\right)}{2\pi }\right)2\pi$.

Figure 2 shows the temporal evolution of P, N and $\varphi$, plotted in a window of duration T, corresponding to the three consecutive pulses of Figure 1 and using the previous choice of random initial conditions. Figure 2a,c show that laser pulses that have a larger switch-on time, defined as the time at which P crosses a fixed level, have also a larger value of the maximum of N and P [9]. Figure 2b shows that $\varphi$ takes values in a range of several multiples of 2$\pi$ during one period. Figure 2b also shows, in a more clear way than in Figure 1, that the phase fluctuations are more important at the beginning and at the end of the pulse. Comparison between Figure 2a,b shows that pulses with a similar evolution of P can have a very different phase evolution (see black and red realizations). In the next section we will focus on the description of the temporal evolution of the phase statistics.

4. Analysis of the Phase Statistics

The dynamical evolution of the variance of the phase, ${\sigma }_{\varphi }^2$, is shown in Figure 2d for the case of random initial conditions and a temporal window of duration $T=$ 1 ns. ${\sigma }_{\varphi }^2\left(t\right)$ has been calculated by averaging over 2 × 10${}^4$ temporal windows. ${\sigma }_{\varphi }^2\left(0\right)>0$ because of our choice of random initial conditions. Large increases of ${\sigma }_{\varphi }^2\left(t\right)$ occur while P is small and dominated by the spontaneous emission noise, that is at the beginning and at the end of the period. While the evolutions of P and $\varphi$ are deterministic and $I>I_{th}$ (0.15 ns < t < 0.5 ns) ${\sigma }_{\varphi }^2\left(t\right)$ oscillates with the frequency of the relaxation oscillations around a value that increases linearly with time, similarly to what was observed by Henry [8]. These oscillations and the linear increase are barely seen in Figure 2d because of the vertical scale determined by the large values of the variance when the laser is switched-off. From 0.5 ns < t < 0.65 ns, while $\varphi$ still has a deterministic evolution, there is a slight decrease of ${\sigma }_{\varphi }^2\left(t\right)$. After that time, both $\varphi$ and P become determined by the spontaneous emission noise. In this way the linear increase of ${\sigma }_{\varphi }^2\left(t\right)$ with t, characteristic of phase diffusion, is observed until the end of the period, as it is seen in Figure 2d.

We now analyze the effect of carrier noise on the statistics of the phase. Figure 3 shows the probability density function (pdf) of $\varphi$ at three different times when the carrier noise is considered (that is, integrating Equation (6)) and when it is neglected (considering instead Equation (5)). This figure has been obtained using the same conditions of Figure 2.

Figure 3 shows that the effect of carrier noise on the statistics of $\varphi$ is very small. In fact, it has been shown that the consideration of noise in the carrier equation is not important during transient regimes [9,33], being only essential in the stationary regime for calculating quantities like the relative intensity noise [6]. Figure 3 also shows the Gaussian distributions of average and standard deviation given by the simulation with carrier noise. It is clear that the Gaussian distribution does not describe well the phase satististics, specially for short times ($t=0.1$ and $t=0.5$ ns). The Gaussian approximation becomes better at longer times ($t=0.9$ ns).

A way of quantifying if the Gaussian distribution is suitable for describing the phase statistics is by calculating moments of $\varphi$ of order higher than 2. Asymmetry and kurtosis coefficient of the simulated data are shown in Figure 4 as a function of time. Both coefficients must vanish if the distribution is Gaussian. Figure 4a shows with black lines the asymmetry, ${\gamma }_r$, and kurtosis, ${\kappa }_r$, coefficients obtained under the same conditions of Figure 2, that is with random initial conditions. While the phase distribution is symmetric (${\gamma }_r\sim 0$), ${\kappa }_r$ is significantly larger than zero. ${\kappa }_r$ decreases fast until it develops a small peak close to the time at which the first relaxation oscillation appears. After that peak it reaches a plateau that finishes when P reaches the spontaneous emission noise level (around $t=$ 0.7 ns). From that time $\varphi$ diffuses and ${\kappa }_r$ monotonously decreases reaching values that are closer to zero at the end of the period (${\kappa }_r=$ 0.65 at $t=0.9$ ns). Figure 4b shows ${\gamma }_r$ and ${\kappa }_r$ when $T=$ 2 ns. In this case $\varphi$ has more time to diffuse when the laser is switched-off and then the Gaussian approximation is better at the end of the period (${\kappa }_r=$ 0.14 at $t=2$ ns).

The reason why $\varphi$ is not Gaussian can be understood by plotting the pdf of $\varphi$ at $t=0$. Figure 5 shows that distribution for the case of $T=1$ ns. The distribution corresponds to a uniform random variable in $[0,2\pi )$. This is because of the way random initial conditions are chosen: Doing the operation $\varphi \left(0\right)=\varphi \left(T\right)-\mathrm{int}\left(\frac{\varphi \left(T\right)}{2\pi }\right)2\pi$ from a broad nearly Gaussian distribution for $\varphi \left(T\right)$ makes $\varphi \left(0\right)$ a uniform random variable, U$(0,2\pi )$. The kurtosis of U$(0,2\pi )$ is $354/5\sim 70.8$. Diffusion of $\varphi$ at the beginning of the period (see Figure 2) makes ${\kappa }_r$ to decrease quickly, but not enough for becoming strictly Gaussian, even at the end of the period.

Of course these results depend on the way initial conditions are chosen. Another way of choosing these values is by considering fixed initial values for $E\left(0\right),$ and $N\left(0\right)$. Figure 4 shows, with red lines, asymmetry and kurtosis coefficients, ${\gamma }_f$ and ${\kappa }_f$, when fixed initial conditions are used. We choose these values in the following way. We first integrate Equations (4) and (6) from arbitrary initial conditions corresponding to below threshold operation in order to find the averaged $<P\left(t\right)>,<N\left(t\right)>,$ and $<\varphi \left(t\right)>$ for $0\le t\le T$. Then we choose $N\left(0\right)=<N\left(T\right)>,$ and $E\left(0\right)=\sqrt{<P\left(T\right)>}(cos<\varphi \left(T\right)>+isin<\varphi \left(T\right)>)$. This election is similar to that considered in [34]. Figure 4 shows that the evolution of ${\gamma }_f$ and ${\kappa }_f$ is very similar to that of ${\gamma }_r$ and ${\kappa }_r$, respectively. ${\kappa }_f>{\kappa }_r$ because the initial delta-like distribution of $\varphi \left(0\right)$ produce larger values of the kurtosis. These differences decrease with t, specially when spontaneous emission dominates the phase evolution: In Figure 4a,b when $t>0.7$ ns ($t>1.2$ ns).

The choice of initial conditions also impacts on the values of the standard deviation as a function of t. In Figure 6a the dynamical evolution of ${\sigma }_{\varphi }$ for both, random and fixed initial conditions, is shown when $T=1$ ns. Again both standard deviations have similar trends but the value for random initial conditions is larger than that obtained for the fixed ones. This is due to the non-zero value of ${\sigma }_{\varphi }\left(0\right)$ obtained with the uniform distribution of $\varphi \left(0\right)$ in contrast to the zero value obtained for fixed initial conditions. Relative differences between both quantities enhance if the speed of QRNG increases as it can be seen in Figure 6b where results obtained for $T=0.4$ ns have been plotted. For instance, ${\sigma }_{\varphi }$ at 0.2 ns is around 20 % smaller for the case of fixed initial conditions.

The dependence of the phase statistics on the way initial conditions are chosen suggests that averages must be done in a different way in order to lose the memory of those initial conditions. We have been considering averages performed in a temporal window with the same duration than the period of the current, T. From now on we will consider longer temporal windows for calculating statistical averages. Figure 7 illustrates the situation found when averages are calculated over a window of duration 2T. Random initial conditions are considered such that $\varphi \left(0\right)=\varphi \left(2T\right)-\mathrm{int}\left(\frac{\varphi \left(2T\right)}{2\pi }\right)2\pi$. Averages have been done over 2 × 10${}^4$ 2T-windows, where T = 1 ns, in order to compare with situations illustrated in previous figures. Figure 7a shows the averaged phase vs. t. The drift towards decreasing values of the phase is similar to that shown in Figure 1c. Standard deviation and variance of the phase are shown in Figure 7b,c, respectively. $<\varphi \left(t\right)>$, ${\sigma }_{\varphi }\left(t\right)$ and ${\sigma }_{\varphi }^2\left(t\right)$ during the second half of the $2T$-window are basically replicas of what was found in the first half. The continuity of $\varphi$ along the 2T-window makes that ${\sigma }_{\varphi }\left(t\right)$ and ${\sigma }_{\varphi }^2\left(t\right)$ monotonously increase. However the situation is different when considering the kurtosis coefficient as Figure 7d shows. In this case, during the second half of the window ${\kappa }_r$ keeps on decreasing towards the zero value. This means that the distribution of the phase keeps on approaching to the Gaussian shape. In fact ${\kappa }_r$ = 0.22 when t = 2 ns.

That approach can be illustrated by plotting the phase pdf at two different times. Figure 8 shows those distributions at times $t=0$ and $t=1.1$ ns. The phase at $t=0$ is a U$(0,2\pi )$ random variable, similarly to Figure 5. The phase at $t=1.1$ ns is approximately Gaussian as it can be seen when comparing with the Gaussian of average and standard deviation obtained from simulations. The kurtosis coefficient in Figure 8b is 0.4. Figure 8b can also be compared with the pdf obtained at t = 0.1 ns in Figure 3b because both distributions correspond to 0.1 ns after switching-on the bias current. The pdf in Figure 3b is not Gaussian while the pdf in Figure 8c is approximately Gaussian. This indicates that in order to have a phase distributed as a Gaussian it is necessary to calculate and average the phase in windows with durations of several modulation periods. In this way the memory of the initial conditions and their distribution is lost.

5. Discussion and Summary

In our study we have considered two types of initial conditions, corresponding to deterministic and random values of the variables. Fixed initial conditions have been considered because they have been used in previous studies of QRNG. They are not the best choice for simulation of these systems because the spontaneous emission noise, that is always present in the system, causes fluctuations in the variables of the system at all times. These include the times at which each period begins, and so initial conditions must be also random, as it is also expected in an experimental realization of the system. We have chosen these random initial values by calculating the phase angle in the $[0,2\pi )$ range that corresponds to the final value in the previous averaging window. Note that the conversion to the $[0,2\pi )$ range is necessary if a calculation of well defined statistical moments of the phase is required. If no conversion is done, not even $<\varphi \left(t\right)>$ could be calculated because $\varphi$ decreases in each averaging window in a magnitude of more than several $2\pi$, as illustrated for instance in Figure 1c.

Deterministic initial conditions and phase averages over windows of T-duration have been recently used for describing the phase statistics [34]. While these conditions can give an approximation to the phase distribution and their statistical moments, our results show that it is necessary to consider averages over windows of several T-duration and random initial conditions for obtaining Gaussian statistics for the phase at the end of the averaging period.

We now briefly discuss the effect of two laser parameters, the non-linear gain and the Auger coefficients, on the standard deviation of the phase. The number of relaxation oscillation peaks increases when the non-linear gain coefficient decreases. The standard deviation of the phase at the end of the modulation period oscillates when changing $I_{on}$ [35]. The number of these oscillations is directly related to the number of relaxation oscillation peaks that are excited. In this way, the main effect of having a small nonlinear gain coefficient is to observe more oscillations of the standard deviation of the phase as a function of $I_{on}$. The effect of the Auger coefficient is also important for describing the standard deviation of the phase. In fact we have shown that the Auger term must be considered in the carrier recombination term for achieving good agreement between experiments and theory [36].

Summarizing, we have theoretically analyzed the phase diffusion in gain-switched semiconductor lasers by performing numerical simulations of the corresponding stochastic rate equations. We have focused on the calculation of the temporal dependence of the statistical moments and distribution of the phase. We have considered several types of initial conditions for the phase. By using the temporal dependence of the kurtosis coefficient we have shown that the phase pdf becomes Gaussian only after the memory of the statistical distribution of the initial conditions is lost. We show that under the typical gain-switching with square-wave modulation used in QRNGs, the time it takes for the phase to become Gaussian is in the ns scale. We finally compared the variance of the phase obtained with random and fixed initial conditions to show that their differences are more important as the modulation speed is increased. This is precisely the situation in which faster generation bit rates are achieved when using QRNGs based on gain-switched laser diodes.

Funding

This research was funded by Ministerio de Economía y Competitividad (MINECO/FEDER, UE), Spain under grant RTI2018-094118-B-C22.

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

Not applicable.

Conflicts of Interest

The author declares no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:

QRNG	Quantum random number generation
DML	Discrete mode laser
PDF	Probability density function
ASE	Amplified spontaneous emission
RNG	Random number generation
LED	Light emitting diode

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures

Figure 1. (a) Photon number, (b) carrier number, and (c) optical phase as a function of time for three consecutive pulses when T= 1 ns.

Figure 2. (a) Photon number, (b) phase, and (c) carrier number dynamical evolution for three different realizations are shown with black, red, and green lines in a temporal window of duration T. (d) Variance of the phase as a function of t. In this figure T= 1 ns and the three realizations are extracted from the time traces of Figure 1.

Figure 3. Probability density function of the phase at three different times in (a) linear, and (b) logarithmic vertical scale. Pdfs obtained with and without noise in the carrier number equation are plotted with red and black solid lines. Gaussian approximations are plotted with blue dashed lines.

Figure 4. Asymmetry and kurtosis coefficients of the phase as a function of time for (a) T = 1 ns, and (b) T = 2 ns. Asymmetry and kurtosis coefficients are plotted with solid and dashed lines, respectively. Results for random and fixed initial conditions are plotted with black and red lines, respectively.

Figure 5. Probability density function of the initial value of the phase for T = 1 ns.

Figure 6. Standard deviation of the phase as a function of time for (a) T = 1 ns, and (b) T = 0.4 ns. Results for random and fixed initial conditions are plotted with black and red lines, respectively.

Figure 7. (a) Average, (b) standard deviation, (c) variance, and (d) kurtosis coefficient of the phase as a function of time for a two-period window with T = 1 ns.

Figure 8. Probability density function of the phase at (a) t = 0 ns, and (b) t= 1.1 ns for a two-period window with T = 1 ns. The Gaussian approximation is plotted with a blue dashed line.