Spontaneous emission rate and phase diffusion in gain-switched laser diodes

Optics & Laser Technology 150 (2022) 107992 Available online 26 February 2022 0030-3992/© 2022 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). Spontaneous emission rate and phase diffusion in gain-switched laser diodes Ana Quirce, Angel Valle * Instituto de Física de Cantabria (CSIC-Univ. Cantabria), Avda. Los Castros s/n, E39005 Santander, Spain ARTICLE INFO Keywords: Laser noise Phase noise Random number generation Semiconductor lasers 2010 MSC: 00-01 99-00 ABSTRACT Quantum random number generation (QRNG) has become a topic of growing interest in recent years due to important applications in cryptography and simulations. Interferometric detection of phase diffusion in gainswitched single-mode semiconductor lasers is one of the main generation techniques. In this paper, we study experimentally and theoretically the phase diffusion in gain-switched discrete mode laser diodes. We derive a stochastic rate equations model for the laser electric field that avoids numerical instabilities that appear when simulating amplitude and phase equations. Measurements are performed in order to extract the parameters of our semiconductor laser. Spontaneous emission rate coupled into the lasing mode is measured as a function of the carrier number for bias currents below threshold. A quadratic dependence is obtained that permits us to evaluate the validity of the linear approximations that have been used to describe laser phase diffusion in QRNG. The good agreement between experiments and theory permits us to give a realistic quantitative description of the dynamical evolution of the phase statistics in this type of QRNGs. 1. Introduction Random numbers are essential for many applications including cryptography, Monte Carlo simulations, massive data processing, stock market prediction, gambling, etc. [1,2]. Typical random number generators (RNG) use software algorithms (pseudorandom number generator, PRNG) or hardware physical devices. Examples of physical processes utilized to generate randomness include: Johnson’s noise, Zener noise, radioactive decay, chaos noise [3,4], stochastic pulse-topulse fluctuation in the supercontinuum [5], and quantum phenomena [1]. Quantum random number generators (QRNGs) are a particular case of physical RNGs in which random numbers are extracted from quantum events. The main advantage of QRNGs is that their inherent quantum mechanical processes are the best guarantee for offering optimum security and privacy while maintaining high performance [2,6]. Quantum optics is the basis of most of the existing QRNGs [2]. Single-photon sources, light emitting diodes or laser diodes are used in these generators. QRNGs based on single-photon detection methods include: generators measuring photon arrival time [7], branching path generators [8], attenuated pulse generators [9], and photon counting generators [10]. A comparison between some of these methods can be found in [11]. Multiphoton QRNGs have also been proposed and demonstrated. Some examples include generators based on amplified spontaneous emission (ASE) signals [12,13], and on quantum vacuum fluctuations [14]. Also special interest has been focused on generators using semiconductor lasers as sources. A first approach measures the phase noise in continuous wave operation [15–17]. A second approach considers fluctuations of the light when the laser bias current is modulated in a large signal regime[18–22]. In these last systems high bit-rate random numbers are obtained by interferometric detection of pulses generated by gain-switched singlemode semiconductor lasers [6,18–20,22]. Fast rates, up to 43 Gbps random bit generation, have been demonstrated [19]. The current applied to a single-mode laser diode is periodically modulated from a well below threshold value to a value above threshold. While the laser is below threshold the phase of the laser fluctuates due to spontaneous emission noise, which is quantum mechanical by nature. Periodic gain switching provides a pulse train in which each laser pulse has a random phase. Phase fluctuations are converted into amplitude fluctuations by using interference effects in an unbalanced Mach–Zehnder interferometer, with a delay that matches the pulse repetition period. Detection and filtering of the amplitude is performed to generate random data with an almost uniform distribution. Advantages of QRNGs based on pulsed gain-switching of single-mode laser diodes include simplicity, low-cost, * Corresponding author. E-mail address: [email protected] (A. Valle). Contents lists available at ScienceDirect Optics and Laser Technology journal homepage: www.elsevier.com/locate/optlastec https://doi.org/10.1016/j.optlastec.2022.107992 Received 29 September 2021; Received in revised form 14 January 2022; Accepted 19 February 2022 Optics and Laser Technology 150 (2022) 107992 2 high-speed, robustness, and operation with flexible clock frequencies [19]. They can be made of commercially available components with high signal levels that permit the use of standard photodetectors. Moreover indium phosphide photonic integrated circuit versions have been demonstrated [6]. The understanding and description of the physical processes underlying the operation of a QRNG are essential for achieving a maximum level of security, specially for quantum key distribution systems [23]. Therefore a QRNG must be precisely characterized and measured. For the case of QRNGs based on laser phase diffusion due to spontaneous emission noise, it is essential a good quantitative description of the dynamical evolution of the statistics of that phase. This description can be made when the complete set of parameters of the stochastic rate equations is known for the specific laser device. In this way we can obtain a realistic dependence of the phase diffusion on the laser and modulation parameters. One of the key ingredients in the Langevin-type rate equations that simulate this type of QRNGs is the spontaneous emission noise term since it is the term responsible for their randomness. The spontaneous emission rate that appears in that term has been assumed to depend linearly on the carrier density [19,22–25]. However a more realistic description for that dependence is that corresponding to the bimolecular recombination [26], that is a quadratic dependence. Also stochastic rate equations for the photon number and optical phase of the laser have been considered [19,22–25,27,28]. Numerical solution of these equations is problematic when the bias current is below threshold because negative photon numbers can eventually appear causing numerical instabilities [25]. Recently, parameter extraction of a single longitudinal mode laser has been performed[27] for describing gain-switching operation of that device. A very good quantitative agreement between numerical and experimental results has been shown for a very wide operation range [27,28]. In this paper, we report an experimental and theoretical study of the phase diffusion in gain-switched semiconductor lasers. We perform the experiments using the single-mode discrete mode laser used in [28]. In the first part we obtain an electric field rate equations model that has the advantages of i) avoiding the above mentioned numerical instabilities, and ii) describing the evolution of the optical phase without any kind of saturation term. The measurement of the spectral linewidth and emitted power as a function of the below threshold bias current, and the use of the Schawlow-Townes formula permit us to find the spontaneous emission rate as a function of the carrier number. We find that the quadratic dependence describes well the experimental results. Also we will evaluate the validity of the linear approximations that have been used to describe the phase diffusion. The paper is organized as follows. In Section 2 we derive the theoretical model. Section 3 is devoted to the extraction of the laser parameters and the analysis of the experimental results. In Section 4 we theoretically analyze the phase diffusion obtained with the extracted parameters and different types of spontaneous emission rates, and finally in Section 5 we discuss and summarize our results. 2. Theoretical model The dynamics of a gain-switched single-mode semiconductor laser can be modelled by using a set of rate-equations. These stochastic differential equations read (in Ito’s sense) [26,27,29] dp dt =[Γvgg(n) 1+∊p−1 τ p]p+ΓRsp(n)+ 2ΓRsp(n)p √Fp(t)(1) dϕ dt = α 2[Γvgg(n)−1 τ p]+ ΓRsp(n) 2p √Fϕ(t)(2) dn dt =I(t) eVa −(an +bn2+cn3)−vgg(n)p 1+∊p(3) where p(t)is the photon density, ϕ(t)is the optical phase in the reference frame corresponding to the resonant frequency at the threshold current [30], and n(t)is the carrier density. In these equations I(t)is the injected current, Va the active volume, e the electron charge, vg the group velocity, g(n)the material gain, ∊ the non-linear gain coefficient, Γ the optical confinement factor, τ p the photon lifetime, Rsp(n)is the rate of the spontaneous emission coupled into the lasing mode, and α is the linewidth enhancement factor. The carrier recombination rate is an + bn2+cn3, where a,b and c are the non-radiative, spontaneous, and Auger recombination coefficients, respectively. We consider a material gain g(n)given by g(n) = dg dn(n−nt), where dg dn is the differential gain and nt the transparency carrier density. The Langevin terms Fp(t)and Fϕ(t)in Eqs. (1)–(2), represent fluctuations due to spontaneous emission, with the following correlation properties, <Fi(t)Fj(t′)>=δijδ(t−t′), where δ(t)is the Dirac delta function and δij the Kronecker delta function with the subindexes i and j referring to the variables p and ϕ. In order to write the corresponding equations for the number of photons inside the laser, P(t), and the number of carriers in the active region, N(t), we make the following change of variables: P=pVp,N= nVa, where Vp is the volume occupied by the photons, to obtain: dP dt =[GN(N−Nt) 1+∊P−1 τ p]P+Rsp(N)+ 2Rsp(N)P √Fp(t)(4) dϕ dt = α 2[GN(N−Nt)−1 τ p]+ Rsp(N) 2P √Fϕ(t)(5) dN dt =I(t) e−(AN +BN2+CN3)−GN(N−Nt)P 1+∊P(6) where GN=Γvgdg dn/Va,Nt=ntVa,A=a,B=b/Va,C=c/V2 a,∊=∊Γ/Va, and where we have used that Γ=Va/Vp. Our departure equations have been written for the photon density (or number), and optical phase. This permits a straightforward comparison with the models that analyze QRNG systems based on gain-switching of single-mode laser diodes [19,22,24,25] The above mentioned systems are such that a large signal modulation of the bias current is applied to the device in such a way that a random evolution of the optical phase is induced by the spontaneous emission noise, specially when the bias current is below the threshold value. Theoretical information on these generators has been mainly obtained by numerical integration of the stochastic rate equations. Integration of these equations is usually performed by using the Euler–Maruyama method for numerical integration of stochastic differential equations [31,32]. However some problems can appear when integrating Eqs. (4)– (6). If the photon number is very small (that is for bias currents close or below threshold) negative photon numbers can be obtained that lead to numerical instabilities due to the square root factor multipliying the noise term in Eqs. (4)–(5). A possible solution could be solving the equations with the constraint that the photon number is non-negative, and so P is equalled to zero when it becomes negative [25]. However, this algorithm would not work because a division by zero would be done in Eq. (5) in order to calculate the phase value at the following integration step. This problem also appears in all numerical algorithms in which the same operation is performed, for instance the Heun’s predictor–corrector algorithm[27]. In any case, numerical solutions of the equations without any type of constraint would be desirable. A different model, in which the variables P and N appearing in the square root terms multiplying the random terms are substituted by their averaged values, has been considered [25]. However the statistics of the variables can be A. Quirce and A. Valle Optics and Laser Technology 150 (2022) 107992 3 different to those obtained with Eqs. (4)–(6) because the stochastic differential equations are different. In order to get an integration algorithm for Eqs. (4)–(6) without the previous problems we begin by observing that the frequency of events in which a negative P appears decreases as the integration time step is decreased. In this way a solution is to decrease that time step until no instabilities are observed. However that situation occurs when the time step is so small that integration for a large number of periods in order to get significant statistics is not practical. The usual solution for this problem has been the integration of the corresponding rate equations for the complex electric field, E, instead of equations for P and ϕ [33]. P does not appear inside the square root factors that multiply the noise terms and hence no instabilities are observed. However, when changing the variables from usual rate equations for E [33] to P and ϕ we obtain that a gain saturation term appears in the equation for ϕ, in contrast to Eq. (5). Taking into account that a non-saturated linear gain term must appear in the equation for the phase [25,27,29] we need to modify the electric field rate equations of [33] in order to be consistent with Eq. (5). We then consider the following rate equations: dE dt =[( 1 1+∊E|2+i α )GN(N−Nt)−1+i α τ p]E 2+ Rsp(N) √ξ(t)(7) dN dt =I(t) e−(AN +BN2+CN3)− GN(N−Nt)E|2 1+∊E|2(8) where E(t) = E1(t) +iE2(t)is the complex electric field and ξ(t) = ξ1(t) +iξ2(t)is the complex Gaussian white noise with zero average and correlation given by <ξ(t)ξ*(t′)>=δ(t−t′)that represents the spontaneous emission noise. The only change with respect to Eq. (7) of [33] is that the saturation term does not multiply i α in our first term of Eq. (7). Applying the rules for the change of variables in the Ito’s calculus[34] to P=E|2=E2 1+E2 2 and ϕ=arctan(E2/E1)in Eqs. (7)–(8), we obtain our initial Eqs. (4)–(6). Therefore the integration of the electric field rate Eqs. (7)–(8) solve the previously mentioned problems because they correspond to our departure equations and they avoid numerical instabilities. In the next section we will describe experiments performed with the single longitudinal mode semiconductor laser that we used in [28]. In [27] measurement of the intrinsic parameters of a similar laser was performed. We will consider most of these parameters for our laser. The numerical values of the parameters are then: GN=1.48 ×104s −1 , Nt= 1.93 ×107,∊=7.73 ×10−8, τ p=2.17ps, α =3,A=2.8×108 s −1 , B= 9.8 s −1 , and C=3.84 ×10−7 s −1 . Simulation and experimental results show not only qualitative but also a remarkable quantitative agreement for a very wide range of gain-switching conditions [27,28]. 3. Experimental results A Discrete Mode Laser (DML) (Eblana Photonics EP1550-0-DM-H19FM) was used in our experiments. The device is a multi-quantum well laser in a ridge waveguide with index perturbations for inducing singlemode operation. The temperature of the device was held constant at 25◦ C by using a temperature controller (Luzwavelabs LDC/E-Temp3). The CW bias current, I, was controlled with the Luzwavelabs LDC/ECurrent200 laser driver. The threshold current of the laser is Ith= 14.14 mA, the emission wavelength is 1546.985 nm at I =30 mA, and the wavelength separation between consecutive longitudinal modes is 1.28 nm. In order to check that the parameters measured in [27] are appropriated for our device we have performed the measurements shown in Fig. 1. In Fig. 1(a) the square of the relaxation oscillation frequency, fr, obtained from relative intensity noise (RIN) spectra as a function of I−Ith, is shown. A linear relation is obtained. The slope of the linear fit, (2.35 ±0.04) GHz 2 /mA, gives a value for the differential gain GN= (1.48 ±0.03) × 104 s −1 [35], very similar to that in [27], 1.47 ×104s −1 . The value for Nt that we consider is slightly different to that in [27] because the threshold currents are slightly different. In order to check the values of A,B, and C we measure the laser linewidth as a function of |I−Ith|for bias currents slightly below threshold with a high resolution (10 MHz) optical spectrum analyzer, BOSA (Aragon Photonics, BOSA 210). We follow the measurement procedure of [35] in order to obtain the width (FWHM) of the optical spectrum, Δ ν . Our experimental results for Δ ν are shown in Fig. 1(b) indicating a linear dependence with the current. This dependence is also obtained from theoretical analysis because when the current is slightly below threshold it is shown that [35]: Δ ν = τ nGN 2 π e|I−Ith|(9) where τ n is the differential carrier lifetime at threshold. Using the slope of the corresponding linear fit, (3.71 ±0.12) GHz/mA, together with Eq. (9) we obtain a value of τ n= (0.25 ±0.01)ns. We can compare with the result obtained from our parameters because τ n is given by τ n= (A+2BNth +3CN2 th)−1 where Nth =Nt+1/(GN τ p) = 5.045 ×107 is the carrier number at threshold [26]. Using our parameters in this expression we obtain τ n=0.24 ns. This value is consistent with that derived from Fig. 1(b) validating in this way our choice of parameters. We note that A,B, and C parameters were extracted in [27] using a different method, mainly based on the measurement of RIN spectra. The Schawlow-Townes law gives the value of the linewidth below threshold and is accurate for the amplified spontaneous emission (ASE) regime [26]. This formula reads Δ ν =Rsp(N) 2 π P(10) where P and N are the averaged steady-state values of the photon and carrier number, respectively. The approximation of Rsp by its value at threshold [35] gives precisely Eq. (9). If instead doing this Fig. 1. (a) Squared relaxation oscillation frequencies as a function of I−I th. (b) Linewidth as a function of |I−Ith|for bias currents below threshold. Linear fits are also plotted with solid lines. A. Quirce and A. Valle Optics and Laser Technology 150 (2022) 107992 4 approximation we use the general Eq. (10) for an extended range of currents below threshold we can measure Δ ν in order to get information about Rsp(N). This measurement is done with an optical spectrum analyzer (OSA, Anritsu MS9710B) with a worse wavelength resolution (0.06 nm) but with a better sensitivity than those of BOSA. The linewidth and the power emitted close to the lasing mode frequency, Pw, as a function of the bias current are shown in Fig. 2. Pw has been measured by centering the ASE spectrum around its maximum and using a wavelength span of 1 nm. Δ ν (Pw) decrease (increase) in the well-below threshold regime as the bias current increases as illustrated in Fig. 2. Their product can be used to evaluate Rsp(N) because Rsp(N) = 2 π Δ ν P and P is proportional to Pw, that is P=kPw where k is the proportionality constant. Fig. 3(a) shows that the product Δ ν Pw is an increasing function of I/Ith. In order to know the dependence of Δ ν Pw as a function of the carrier number we can use Eq. (6) to calculate the steady state value of N corresponding to I. In the well-below threshold regime the stimulated emission term can be neglected and N is obtained by solving the cubic equation: AN+BN2+CN3−I/e=0, with the A,B,C corresponding to our laser. Fig. 3(b) shows the real root of this equation as a function of the current. Using this dependence in the experimental data shown in Fig. 3(a) we can obtain Δ ν Pw as a function of N. This dependence is shown in Fig. 3 (c) and would give Rsp =Rsp(N), provided k is known, because Rsp(N) = 2 π kΔ ν Pw(N). A linear dependence of Rsp on N is typically used in rate equation models. A more rigorous treatment [26] considers a quadratic dependence: Rsp(N) = βBN2 where β is the fraction of spontaneous emission coupled into the lasing mode. Assuming this dependence we get Δ ν Pw(N) = βBN2/(2 π k), so we can get the value of k if Δ ν Pw is linearly fit vs N2 and if the value of β measured for the laser is used (β=5.3× 10−6) [27]. The linear fit, shown in Fig. 4(a), has a slope s= (7.36 ± 0.04) × 10−13 HzW, and so k=βB/(2 π s) = (1.12 ±0.01) × 107 W −1 . We now analyze the validity of the linear approximation for Rsp(N). The estimation of k permits us to calculate Rsp(N) = 2 π kΔ ν Pw(N). Values of Rsp(N)for our experimental data are shown in Fig. 4(b). Also two linear fits are plotted in that Figure. The first one considers a dependence Rsp(N) = β1(N−N0)while the second one is a linear fit in which crossing the origin is imposed, that is Rsp(N) = β2N. The values of the fit parameters are: β1= (4900 ±200)s −1 , N0= (2.39 ±0.17) × 107, and β2= (2300 ±100)s −1 . The β2 parameter has been calculated as the slope of the regression straight line obtained when fixing the intercept of that line to the zero point, (N,Rsp) = (0,0). We note that the dependence Rsp(N) = βN, where β is constant, is frequently used in the literature [33,36,37] being usually considered when modelling QRNG in gain Fig. 2. (a) Full width at half maximum of the optical spectrum and (b) power as a function of the bias current relative to threshold. Fig. 3. (a) Linewidth-power product as a function of I/I th (b) Calculated carrier number relative to threshold vs I/Ith, and (c) Linewidth-power product as a function of N/Nth. Fig. 4. (a) Linewidth-power product as a function of N2 (b) Rate of spontaneous emission coupled into the lasing mode as a function of N. Results for linear fits are plotted with red and blue solid lines. Theoretical results obtained for the quadratic dependence are shown with black solid line. A. Quirce and A. Valle Optics and Laser Technology 150 (2022) 107992 5 switched laser diodes [6,19,22,24,25]. Results from Fig. 4(b) show that while Rsp(N) = β1(N−N0)seems a reasonable approximation, considering Rsp(N) = β2N does not seem in principle adequate in order to get a proper quantitative agreement between experimental and theoretical results. The better description obtained when using Rsp(N) = β1(N−N0)can be understood if the Taylor expansion of N2 around Nth is performed: N2∼N2 th +2Nth(N−Nth). In this way Rsp(N)= βBN2∼βB(2NthN−N2 th)=2βBNth(N−Nth 2). The parameters of the linear dependence are then β1=2βBNth and N0=Nth/2. The evaluation of these parameters using those of the model gives β1=5240s −1 and N0=2.52 ×107, values that are close to those obtained with the linear fit with two parameters. Fig. 4(b) also shows the values obtained with Rsp(N) = βBN2. As expected, the linear fit with two parameters becomes a good approximation of the quadratic rate of spontaneous emission as N approaches Nth. In the next section we will discuss the validity of the linear approximations for describing phase diffusion. 4. Theoretical results In this section we will theoretically analyze the dynamical evolution of the statistics of the relevant variables with a special focus on the optical phase. We consider an injected bias current that follows a squarewave modulation of period T with I(t) = Ion during T/2, and I(t) = Ioff during the rest of the period. Fig. 5(a), (b), and (c) show the dynamics of the photon number, optical phase and carrier number, respectively, when Ion =30mA, Ioff =7 mA, and T=1 ns for three consecutive periods. The laser is switched-off with a bias current well below threshold, close to Ith/2, in order to get randomness of the phase due to the spontaneous emission noise. Integration of Eqs. (7)–(8) with Rsp(N) = βBN2 is performed by using the Euler numerical scheme [31,32] with an integration time step of 0.01 ps. We integrate these 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. All our results will be plotted in one time window of duration T. Fig. 5(a) shows P in a logarithmic scale in order to better appreciate the fluctuations. At first P is random and determined by the spontaneous emission events. After the switch-on, P develops relaxation oscillations and begins to decrease at t=0.5 ns that is when Ioff is applied. Around t=0.7 ns, P reach again random values with an average similar to that found at the beginning of the period. The optical phase for those three realizations is shown in Fig. 5(b). ϕ is calculated at each integration step from E1 and E2 in such a way that it is a continuous function of time within each period. Fig. 5(b) shows that ϕ decreases after the bias current changes to Ioff , in such a way that the final value at the end of the period is much smaller than 0 (more than several multiples of 2 π ). The initial value of the phase at the beginning of the period, ϕ0, is taken as the value of the phase at the end of the previous period, ϕT, but in the [0, 2 π )range, that is ϕ0=ϕT– int(ϕT 2 π )2 π . Note that the conversion to the [0,2 π )range is necessary if a calculation of well defined statistical moments of the phase is required. If no conversion is done, not even the averaged phase, <ϕ(t)>, could be calculated because ϕdecreases in each period in a magnitude of more than several 2 π . For large values of P the noise term in Eq. (5) is negligible, so ϕ evolves in a deterministic way. However, for small values of P the noise term in Eq. (5) dominates the other term of that equation and so ϕ mainly evolves in a random way. Therefore, similarly to P, fluctuations of ϕ are more important at the beginning and at the end of the pulse. The carrier number evolution is shown in Fig. 5(c). N displays relaxation oscillations before beginning to decrease after t=0.5 ns. This decrease is monotonous since Ioff is smaller than the threshold current. The results included in Fig. 5 have been obtained by integrating Eqs. (7)–(8). We have also integrated the equations for the photon number and optical phase, Eqs. (4)–(6), with different integration time steps, Δt, under the conditions of Fig. 5. In these integrations we have recorded the percentage of integration time steps in which P becomes negative. These percentages are 2.3, 0.9 and 0.1% for Δt=0.1, 0.01, and 0.001 ps, respectively. These results indicate that even with a very small Δt, like 0.001 ps, problems are found in a significant fraction of integration steps when using Eqs. (4)–(6). The quantity that is relevant for determining the performance of this system as a QRNG is the standard deviation of the phase. This is shown as a function of time in Fig. 5(d). σ ϕ(t)has been calculated by averaging over 5 ×104 periods. σ ϕ(0)>0 due to our choice of random initial conditions. Large increases of σ ϕ(t), corresponding to phase diffusion, occur at the beginning and at the end of the period because P has small values determined by the spontaneous emission noise. Phase diffusion is characterized by a linear increase of σ 2 ϕ with t. We can check this dependence by fitting the values of σ 2 ϕ from t =0.75 ns to t =1 ns: we obtain a linear fit with a correlation coefficient of 0.999. During the deterministic evolution of ϕ (between 0.2 and 0.6 ns, approximately) σ 2 ϕ oscillates around a value that increases linearly with t [38]. The frequency of these oscillations is the relaxation oscillation frequency[38]. These oscillations can be barely seen in Fig. 5(b) because their amplitude is small compared to the value of σ 2 ϕ(T). From now on we will evaluate the suitability of the different linear spontaneous emission rate terms considered in the previous section for Fig. 5. (a) Photon number, (b) phase, and (c) carrier number evolution for three consecutive periods are shown with black, red, and green lines. (d) Standard deviation of the phase as a function of t. In this figure I off =7 mA and T =1 ns. A. Quirce and A. Valle Optics and Laser Technology 150 (2022) 107992 6 describing the transient statistics of our system. We will compare results obtained with Rsp,1=β1(N−N0), Rsp,2=β2N, and Rsp,3=β1N with the results obtained with Rsp =βBN2, that is the best available description of that rate. The advantage of Rsp,1 and Rsp,3 with respect to Rsp,2 is that β1 and N0 have simple expressions in terms of the parameters of the model (β1=2βBNth,N0=Nth/2) in contrast to β2 that must be calculated from a fitting of the experimental results. We have considered Rsp,3 because it corresponds to the linear dependence of Rsp on N for which the parameter, β1, has a simple expression in terms of the parameters of the model. Figs. 6(a), 6(b) and 6(c) show the dynamical evolution of the averaged photon number, phase and carrier number, respectively, considering the previously mentioned dependences for Rsp(N), when Ioff has been increased to 12 mA. This bias current is below the threshold value and is enough for making the phase to diffuse at the beginning and the end of the period as it can be seen in Fig. 6(d) in which σ ϕ(t)has been plotted. Fig. 6 shows that the linear approximations are good for describing the statistics with the exception of Rsp,3. It would be expected that the best description is that corresponding to Rsp,1 because it corresponds to the linear expansion of βBN2 around Nth. Ioff is close to threshold in Fig. 6 (0.85 Ith), so the deviations of N with respect to Nth are not large. We can confirm that it is the best approximation by making zooms of Fig. 6(b) and Fig. 6(d) that are shown in the insets of these figures. The inset of Fig. 6(d) also shows in a much clearer way the oscillations of σ ϕ(t)at the relaxation oscillation frequency. The reason why Rsp,3 gives the worst description is that it overestimates the spontaneous emission rate as it can be seen in Fig. 4(b): Rsp,3 corresponds to a straight line with β1 slope and crossing through the origin. This overestimation can also be seen in Fig. 6(a) in which <P>obtained with Rsp,3 has the largest values at the beginning and the end of the period. The increase of the speed of QRNGs based on gain-switched lasers is obtained by decreasing the modulation period. For very short periods Ioff can not be very small because then laser pulses with significant power would not be fired. In this sense it is interesting to check again the validity of the linear approximations when T is decreased. Fig. 7(a) shows σ ϕ(t)under the same conditions of Fig. 6 but decreasing T to 0.3 ns. The best linear approximation is that corresponding to Rsp,1. The relative errors of σ ϕ(T)for Rsp,1,Rsp,2 and Rsp,3 are 0.2%, 4% and 19%, respectively. These relative errors are smaller when Ioff decreases in order to increase σ ϕ as it is shown in Fig. 7(b) where results for T=0.3 ns and Ioff =7mA. They are 0.15%, 0.9% and 6%, for Rsp,1,Rsp,2 and Rsp,3, respectively. In order to gain insight on the gain-switched operation we have done simulations under various operating conditions. Fig. 8 shows σ ϕ(T)as a function of Ioff ,Ion, and T, taking as a reference the situation illustrated in Fig. 5. In this way we analyze the dependence of the fluctuations of the phase on bias current, pulse current amplitude and pulse frequency. Fig. 8(a) shows that σ ϕ(T)increases monotonously as Ioff decreases since the strength of the phase fluctuations is more important at the beginning and the end of the period (this is also illustrated when comparing Fig. 5 and Fig. 6). The situation is not monotonous when analyzing σ ϕ(T)as a function of Ion because there are some oscillations, as it can be seen in Fig. 8(b). The origin of these oscillations is related to the relaxation oscillations. Our gain saturation coefficient is such that only two clear relaxation oscillation peaks appear before reaching the steady state. The maximum value of σ ϕ(T)at Ion=19 mA corresponds to a situation in which just one peak is emitted with its maximum appearing at T/2. The minimum value is obtained at Ion=23 mA. It corresponds to the situation in which the second relaxation oscillation spike is beginning to appear. The maximum at Ion=26 mA corresponds to the situation in which there is no emission after the second fully developed relaxation oscillation peak. The last minimum, at Ion=30 mA corresponds to the situation at which the photon number at T/2 is close to its steady state value, as it can be seen in Fig. 5. Fig. 8(c) shows the variation of σ ϕ(T)as a function Fig. 6. (a) Averaged photon number, (b) averaged phase, (c) averaged carrier number, and (d) standard deviation of the phase as a function of t. Results for different dependences of R sp(N)are plotted with lines of different colours, as indicated in the legend. In this figure Ioff =12 mA and T =1 ns. Fig. 7. Standard deviation of the phase as a function of time for T=0.3 ns and (a) I off =12 mA, and (b) Ioff =7 mA. Results for different dependences of Rsp(N) are plotted with lines of different colours, as indicated in the legend. A. Quirce and A. Valle Optics and Laser Technology 150 (2022) 107992 7 of T. For small values of T there are some oscillations that can be attributed to the relaxation oscillations of the laser. When T>0.65ns the increase of σ ϕ(T)is monotonous: longer values of the time for which the laser is switched off produce larger values of the fluctuations of the phase. 5. Discussion and conclusions The best description for the spontaneous emission rate is that corresponding to a quadratic dependence on the carrier number[26]. The considered linear approximations have advantages and disadvantages. Rsp,1 is the best approximation when the carrier number is close to Nth. Simulation results obtained with Rsp and Rsp,1 are similar when N is close to Nth because Rsp,1 is an approximation to Rsp that only works well close to Nth because it is derived from a Taylor expansion of Rsp(N)around Nth. However, the use of Rsp,1 in simulations gives unphysical results when the bias current is small enough to obtain N smaller than Nth/2 because we would obtain a negative spontaneous emission rate since Rsp,1= 2βBNth(N−Nth/2). This also produces numerical problems when integrating Eqs. (4) and (5) since Rsp,1 appears within a square root sign. That problem does not appear for the other considered linear dependences. However, Rsp,3=2βBNthN has the problem of a very large overestimation of the spontaneous emission rate that causes very large values of <P>when the spontaneous emission noise dominates the evolution (see Fig. 6). The use of Rsp,2=β2N has a different problem: β2 has no expression in terms of the parameters of the model and so, no direct physical meaning can be assigned to β2. Also β2 must be obtained from an extra fitting, that shown in Fig. 4. Concerning the description of the phase fluctuations, Fig. 6, and Fig. 7 show that the best (worst) linear description is given by Rsp,1 (Rsp,3). We have discussed three different expressions for the linear dependence of the spontaneous emission rate because they are used in the literature when modelling gain switching in laser diodes. Comparison of these dependences with the quadratic one is of interest for knowing the validity of the linear approximations. The simulation results relative to the phase statistics, mainly average and standard deviation of the phase, have not been compared with experiments, so direct validation of these results has not been done in this work. This will be the subject of future work. Our approach for obtaining a good description of the phase evolution has been to compare experimental and theoretical results in the static case for extracting the model parameters. Using these parameters in the rate equation model we have obtained the time evolution of the phase statistics. Most of the results shown in our work correspond to long current pulses in such a way that more than one relaxation oscillations peak is observed. We have chosen these conditions because they correspond to a typical operation of commercial QRNGs based on gain-switching of laser diodes. Summarizing, we have analyzed theoretically and experimentally the phase diffusion appearing in gain-switched single-mode semiconductor lasers. This study finds application in multi-photon QRNGs based on interferometric detection of pulses with random phases because it contributes to a better statistical characterization of the entropy source of these generators. We have obtained an electric field stochastic rate equations model with the advantage of avoiding the numerical instabilities that appear for the below threshold operation, that is the regime in which the randomization is achieved. We have used [27] for an initial estimation of the laser parameter values. Further measurements have been performed in order to get a complete set of parameters of the model for our laser. Using the Schawlow-Townes formula, that gives a good description of the below threshold situation, we have experimentally obtained the spontaneous emission rate coupled into the lasing mode as a function of the carrier number. We have compared the obtained quadratic dependence with the linear approximations that have been used to describe laser phase diffusion in QRNG, discussing advantages and disadvantages of the linear approach. We have calculated the dynamical evolution of the statistics of the optical phase. A realistic quantitative description of the variance of the optical phase has been obtained because it relies in a model for which the extraction of the laser parameters has been performed. CRediT authorship contribution statement Ana Quirce: Methodology, Software, Validation, Investigation, Visualization. Angel Valle: Conceptualization, Methodology, Software, Formal analysis, Investigation, Writing – original draft, Writing – review & editing, Funding acquisition. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgment This research was supported by the Ministerio de Ciencia e Innovaci´ on, RTI2018-094118-B-C22 MCIN/AEI/10.13039/ 501100011033/ FEDER “Una manera de hacer Europa”. A. Quirce acknowledges financial support from Beatriz Galindo program, Ministerio de Ciencia, Innovaci´ on y Universidades (Spain). Angel Valle would like to thank Marcos Valle Mi˜ n´ on for his help in the calculations. References [1] M. Stipˇ cevi´ c, Ç.K. Koç, True random number generators, in: Open Problems in Mathematics and Computational Science, Springer, 2014, pp. 275–315. [2] M. Herrero-Collantes, J.C. Garcia-Escartin, Quantum random number generators, Rev. Mod. Phys. 89 (1) (2017) 015004. [3] M. Sciamanna, K.A. Shore, Physics and applications of laser diode chaos, Nature Photon. 9 (3) (2015) 151–162. [4] P. Li, Y. Guo, Y. Guo, Y. Fan, X. Guo, X. Liu, K. Li, K.A. Shore, Y. Wang, A. Wang, Ultrafast fully photonic random bit generator, J. Lightwave Technol. 36 (12) (2018) 2531–2540. Fig. 8. Standard deviation of the phase at the end of the period as a function of (a) I off (for Ion =30 mA, T=1 ns), (b) Ion(for Ioff =7 mA, T=1 ns), and (c) T (for Ion =30 mA, Ioff =7 mA). A. Quirce and A. Valle Optics and Laser Technology 150 (2022) 107992 8 [5] P. Li, K. Li, X. Guo, Y. Guo, Y. Liu, B. Xu, A. Bogris, K.A. Shore, Y. Wang, Parallel optical random bit generator, Opt. Lett. 44 (10) (2019) 2446–2449. [6] C. Abellan, W. Amaya, D. Domenech, P. Mu˜ noz, J. Capmany, S. Longhi, M. W. Mitchell, V. Pruneri, Quantum entropy source on an inp photonic integrated circuit for random number generation, Optica 3 (9) (2016) 989–994. [7] M. Stipˇ cevi´ c, B.M. Rogina, Quantum random number generator based on photonic emission in semiconductors, Rev. Sci. Instrum. 78 (4) (2007) 045104. [8] T. Jennewein, U. Achleitner, G. Weihs, H. Weinfurter, A. Zeilinger, A fast and compact quantum random number generator, Rev. Sci. Instrum. 71 (4) (2000) 1675–1680. [9] W. Wei, H. Guo, Bias-free true random-number generator, Opt. Lett.s 34 (12) (2009) 1876–1878. [10] H. Fürst, H. Weier, S. Nauerth, D.G. Marangon, C. Kurtsiefer, H. Weinfurter, High speed optical quantum random number generation, Opt. Exp. 18 (12) (2010) 13029–13037. [11] T. Durt, C. Belmonte, L.-P. Lamoureux, K. Panajotov, F. Van den Berghe, H. Thienpont, Fast quantum-optical random-number generators, Phys. Rev. A 87 (2) (2013) 022339. [12] A. Argyris, E. Pikasis, S. Deligiannidis, D. Syvridis, Sub-tb/s physical random bit generators based on direct detection of amplified spontaneous emission signals, J. Lightwave Technol. 30 (9) (2012) 1329–1334. [13] Y. Guo, Q. Cai, P. Li, Z. Jia, B. Xu, Q. Zhang, Y. Zhang, R. Zhang, Z. Gao, K.A. Shore, et al., 40 gb/s quantum random number generation based on optically sampled amplified spontaneous emission, APL Photon. 6 (6) (2021) 066105. [14] Y. Shen, L. Tian, H. Zou, Practical quantum random number generator based on measuring the shot noise of vacuum states, Phys. Rev. A 81 (6) (2010) 063814. [15] H. Guo, W. Tang, Y. Liu, W. Wei, Truly random number generation based on measurement of phase noise of a laser, Phys. Rev. E 81 (5) (2010) 051137. [16] B. Qi, Y.-M. Chi, H.-K. Lo, L. Qian, High-speed quantum random number generation by measuring phase noise of a single-mode laser, Opt. Lett. 35 (3) (2010) 312–314. [17] F. Xu, B. Qi, X. Ma, H. Xu, H. Zheng, H.-K. Lo, Ultrafast quantum random number generation based on quantum phase fluctuations, Opt. Exp. 20 (11) (2012) 12366–12377. [18] M. Jofre, M. Curty, F. Steinlechner, G. Anzolin, J. Torres, M. Mitchell, V. Pruneri, True random numbers from amplified quantum vacuum, Opt. Exp. 19 (21) (2011) 20665–20672. [19] C. Abell´ an, W. Amaya, M. Jofre, M. Curty, A. Acín, J. Capmany, V. Pruneri, M. Mitchell, Ultra-fast quantum randomness generation by accelerated phase diffusion in a pulsed laser diode, Opt. Exp. 22 (2) (2014) 1645–1654. [20] Z. Yuan, M. Lucamarini, J. Dynes, B. Fr¨ ohlich, A. Plews, A. Shields, Robust random number generation using steady-state emission of gain-switched laser diodes, Appl. Phys. Lett. 104 (26) (2014) 261112. [21] D.G. Marangon, A. Plews, M. Lucamarini, J.F. Dynes, A.W. Sharpe, Z. Yuan, A. J. Shields, Long-term test of a fast and compact quantum random number generator, J. Lightwave Technol. 36 (17) (2018) 3778–3784. [22] R. Shakhovoy, V. Sharoglazova, A. Udaltsov, A. Duplinskiy, V. Kurochkin, Y. Kurochkin, Influence of chirp, jitter, and relaxation oscillations on probabilistic properties of laser pulse interference, IEEE J. Quantum Electron. 57 (2) (2021) 1–7. [23] K. Nakata, A. Tomita, M. Fujiwara, K.-I. Yoshino, A. Tajima, A. Okamoto, K. Ogawa, Intensity fluctuation of a gain-switched semiconductor laser for quantum key distribution systems, Opt. Exp. 25 (2) (2017) 622–634. [24] B. Septriani, O. de Vries, F. Steinlechner, M. Gr¨ afe, Parametric study of the phase diffusion process in a gain-switched semiconductor laser for randomness assessment in quantum random number generator, AIP Adv. 10 (5) (2020) 055022. [25] R. Shakhovoy, A. Tumachek, N. Andronova, Y. Mironov, Y. Kurochkin, Phase randomness in a gain-switched semiconductor laser: stochastic differential equation analysis, arXiv preprint arXiv:2011.10401. [26] L.A. Coldren, S.W. Corzine, M.L. Mashanovitch, Diode Lasers and Photonic Integrated Circuits, vol. 218, John Wiley & Sons, 2012. [27] A. Rosado, A. P´ erez-Serrano, J.M.G. Tijero, A. Valle, L. Pesquera, I. Esquivias, Numerical and experimental analysis of optical frequency comb generation in gainswitched semiconductor lasers, IEEE J. Quantum Electron. 55 (6) (2019) 1–12. [28] A. Quirce, A. Rosado, J. Díez, A. Valle, A. P´ erez-Serrano, J.-M.G. Tijero, L. Pesquera, I. Esquivias, Nonlinear dynamics induced by optical injection in optical frequency combs generated by gain-switching of laser diodes, IEEE Photonics J. 12 (4) (2020) 1–14. [29] N. Schunk, K. Petermann, Noise analysis of injection-locked semiconductor injection lasers, IEEE J. Quantum Electron. 22 (5) (1986) 642–650. [30] K. Petermann, Laser Diode Modulation and Noise, vol. 3, Springer Science & Business Media, 1991. [31] H. Risken, Fokker-planck equation, in: The Fokker-Planck Equation, Springer, 1996. [32] P.E. Kloeden, E. Platen, Stochastic differential equations, in: Numerical Solution of Stochastic Differential Equations, Springer, 1992. [33] S. Balle, N.B. Abraham, P. Colet, M. San Miguel, Parametric dependence of stochastic frequency variations in the gain switching of a single-mode laser diode, IEEE J. Quantum Electron. 29 (1) (1993) 33–41. [34] C.W. Gardiner, et al., Handbook of Stochastic Methods, vol. 3, Springer, Berlin, 1985. [35] P. P´ erez, A. Valle, I. Noriega, L. Pesquera, Measurement of the intrinsic parameters of single-mode vcsels, J. Lightwave Technol. 32 (8) (2014) 1601–1607. [36] G.P. Agrawal, N.K. Dutta, Semiconductor Lasers, Springer Science & Business Media, 2013. [37] C. Mirasso, A. Valle, L. Pesquera, P. Colet, Simple method for estimating the memory diagram in single mode semiconductor lasers, IEE Proc.-Optoelectron. 141 (2) (1994) 109–113. [38] C. Henry, Phase noise in semiconductor lasers, J. Lightwave Technol. 4 (3) (1986) 298–311. A. Quirce and A. Valle