Kyungduk Kim

Stefan Bittner

Yongquan Zeng

Stefano Guazzotti

Qi Ortwin Hess

Jie Wang

Hui Cao
email: hui.cao@yale.edu

Massively parallel ultrafast random bit generation with a chip-scale laser

come   L'archive ouverte pluridisciplinaire

tions at different spatial locations impede independent parallel bit stream generation. Spectral demultiplexing of amplified spontaneous emission (22,23) or heterodyning chaotic laser emission (9) are employed for parallel RBG with rates up to hundreds of Gb/s per channel. So far, such spectral-domain parallel RBG has been demonstrated with less than 10 channels.

We demonstrate a method that enhances the random bit rate in a single channel and also provides hundreds of channels for simultaneous generation of independent bit streams. The spatio-temporal interference of many lasing modes is used to generate picosecond-scale emission intensity fluctuations in space, so as to massively produce ultrafast random bit streams in parallel. This is achieved by tailoring the geometry of a broad-area semiconductor laser to vastly increase the number of transverse lasing modes, suppressing characteristic dynamical instabilities such as filamentation. Specifically, we design a chip-scale laser diode to enable a large number of spatial modes lasing simultaneously with incommensurate frequency spacings, so that their interference patterns are complex and aperiodic. Spontaneous emission adds stochastic noise to make the intensity fluctuations unpredictable and non-reproducible.

A conventional broad-area edge-emitting semiconductor laser has a stripe geometry with two flat facets (Fig. 1A). Characteristically, lasing occurs only in the low-order transverse modes. Nonlinear interactions between the light field and the gain material entail irregular pulsation and filamentation (24) (Fig. 1B). The spatio-temporal correlation function of the intensity fluctuations (25,26) C(∆x, ∆t) reveals non-local correlations in space and time (Fig. 1C).

On one hand, long-range temporal correlation reflects memory which degrades the quality of random bits generated at one spatial location. On the other hand, due to long-range spatial correlation, the random bit streams generated at different locations are not completely independent, thus impeding parallel RBG (21).

To achieve massively parallel ultrafast RBG, we enhance the number of transverse lasing modes by increasing the cavity width and curving the end facets (Fig. 1D), effectively sup-pressing modulational instabilities. High-order transverse modes are well confined inside such a cavity, and their optical gain is enhanced by tailoring the top metal contact shape (26). The number of transverse lasing modes is maximized by fine-tuning the cavity geometry (27). Lasing on small length-scales of transverse-wavelengths of high-order modes prevents lensing and self-focusing effects (28) that would normally cause filamentation and instabilities (Fig. 1E). In turn, the absence of filaments and pulsations eliminates long-range spatio-temporal correlations in the lasing intensity (Fig. 1F). It is this dramatically shortening of the correlation lengths in space and time, that paves the ground for a significant increase in the number of independent spatial channels for parallel RBG, and a great enhancement of the RBG rate of every individual spatial channel.

With lasing instabilities suppressed, the dynamic variations of the emission intensity are orchestrated by the interference of lasing modes with different frequencies. The characteristic time scale of such intensity fluctuations is inversely proportional to the spectral width of the total emission, and is about 1 ps for the GaAs quantum-well laser (26). We show the spatio-temporal beat pattern of the intensity emitted at one laser facet in Fig. 2A. The temporal correlation length is determined by the full-width-at-half-maximum (FWHM) of C(∆x, ∆t) in time (Fig. 2B). Its value of 2.8 ps is limited by the temporal resolution of our detection (26).

Thanks to the ultrafast dynamics of lasing intensity, the power spectrum is extremely broad (Fig. 2C). Its bandwidth, which contains 80% of the radio-frequency (RF) power, is 315 GHz.

For comparison, the numerically-simulated spectrum is even broader with a bandwidth of 632 GHz (26). After accounting for the temporal resolution of photodetection, the simulated RF spectrum matches the measured one (Fig. 2C). This agreement confirms that the ultra-broad power spectrum results from the interference of many transverse and longitudinal modes.

For spatial-multiplexing of RBG, the number of independent parallel channels depends on the spatial correlation length of the lasing emission. Now with non-local correlations removed, the local correlation length estimated from the spatial FWHM of C(∆x, ∆t) is 1.5 µm (Fig. 2B), which is limited by the spatial resolution of our detection. Without the finite experimental resolution, our simulation gives a correlation length of 0.5 µm (26). Thanks to this extremely short spatial correlation length, hundreds of independent spatial channels are available for parallel RBG.

As the transverse mode frequency spacing in our cavity design is incommensurate to the longitudinal mode spacing, the spatio-temporal interference pattern cannot repeat itself (26). Moreover, the spontaneous emission, generated by quantum fluctuations, constantly feeds stochastic noise into the lasing modes, making their beat pattern unpredictable and irreproducible.

To generate random bits, we divide the laser end facet into 1-µm-wide spatial channels. Due to the restricted field-of-view of our collection optics, only 254 spatial channels are recorded simultaneously, which is about half of the number possible with complete collection of emission. We sample the emission intensity at every spatial channel in τ = 1.46 ps long intervals (sampling rate 683 GHz, left inset of Fig. 3A) (26). The emission intensity integrated over one sampling period I n has an asymmetric probability density function (PDF), which would yield biased bits (26). We adopt the procedure from Ref. (4) to acquire the differential intensity ∆I n = I n+4 -I n , which has a symmetric PDF (Fig. 3A). The differential intensity ∆I n is digitized to 6 bits (Fig. 3A), and three least significant bits (LSBs) are used for RBG (4). All 8 combinations for 3 LSBs have almost equal probability (Fig. 3B). We remove the residual bias by performing exclusive-OR (XOR) on two bit streams from distant spatial channels (26), which reduces the number of parallel bit streams to 127. Figure 3C reveals the correlation between bits in a single bit stream reaches the limit 1/ √ N given by the bit stream length N . This leads to a single-channel bit generation rate of 2 Tb/s, which is twice the current record with off-line post-processing (13)(14)(15).

We evaluate the quality of the generated random bits with two standard statistical test suites -NIST SP 800-22 and Diehard (26). Fig. 3E provides the test results for 127 parallel bit streams.

95 of them pass all NIST tests, yielding a pass rate of 75%. Compared to the pass rates reported previously for pseudo-RBG and physical RBG, our pass rate is within the acceptable range for reliable random bit generators (16,[START_REF] Lihua | Tenth International Conference on Computational Intelligence and Security[END_REF]. In addition, we perform the Diehard tests on all bit streams (26). The average pass rate over 10 separate tests with different data sets is 93%, comparable to the pass rate of pseudo-RBG.

To investigate the effect of photodetection noise on RBG, we define the signal-to-noise-ratio (SNR) as the bin size s for digitization of ∆I n divided by 2σ of the background fluctuation in a channel. With the SNR significantly higher than 1 for all spatial channels, the random bits are generated predominately by the laser emission, and a numerical estimation of the noise contribution is presented in (26). Figure 3F shows that the percentage of all 188 subtests that every bit stream passes is uncorrelated with the SNR for the pair of spatial channels XOR'ed to create it, indicating the level of detection noise does not affect the random bit quality.

To confirm the absence of correlations among the parallel bit streams, Fig. 3D shows that the mutual information (MI) between any pair of bit streams (26) is as small as the MI of uncorrelated bit streams from different lasers. Moreover, to exclude short-term correlations, we combine odd bits from one stream and even bits from another to generate new sequences. The NIST tests of such combined bit sequences yield a pass rate of 72% -73% (26), demonstrating that all the original parallel bit streams are truly independent.

All these test results certify the randomness of our parallel random bits generated at a cumulative rate of 2 Tb/s ×127 = 254 Tb/s. The very high RBG rate that we obtain indicates an enormous amount of entropy created by our laser. To establish its physical origin, we consider a simple model including only the interference of transverse and longitudinal lasing modes and spontaneous emission noise (26). Using the Cohen-Procaccia algorithm (1), we estimate the entropy rate h CP for a bit stream generated from the simulated intensity fluctuations of a single spatial channel. Figure 4A shows the convergence of h CP for different embedding dimensions d. Both the interference of a large number of lasing modes and the spontaneous emission noise contribute to entropy generation (26). Due to stochastic intensity fluctuations, h CP increases linearly with the number of digits N digit . The fact that h CP reaches the information theoretical limit h 0 (1,26) indicates that the maximal possible bit rate for a single channel is achieved.

This rate exceeds the experimentally obtained value, due to the limited temporal resolution and dynamic range of the photodetector.

To determine how many independent spatial channels are available for parallel RBG, we investigate the effective spatial degrees of freedom (DoFs) of the emission pattern of our laser.

Intuitively, the number of spatial DoFs in the total intensity pattern is expected to be 2M , where M is the number of transverse lasing modes (26). However, the mode amplitudes are not uniformly distributed due to gain competition and saturation, effectively reducing the spatial DoFs.

Applying the Karhunen-Loeve decomposition to the intensity pattern I(x, t), we compute the Shannon entropy H to obtain the complexity as a function of the number of transverse modes M (26, 30). In Fig. 4B, the number of effective DoFs 2 H grows linearly with M , but with a slope smaller than 2. By maximizing M with our cavity design, the maximal number of spatial channels is available for parallel RBG. Keeping only three LSBs after digitizing the emission intensity further reduces the spatial correlation length, and the number of independent channels is thus further increased (26).

In this proof-of-concept experiment, we demonstrate parallel RBG in 127 independent channels with a rate of 2 Tb/s per channel. Real-time streaming of parallel random bits to a computer by conducting the post-processing (including XOR) "on the fly" remains a major technological challenge (16,17).

Besides the application of RBG, the extraordinary spatio-temporal complexity of our laser facilitates rich, diverse dynamical behavior, which can be remarkably tailored via the cavity geometry. By varying the spatial structure of cavity modes and tuning their characteristic length scale, we could effectively manipulate their nonlinear interactions with the gain medium to create deterministic spatio-temporal structures on demand. Such ability of controlling the number of active modes and their nonlinear interactions promotes our laser as a model system to study many-body phenomena and for harvesting spatio-temporal quantum fluctuations. Since our laser possesses a variety of temporal and spatial scales, it may also be useful for studying optical turbulence with high Reynolds numbers. While having a high-dimensional phase space with a complex landscape, our laser is yet compact and may be employed for reservoir computing and for creating physical unclonable functions (PUFs). The simulated power spectrum (black) is broader, but becomes narrower when the temporal resolution of our detector is taken into account (blue), in agreement with the measured one (red). The mutual information between bit streams in two channels with varying separation (Intra-cavity) is equal to that between the streams from two independent lasers (Inter-cavity). (E) The NIST SP800-22 test results include 15 kinds of statistical tests, in total 188 subtests for 127 parallel bit streams. The green (red) color denotes one stream passing (failing) one test. 95 bit streams pass all subtests, yielding a pass rate of 75%, which is considered acceptable for a reliable RBG. (F) The percentage of all 188 subtests that every bit stream passes (green) is uncorrelated with the SNR s/2σ of the corresponding pair of spatial channels (red dashed and dotted lines). Massively parallel ultrafast random bit generation with a chip-scale laser: Supplementary Material

Materials and Methods

• Device fabrication

• Device testing

• Spatio-temporal cross-correlation

• Radio-frequency (RF) spectrum

• Random bit generation

Materials and methods

Device fabrication

We fabricate edge-emitting semiconductor lasers with a commercial GaAs/AlGaAs quantum well (QW) epiwafer (Q-Photonics QEWLD-808). The cavities with curved facets are defined

Device testing

The diode laser is mounted on a copper plate. A tungsten needle (Quater Research H-20242) is placed on the top contact for electric current injection. To reduce sample heating, we use a diode driver (DEI Scientific, PCX-7401) to generate 150-ns-long current pulses at a repetition rate of 7 Hz.

In Fig. 1B, the electric current injected to the stripe laser is 600 mA, where the lasing threshold current is 320 mA. In Figs. 1E, 2, and 3, the pump currents are set at two times the lasing thresholds, where the threshold currents are 400 mA, 600 mA, and 800 mA, respectively. Because the facet is curved, the regions near the boundary of the field of view are out of focus when imaged onto the entrance slit of the streak camera. Thus the measured intensity of emission from these regions is lower than that from the central region of the facet.

As schematically shown in

In the streak images shown in Figs. 1B&E, the time window is 5-ns-long, and the temporal resolution is about 12 ps. Shortening the measurement time window to 0.5 ns in Fig. 2A increases the temporal resolution to 1.2 ps.

Due to the finite temporal measurement range of the streak camera, we concatenate one million streak images from consecutive pulses to obtain a 0.5-ms-long time trace for each spatial position. This process does not artificially increase the entropy generation rate according to our numerical simulation detailed in the supplementary text.

The emission spectrum is measured by an imaging spectrometer (Acton SP300i) equipped with an intensified CCD (ICCD) camera (Andor iStar DH312T-18U-73).

Spatio-temporal cross-correlation

We calculate the correlation function C(∆x, ∆t) of the spatio-temporal emission intensity pat-

tern I(x, t), C(∆x, ∆t) = δI(x 0 + ∆x, t + ∆t)δI(x 0 , t) t δI 2 (x 0 + ∆x, t) t δI 2 (x 0 , t) t (S1)
where δI(x, t) = I(x, t) -I(x, t) t represents the temporal fluctuation of emission intensity at the transverse position x on the end facet, and x 0 denotes the center of the facet (25). In

Radio-frequency (RF) spectrum

We calculate the RF spectrum of the emission intensity by performing the Fourier transform of a time trace at the spatial location x = 0. In Fig. 2C, the measured RF spectra of the lasing emission and background noise are averaged over 10 measurements. To be comparable to the experimental data, the simulated RF spectra are averaged over 10 realizations with different random phases (see the section Numerical modeling).

Random bit generation

The lasing emission intensity in every spatial channel is temporally sampled at intervals of length τ (left inset of Fig. 3A). Faster temporal sampling will give a higher bit rate, but adjacent bits are correlated if the sampling period τ becomes shorter than the correlation time of the fluctuations. The optimal sampling period of τ is found at 1.46 ps (see Fig. S8 and related discussion).

For spatial multiplexing, a smaller channel width ∆x yields a larger number J of spatial channels. However, correlated bit streams will be produced if ∆x is shorter than the spatial correlation length of the lasing emission. A larger width also increases the signal in every channel. We find the optimal width ∆x = 1 µm (see Fig. S13 and related discussion).

In a 2D streak image of the spatio-temporal intensity of the laser emission, each spatial channel of width ∆x = 1 µm contains 4 spatial pixels. The temporal sampling period τ = 1.46 ps corresponds to 3 temporal pixels. We sum the intensities of 4 × 3 = 12 pixels in a streak image to obtain one intensity value I n .

The streak camera background counts fluctuate on the spatial scale of a single pixel. The noise contribution to RBG is determined by its fluctuation, not its mean value which cancels out in the differential intensity. The standard deviation σ of the differential background count (with the laser turned off) is 1.1. When summing over 12 pixels, the standard deviation is increased √ 12 = 3.46 times to σ = 3.9.

Since σ is much smaller than the bin size s = 54 of differential intensity ∆I n in Fig. 3A, the random bits are determined predominately by the laser emission. However, if the value of ∆I n is close to the boundary of one bin, the noise could alter the bit extraction. This is estimated to happen for 3.3% of the bits in the spatial channel shown in Fig. 3A. The percentage can be reduced by increasing the signal strength with better collection of the laser emission.

To further remove any residual bias in the parallel bit streams from 254 spatial channels, we divide the 254 spatial channels into two groups, and perform an XOR between the j-th spatial channel in the first group and the (j + 127)-th spatial channel in the second group, where j is an integer ranging from 1 to 127. Thus the total number of parallel bit streams is 127.

Mutual information

We calculate the mutual information between a pair of bit streams, Y i and Y j , generated in parallel,

h ij = Y i ,Y j p(Y i , Y j ) log 2 p(Y i , Y j ) p(Y i )p(Y j ) . ( S2
)
Here p(Y i ) is the probability density function (PDF) of a random bit stream Y i , p(Y i , Y j ) is the joint PDF of two random bit streams Y i and Y j , and ... denotes an average over all pairs of channels i and j with a constant distance.

NIST tests of random bits

The NIST SP800-22 Random Bit Generator test suite consists of 15 different kinds of statistical tests, some of which include subtests (31). Each test returns a single or multiple p-values.

When the p-value exceeds a significance level of α = 0.01, the bit stream is considered random, as recommended by NIST. For k bit streams, we examine if they pass or fail each statistical test. The pass proportion should be within (1 -α) ± 3 α(1 -α)/k. For each spatial channel, we use k = 1000 bit sequences, each having 2 20 bits, in total over 10 9 bits. As an example, Fig. S15 shows the test results for one bit stream. In panel A, the pass proportions are all above the criterion indicated by the red line.

For a good random bit generator, the p-values from the k bit streams should be uniformly distributed. The composite P-value (p-value of the p-values) is a measure of uniformity of the p-values. The distribution of p-values is considered as uniform when the composite P-value is larger than a significance level of 10 -4 . The example in Fig. S15B shows that the composite P-values of all subtests are above the significance level indicated by the red line.

Diehard tests of random bits

The Diehard battery of randomness tests consists of 18 different kinds of statistical tests (32).

Each test returns a single or multiple p-values. Some of the tests return a large number of p-values, and a composite p-value is calculated by the Kolmogorov-Smirnov (KS) test to determine if the p-values are uniformly distributed in [0,1). A bit sequence is considered random when p-values from all the tests are within the 95% confidence interval of [0.0001, 0.9999] (33).

These tests require about 100 Mbit per channel, and we acquire sufficient data to perform the tests 10 times. The pass rates for all bit streams are listed in Table S1.

Numerical modeling

The cavity resonances are calculated with the eigenmode solver module of COMSOL Multiphysics. To investigate the mode competition for optical gain, we calculate the lasing modes with the single pole approximation steady-state ab-initio lasing theory (SPA-SALT) (34,41,42).

To maximize the number of transverse lasing modes, we optimize the cavity geometry: the ratio of cavity length to width L/W = 1.41, and the radius of the end facets R/L = 0.58 (27). The refractive index n is set to 3.37. The total emitted field is a sum of fields in numerous transverse and longitudinal modes with frequencies within the GaAs QW gain spectrum. The intensity pattern at one facet can be written as

I(x, t) = M -1 m=0 q A m,q e iφm,q(t) ψ m (x)e i2πνm,qt 2
, ( S3
)
where ν m,q is the frequency of a mode with the transverse index m and longitudinal index q, ψ m (x) represents its transverse field profile on the end facet, and A m,q and φ m,q denote its global amplitude and phase. Due to spontaneous emission, φ m,q (t) fluctuates in time, following a Wiener process,

∆φ(t) = φ(t + ∆t) -φ(t) = √ 2πδν∆t Z(t), ( S4
)
where δν is the lasing mode linewidth, ∆t is a discrete time step, and Z(t) is a normaldistributed random number with a standard deviation of 1. For each mode, the initial phase φ(t = 0) is randomly chosen in the interval [0, 2π) with a uniform probability density.

Entropy rate calculation

Using

H d ( , τ ) = - 1 N N j=1 log 2 [f j ( )]. ( S5
)
The Cohen-Procaccia entropy rate estimate is then obtained by minimum of 0, and a maximum of 300. The information theoretical limit h 0 is given by (1,36), From the intensity fluctuation δI(x, t) = I(x, t) -I(x, t) t , the spatial covariance matrix C ab = δI(x a , t)δI(x b , t) t is constructed and its eigenvalues λ α are computed [START_REF] Hess | [END_REF]. λ α is sorted from high to low with the index α, and it reflects the amplitude of the corresponding eigenmode in I(x, t). The value of λ α has a sudden drop at α = 2M , where M is the number of transverse lasing modes. Hence, the spatial degrees of freedom is 2M , where the factor 2 stems from the independent degrees of freedom in the amplitude and phase of the field of one mode.

h CP ( , τ, d) = τ -1 [H d ( , τ ) -H d-1 ( , τ )]. (
h 0 = min(τ -1 , 2f BW ){N digit -D KL [p(I)||u(I)]} (
Applying the Karhunen-Loeve decomposition to the simulated intensity pattern I(x, t), we quantify the spatial complexity by the Shannon entropy of the eigenvalues (30)

H = - α p α log 2 p α , (S8)
where p α = λ α / α λ α is the normalized eigenvalue. The effective degree of freedom 2 H is less than 2M due to the different amplitudes of individual eigenmodes. The lasing threshold current is 570 mA. The inset is a double-logarithmic plot of the LI curve, exhibiting the S-shape that is a hallmark of lasers. (B) Time-integrated spectrum of lasing emission at a pump current of 1200 mA (gray solid line). The FWHM is 1.3 nm. The spectrum S(λ) is fitted by a Gaussian function (red dashed line).

i.e., electric field parallel to the cavity plane. Hence the emission spectrum S(λ) contains 8 longitudinal modal groups within its FWHM.

Lasing dynamics

In a conventional broad-area edge-emitting semiconductor laser, lasing occurs only in the lower order transverse modes because the higher order ones experience stronger diffraction losses and less gain. Furthermore, modulational instabilities induced by nonlinear interactions of the lasing modes with the gain material entail irregular pulsation and filamentation (25,[START_REF] Ohtsubo | Semiconductor lasers: stability, instability and chaos[END_REF]. As shown in Fig. 1B, the emission from a wide-stripe GaAs quantum well (QW) laser is spatially concentrated at multiple locations, forming wire-like streaks called filaments. They originate from carrier depletion in a transverse region of high lasing intensity, which leads to a local increase in the refractive index. The ensuing lensing effect causes self-focusing of the optical field, which evolves into a filament with a transverse size of several microns. The filaments are inherently unstable, and their intensities oscillate irregularly on a sub-nanosecond time scale.

Since the spatio-temporal scales of these self-organized structures are determined by intrinsic properties of the active gain medium (24, 40), they fundamentally limit the parallel RBG rate.

To achieve massively parallel ultrafast RBG, we greatly increase the number of transverse lasing modes and eliminate the long-range spatio-temporal correlations of the emission intensity by tailoring the cavity geometry. In a cavity of increased width and curved facets, high-order transverse modes are effectively confined. The spatial overlap of these modes with the injected carriers is increased by adapting the shape of the top metal contact. Fine tuning of the facet curvature maximizes the number of transverse lasing modes (27). The higher order transverse modes have a transverse wavelength λ t ∼ 1µm, which reduces the lateral width of spatial holes burnt in the carrier density. The resulting refractive index changes induce optical lenses that are too small to focus light and create filaments. In addition, the spatial modulations of the refractive index on such short scales supersede and disrupt the large lenses induced by lower order transverse modes and thus prevent filamentation (28). Consequently, long-range spatiotemporal correlations disappear in C(∆x, ∆t) (Fig. 1F).

Numerical modeling 2.1. Cavity resonances

We calculate the resonances of the passive cavity with the COMSOL eigenfrequency solver module. Since the laser emission is purely TE polarized, we compute only TE-polarized modes which are the solutions of the scalar Helmholtz equation index n is equal to 3.37. We calculate the resonant modes with wavelengths around 800 nm which is similar to the measured laser wavelength (Fig. S2B). However, since the smaller cavity has a larger mode spacing, we extend the wavelength range of simulation to 795-805 nm in order to increase the number of modes, covering 4 longitudinal mode groups.

[∇ 2 + k 2 n 2 (x, y)]H z (x, y) = 0 (
Figure S3A shows the quality (Q) factors and wavelengths of the passive cavity modes. The transverse mode numbers are indicated. In Fig. S3B, the modes are arranged in terms of their transverse and longitudinal mode indices. Although its size is relatively small, the simulated cavity exhibits about 50 transverse modes. The number of transverse modes scales linearly with the cavity size (27).

Lasing modes

Mode competition for optical gain tends to reduce the number of transverse lasing modes.

We calculate the lasing modes using the single pole approximation steady-state ab-initio laser theory (SPA-SALT) (34,41,42), taking gain saturation fully into account. We assume a spatially uniform distribution of the pump and a flat gain spectrum within the wavelength range of 795-805 nm.

The modes that lase in the simulation are marked by red squares in Fig. S3. When the pump is two times the lasing threshold, the number of transverse lasing modes M is 46. Hence, most of the transverse modes in the cavity can lase in spite of gain competition.

Many-mode interference

With spatio-temporal instabilities and filamentation suppressed in our laser, the lasing modes correspond to the high-Q resonances of the passive cavity. To simulate a large cavity of size equal to the fabricated one (L = 600 µm), we compute the mode frequencies ν m,q = c/λ m,q with the analytical expression for Hermite-Gaussian modes [START_REF] Siegman | Lasers[END_REF],

ν m,q = c 2nL q + 1 π m + 1 2 arccos(g) (S10)
where g ≡ 1 -L/R is the cavity stability parameter, which is g = -0.74 for the lasers considered here.

The longitudinal mode spacing is ∆ν q = ν m,q+1 -ν m,q = c/2nL, and the transverse mode spacing is ∆ν m = ν m+1,q -ν m,q = (c/2nL) arccos(g)/π. Their ratio is ∆ν m /∆ν q = arccos(g)/π. It is an irrational number for g = -0.74, making the longitudinal mode spacing incommensurate with the transverse mode spacing. This result is confirmed by our numerical simulation of cavity resonances with COMSOL. We note that the frequency spacings are calculated for the modes of the passive cavity. Experimentally with electric current injection, the spatially inhomogeneous distribution of carriers can cause non-uniform changes of the refractive index, making the lasing mode frequencies deviate from the passive cavity mode frequencies. Nevertheless, the experimentally measured RF spectrum of the lasing emission in Fig. 2C remains broad and smooth, similar to the numerically simulated spectrum of the passive cavity. Due to the presence of hundreds of lasing modes with many different frequency spacings between them, there is an almost continuous distribution of time scales in their interference pattern, so even if some of the time scales were by chance commensurate, the effect would be lost in the sea of other incommensurate ones. Therefore, we ignore such possible deviations in the following calculation.

Summing over all lasing modes, the emission intensity at one spatial position, e.g., the center of the end facet x 0 , is

I(t) = M -1 m=0 q A m,q e i[2πνm,qt+φm,q(t)] 2 (S11)
with A m,q approximated by S(λ m,q ), where S(λ m,q ) is the fit of the measured emission spectrum in Fig. S2B. The total number of transverse modes is M = 200 as in the experiments (27).

To account for the spontaneous emission, we introduce a stochastic fluctuation to the phase φ m,q (t) of each lasing mode [START_REF] Gallion | [END_REF]. The mode linewidth is set to δν = 100 MHz, which is typical for a GaAs/AlGaAs QW edge-emitting laser (45,46). The corresponding coherence time, given by the inverse of the linewidth, is about 10 ns. We set the discrete time step ∆t to 0.1 ps. As shown in Fig. S4A, the phase of each mode undergoes a random walk. The optical spectrum of a single lasing mode, calculated via the temporal Fourier transform of its field, exhibits a Lorentzian-shaped line with FWHM equal to δν as shown in Fig. S4B.

Single-channel dynamics

Figure S5A shows a portion of the simulated time trace I(t) for a single spatial channel. The sampling period is set to 0.5 ps, corresponding to the temporal pixel size of our streak camera. We compute the temporal correlation function for I(t) defined as

C(∆t) = δI(t)δI(t + ∆t) t δI 2 (t) t , ( S12
)
where δI(t) = I(t) -I(t) t is the intensity fluctuation around the mean I(t) t . As shown in Fig. S5B, the half-width-at-half-maximum (HWHM) of the correlation function is 0.5 ps. This is smaller than the measured correlation width (HWHM) of 1.4 ps due to the finite temporal resolution of the streak camera.

To take into account the temporal resolution, we convolve the simulated time trace with a temporal point spread function (PSF). The PSF of the streak camera is approximated by a Lorentzian function with a FWHM of 1.2 ps. The convolution smoothens the time trace (Fig. S5A), and increases the temporal correlation width to 1.4 ps, in agreement with the experimental value (Fig. S5B).

Power spectra

To simulate the radio frequency (RF) spectra shown in Fig. 2C, we calculate the Fourier transform of the simulated intensity I(t). The time trace is 500-ps-long, and the sampling period is 0.5 ps. The power spectra are averaged over 10 time traces with different realizations of phase noise to be consistent with the measurements.

To understand how the broad RF spectrum is formed, we vary the number of transverse modes M while keeping the number of longitudinal mode groups constant. This corresponds to changing the width of a laser cavity while keeping its length fixed. With only the fundamental transverse mode M = 1 (Fig. S6A), the RF spectrum features multiple peaks separated by the longitudinal mode spacing (free spectral range) ∆ν FSR = c/2nL. With increasing M , each peak becomes a group of peaks that originate from the beating of transverse modes. For the cavity with g = -0.74 considered here, the transverse mode spacing is incommensurate with the free spectral range because arccos(g)/π in Eq. (S10) is an irrational number. As a result, the additional beating frequencies that appear when increasing the number of transverse modes eventually fill the entire frequency range. With M = 200 transverse modes in Fig. S6C, the power spectrum becomes continuous and featureless, in agreement with the experimental data in Fig. 2C. The dense packing of the power spectrum increases the entropy generation as will be shown in the next subsection.

Entropy rate

To understand how the number of transverse modes M and the spontaneous emission noise affect entropy generation, we vary M and turn on/off the phase noise when calculating the entropy rate. The number of longitudinal mode groups is fixed at 8 as before. We calculate the original entropy created by multimode interference using a time trace of intensity without accounting for the finite temporal measurement resolution. The length of the simulated time traces is 1.5

µs. With only one transverse mode and without spontaneous emission noise, the intensity trace is periodic in time (Fig. S7A). The periodic modulation results from the temporal beating of the longitudinal modes with equal frequency spacing. Fig. S7D shows the Cohen-Procaccia entropy rate h CP versus the number of digits N digit for embedding dimension d = 3. As N digit increases, h CP first rises then drops, as the intensity trace repeats itself in time and no additional entropy is created eventually. With M = 3 transverse modes, the intensity trace in Fig. S7B exhibits more complex and aperiodic modulations. Since the transverse and longitudinal mode spacings are incommensurate, the temporal beating of 24 modes produces an intensity trace that will never repeat itself. Consequently, h CP keeps increasing with N digit , first linearly then sublinearly in Fig. S7D. Adding the spontaneous emission noise (Fig. S7C) contributes to entropy generation as can be seen in a further increase of h CP in Fig. S7D. Moreover, h CP grows linearly with N digit as a result of the stochastic noise.

We note that the entropy rate is an average quantity obtained from a long time trace. Therefore, it cannot certify the randomness in any short time window, which is required to pass the NIST and Diehard tests. The details of randomness verification of our experimentally generated bit streams are presented in section 4.

Random bit generation

Temporal sampling rate

The black curve in Fig. S8A is the emission intensity at one spatial position of the laser facet that is recorded by the streak camera. The temporal pixel size of the streak camera corresponds to 0.49 ps, and the sampling points are denoted by black dots. The temporal correlation function (Eq. S12) for this trace in Fig. S8B reveals significant correlations among neighboring sampling points. To create independent bits, we must choose a longer sampling period τ , which reduces the sampling rate and hence bit generation rate, however. In order to find the optimal sampling period, we calculate the Cohen-Procaccia entropy rate per sample for the experimental bit streams generated with different τ . The sampled intensity is binned into 2 N digit equally spaced intervals to create N digit bits. As shown in Fig. S8C, the entropy rate per sample first increases with the sampling period τ , then levels off for τ ≥ 1.46 ps. This trend is similar for different N digit . Hence we set the sampling period to be 3 times the streak camera pixel size, τ = 1.46 ps, to extract the maximal entropy per sample with the fastest possible sampling rate.

Intensity PDF

Figure S9A shows the probability density function (PDF) of the sampled intensity I n in one spatial channel. The PDF exhibits an exponentially decaying tail, which is a characteristic of Rayleigh speckle pattern. The PDF is not a perfect exponential function due to the finite spatial produces a bit stream with long-range correlations, which also degrade the random bit quality.

We choose m = 4 as an optimal sample distance. Figure S9B shows the PDF of the differential intensity from experimental data. It is fit well by a Gaussian function, which leads to the equal probability of 8 possibilities for 3 LSBs in Fig. 3B.

Suppressing temporal correlations

As seen in Fig. S8B, even after choosing a sampling period of τ = 1.46 ps, the correlation of adjacent samples is not completely eliminated. RBG demands negligible correlation between successive bits. Here digitization and post-processing play a crucial role in removing the remaining correlation. During the analog-to-digital conversion (ADC), the measured intensity is transformed to N digit = 6 bits, and only the three least significant bits (LSBs) are kept. In Fig. S10, we compare the temporal correlation function [Eq. (S12)] of the bit stream to that of the original intensity trace, averaged over all channels. For long time lags, keeping only the LSBs reduces the correlation below 10 -3 , which is the lower limit given by the finite length of the bit stream (2 20 samples). For short time lags (Fig. S10B), the correlation for the LSBs black dotted line in Fig. S11A) and are thus out of focus. Consequently, the emission intensity drops rapidly near the two edges of the FoV (red curve in Fig. S11B).

In Fig. S11B, we compare the measured emission profile to the ideal one obtained from numerical simulation. Using SPA-SALT, we calculate the lasing modes in two small cavities of lengths L = 20, 40 µm with the same geometry as the laser cavities in the experiment.

After rescaling, the calculated intensity distribution on the curved facet is identical for the two cavities. Hence, the emission profile is universal and scales linearly with the cavity size. The FWHM of the measured profile is about half of the simulated one. Hence, only half of the available spatial channels are used for random bit generation in the current experiment.

Spatial correlation length

The spacing between independent spatial channels is determined by the spatial correlation length of the emission intensity. To find its relation to the transverse wavelength of the lasing modes, we simulate the spatio-temporal intensity pattern of the lasing emission from a small cavity with L = 40 µm, W = 28.2 µm and R = 23 µm (see materials and methods). The number of transverse lasing modes M at a pump level of two times the lasing threshold is 46 (see Fig. S3). Figure S12A is a portion of the spatio-temporal intensity pattern at one cavity facet with full spatial resolution. In the experiment, the numerical aperture (NA) of the imaging op-tics is 0.4, and the spatial resolution of the streak camera is about 0.5 µm. The convolution with the spatial point spread functions (PSF) of imaging optics and streak camera enlarges the spatial speckle grains (Fig. S12B). In Fig. S12C, we show the correlation of the intensity fluctuations δI(x, t) = I(x, t) -I(x, t) t at spatial locations separated by ∆x,

C(∆x) = δI(x, t)δI(x + ∆x, t) t δI 2 (x, t) t δI 2 (x + ∆x, t) t x . ( S13
)
The HWHM of C(∆x) is 0.25 µm for the fully resolved pattern in Fig. S12A. The average speckle grain size (FWHM) is 0.5 µm, which is close to half of the transverse wavelength of the highest-order transverse lasing mode. The limited spatial resolution of imaging optics and streak camera enlarges the speckle grain size to 1.5 µm, thus reducing the number of independent channels.

Suppressing spatial correlations

Digitization and post-processing reduce the spatial correlation length, similarly to the suppression of temporal correlation in Fig. S10. In Fig. 3D, the mutual information of bit streams from two spatial channels becomes negligible when their separation ∆x exceeds 1 µm. However, the original intensity pattern in Fig. 2 shows a spatial correlation extending over a distance of 1.5

µm. In Fig. S13 we compare the mutual information (MI) between two bit streams produced experimentally for three different cases: (i) thresholding N digit = 1, the simplest bit-extraction scheme; (ii) keeping 3 LSBs from analog-to-bit conversion with N digit = 6; and (iii) conducting XOR of (ii) with a bit stream from a spatially distant channel. In comparison to (i), the MI between neighboring channels (with 1 µm spacing) is reduced by five orders of magnitude in (ii). Moreover, the MI is further reduced at short-range in (iii). It stays at the residual level of 10 -6 and becomes independent of the channel separation. all channels are statistically independent and their residual MI is a result of the finite bit stream length.

Effect of noise on binning

For every spatial channel, the entire binning range of its differential intensity ∆I n , which is set to 8.2 times the standard deviation of ∆I n , is divided into 2 6 = 64 equally sized bins. If the value of ∆I n is close to the boundary of one bin, the noise can alter the bit extraction.

To examine how strong this effect is, we simulate the streak camera noise and add it to an experimental intensity trace. We compare the bit stream generated by the noise-altered trace to the original one. With 3 LSBs extracted from 6 digits, 3.3% of all the bits in the spatial channel shown in Fig. 3A are altered by the addition of noise. For all 254 spatial channels, the average bit error rate is 4.1 ± 1.2%. This percentage can be reduced by increasing the signal strength with better collection of the laser emission.

Concatenating time traces

Since the temporal measurement range of our streak camera is limited, it is impossible to measure a long continuous time trace. Instead we make separate measurements and concatenate the time traces. Since this process could potentially increase the randomness, we check numerically whether the entropy generation rate is changed by it.

Using Eq. (S11) to simulate the many-mode interference, we obtain a 1.5-µs-long time trace of the emission intensity in a single spatial channel. We repeat this process with different random phases for each mode, and obtain 3000 time traces. Short segments of three such traces are shown in Fig. S14A. Then we extract 500-ps-long segments in consecutive time windows from each trace, and concatenate these segments together for a trace of length 1.5 µs (Fig. S14B).

We calculate the Cohen-Procaccia entropy rate of the original trace and the concatenated trace. Figure 3E shows that 75% of the parallel random bit streams can pass the entire NIST test.

Considering the statistical nature of the NIST tests, the pass rate was previously evaluated for pseudo-random number generators and physical random number generators. In Ref. (47), the NIST test was applied to 100 different sample data sets (each set consisting of 1000 segments of 1 Mbit length), and a pass rate of 41%-56% was obtained for various representative pseudo-

Figure 1 :

 1 Figure 1: Reducing spatio-temporal correlations of the lasing emission. (A) A wide-stripe edge-emitting semiconductor laser with planar facets supports only low-order transverse modes with the typical spatial profile shown (not to scale). (B) Emission intensity I(x, t) at one facet of a 100-µm-wide, 1000-µm-long GaAs quantum-well laser, measured by a streak camera, features filamentation and irregular pulsations. (C) The spatio-temporal correlation function C(∆x, ∆t) of the emission intensity in (B) reveals long-range spatio-temporal correlations. (D) Our specially-designed laser cavity with curved facets confines high-order transverse modes. The spatial intensity distribution of an exemplary high-order transverse mode is plotted. (E) The measured spatio-temporal trace of the lasing emission from our cavity of length 400 µm, width 282 µm and facet radius 230 µm is free of micron-sized filaments and GHz oscillations as seen in (B). (F) The spatio-temporal correlation function C(∆x, ∆t) of the emission intensity in (E) shows no long-range spatio-temporal correlations.

Figure 2 :

 2 Figure 2: Ultrafast beating of lasing modes. (A) The lasing emission at one facet of our 600µm-long cavity exhibits a spatio-temporal interference pattern. The white solid curve is the temporal intensity fluctuation at position x = 0 (white dashed line). (B) The correlation function C(∆x, ∆t) of the measured lasing emission intensity in (A) gives the spatial and temporal correlation widths of 1.5 µm and 2.8 ps. (C) The power spectrum of the emission intensity in (A) at x = 0 (red) is much higher than the background (gray) with the laser turned off. The simulated power spectrum (black) is broader, but becomes narrower when the temporal resolution of our detector is taken into account (blue), in agreement with the measured one (red).

Figure 3 :

 3 Figure 3: Parallel random bit generation and evaluation. (A) Left inset is a segment of intensity time trace of a single spatial channel (red line), sampled at intervals τ = 1.46 ps (red dots). The blue curve is the background count. Main panel shows the PDF of the differential intensity, ∆I n = I n+4 -I n , which is digitized to 6 bits by binning the range [-1740, 1740] counts (vertical red dashed lines) into 64 equally spaced intervals. 3 LSBs are taken from each sample. The gray scale of the bars represents their 8 combinations. Right inset is the PDF of the differential background count with a standard deviation σ = 3.9 counts much smaller than the bin size s = 54 for ∆I n . (B) The probability for all 8 combinations of 3 LSBs is almost equal. (C) A bit stream with length N = 2 20 has a bit correlation (red squares) around the lower limit 1/√ N (black line). Inset is zoom on short delay times. (D) The mutual information between bit streams in two channels with varying separation (Intra-cavity) is equal to that between the streams from two independent lasers (Inter-cavity). (E) The NIST SP800-22 test results include 15 kinds of statistical tests, in total 188 subtests for 127 parallel bit streams. The green (red) color denotes one stream passing (failing) one test. 95 bit streams pass all subtests, yielding a pass rate of 75%, which is considered acceptable for a reliable RBG. (F) The percentage of all 188 subtests that every bit stream passes (green) is uncorrelated with the SNR s/2σ of the corresponding pair of spatial channels (red dashed and dotted lines).

Figure 4 :

 4 Figure 4: Information capacity of spatio-temporal intensity pattern of simulated lasing emission. (A) The entropy rate h CP in a single spatial channel converges for different embedding dimensions d, and reaches the information theoretical limit h 0 (black dashed line), indicating the maximum possible RBG rate is reached. (B) The effective number of spatial degrees of freedom 2 H in the emission pattern grows linearly with the number of transverse lasing modes M . The black solid line is a linear fit. The blue dashed line is the naively expected number of spatial degrees of freedom 2M , which exceeds 2 H .

2 .

 2 by photolithography and etched by an inductively coupled plasma reactive ion etcher. The etch depth is 4 µm. The top metal contacts are fabricated by photolithography, Ti/Au deposition, and lift-off (see Ref. (27) for details). A scanning electron microscope (SEM) image of a fabricated device is presented in Fig. S1B. The cavity length L varies from 400 µm to 800 µm. The transverse width W scales with L as W = L/ √ The radius of the curved facets R is also proportional to L, so that the cavity stability factor g = 1 -L/R is fixed to -0.74. In Figs. 1E, 2, and 3, the cavity lengths are L = 400 µm, 600 µm, and 800 µm, respectively. The wide-stripe lasers with planar facets (g = 1) in Fig. 1B are made by wafer cleaving.

Figure S1 :

 S1 Figure S1: Schematic of a chip-scale parallel random bit generator. (A) Sketch of a manymode electrically-pumped semiconductor laser with curved end facets. Its emission is measured by two arrays of photodetectors for parallel RBG. (B) Top-view SEM image of a fabricated laser with cavity length L = 400 µm, width W = 282 µm, and end facet radius R = 230 µm. (C) Sketch of our experimental setup. The lasing emission at one end facet is imaged onto a streak camera that records the temporal fluctuation of the emission intensities at many spatial locations simultaneously.

Figs

  Figs. 1C&F, C(∆x, ∆t) is averaged over 10 measurements.

  the Cohen-Procaccia algorithm (1, 35), we compute the entropy rate h CP as a function of the bin size for intensity digitization and the temporal sampling period τ . For a time trace of emission intensity I(t) at a single position of the laser facet, we construct d-dimensional data sets by introducing time delays: I 1 = I(t), I 2 = I(t + τ ), ... , I d = I(t + (d -1)τ ). Then we randomly select N reference points in the d-dimensional space. For each reference point j, we compute f j ( ), the fraction of other points within a d-dimensional box of width . The d-dimensional pattern entropy estimate is given by

  S6) Here = (I max -I min )/2 N digit , with I max and I min being the maximum and minimum intensities. In Fig. 4A, h CP is plotted as a function of N digit for different d. The time trace of emission intensity in a single channel is numerically calculated with the cavity parameters identical to the experimental ones (L = 600 µm, W = 424 µm, R = 345 µm). There are 8 longitudinal mode groups within the emission spectrum. The number of transverse modes is M = 200. The optical linewidth of each mode is δν = 100 MHz. The sampling period is τ = 1.5 ps. The temporal range is 1.5 µs, yielding a time trace of 10 6 samples. To have the same dynamic range as the experimental data, the intensity values are rounded to an integer number with a mean of 60, a

  S7) where f BW denotes the signal bandwidth, and D KL [p(I)||u(I)] = I p(I) log 2 [p(I)/u(I)] is the Kullback-Leibler divergence (37) between the intensity PDF p(I) and the uniform PDF u(I) within the same range of digitization.Karhunen-Loeve decompositionWe perform a Karhunen-Loeve decomposition of the simulated lasing intensity pattern I(x, t) in a cavity of length L = 40 µm, width W = 28.2 µm, radius of end facets R = 23 µm and refractive index n = 3.37. The mode frequencies ν m,q and spatial field profiles ψ m (x) are obtained from the COMSOL calculation of cavity resonances. The amplitudes A m,q of individual lasing modes are calculated with the SPA-SALT. Their phases φ m,q are random numbers in the range of [0,2π).

Figure S2 :

 S2 FigureS2: Many-mode laser characteristics. (A) LI curve: total emission power (from both facets of the laser) as a function of the electric current injected into the 600-µm-long cavity. The lasing threshold current is 570 mA. The inset is a double-logarithmic plot of the LI curve, exhibiting the S-shape that is a hallmark of lasers. (B) Time-integrated spectrum of lasing emission at a pump current of 1200 mA (gray solid line). The FWHM is 1.3 nm. The spectrum S(λ) is fitted by a Gaussian function (red dashed line).

Figure S3 :

 S3 Figure S3: Cavity resonances and lasing modes. (A) Calculated Q factors and wavelengths of the resonant modes (black open squares) in a cavity with length L = 40 µm, width W = 28.2 µm, curvature radius of the end facets R = 23 µm, and effective refractive index n = 3.37. The transverse mode number m is written next to each mode. The longitudinal mode spacing (free spectral range) ∆λ q and the transverse mode spacing ∆λ m are indicated by arrows. The modes that lase for pumping at two times the lasing threshold are marked by red solid squares. (B) The cavity resonances (black open squares) and lasing modes (red solid squares) shown in (A) are arranged in terms of their longitudinal and transverse mode numbers. Almost all transverse modes lase even in the presence of gain competition.

Figure S4 :

 S4 Figure S4: Modeling of phase noise. (A) Simulated random phase drift in time for two lasing modes. The inset shows the phase change ∆φ t after a single time step ∆t. (B) The simulated optical spectrum of a lasing mode is fitted by a Lorentzian function with a FWHM of 96 MHz, which is close to the mode linewidth δν = 100 MHz.

Figure S5 :

 S5 Figure S5: Temporal fluctuations and correlation function. (A) The time trace of numerically simulated emission intensity in a single spatial channel (blue dotted line) is smoothed by convolution with the PSF of the streak camera (red solid line). The symbols represent sampling points separated by 0.5 ps. (B) Temporal correlation functions for the two traces in (A). The HWHM values are indicated. The black dashed line is the experimentally measured temporal correlation function.

Figure S6 :

 S6 Figure S6: Simulated power spectrum. Fourier transform of the simulated emission intensity with the number of transverse modes (A) M = 1, (B) M = 6, and (C) M = 200. The number of longitudinal mode groups is fixed at 8. The power spectrum becomes more and more densely packed with increasing M .

Figure S7 :

 S7 Figure S7: Cohen-Procaccia entropy rate of simulated time traces. (A) A portion of the simulated intensity time trace I(t) with only the fundamental transverse modes (M = 1), and without spontaneous emission noise. (B) I(t) with M = 3 and without spontaneous emission noise. (C) I(t) with M = 3 and with spontaneous emission noise (δν = 100 MHz). (D) Entropy rate h CP for the time traces in (A)-(C) with embedding dimension d = 3. The sampling period τ is 1.5 ps.

Figure S8 :Figure S9 :

 S8S9 FigureS8: Selection of temporal sampling period. (A) Measured time trace of the emission intensity with a sampling period of τ = 0.49 ps (black), corresponding to a single pixel of the streak camera. To reduce the temporal correlation of sampled intensity, τ is increased to 1.46 ps (red). The intensity is integrated within the sampling period τ to increase the signal-to-noise ratio. (B) The temporal correlation functions for the two traces in (A) as a function of the time lag in units of the sampling period τ show a faster decay for the longer sampling period τ . (C) Cohen-Procaccia entropy rate of a measured time trace for different sampling periods τ and numbers of digits N digit . The embedding dimension d is set to 3.

Figure S10 :Figure S11 :

 S10S11 Figure S10: Temporal correlation of bits in a single channel. (A) The magnitude of the temporal correlation function |C(∆t)| of the measured intensity trace sampled with τ = 1.46 ps (black circles), after keeping only the 3 LSBs (blue triangles), and after performing the XOR operation in addition (red squares). The correlation functions are averaged over spatial channels. The number of data points in time is 2 20 . (B) The magnified view for short time lags shows the quick decay of correlations for the three LSBs. The small peak at 5.84 ps (= 4τ ) is attributed to the subtraction of the sampled intensity ∆I n = I n+4 -I n .

Figure S12 :

 S12 Figure S12: Spatial correlation of simulated intensity fluctuations. (A) The fully resolved spatio-temporal pattern of emission intensity on the end facet of the laser cavity features fine speckles. (B) The intensity pattern after convolution of (A) with the spatial PSFs of the imaging optics (NA = 0.4) and the streak camera (Lorentzian with FWHM = 0.47 µm) shows an increase of the spatial speckle size. (C) The spatial correlation functions C(∆x) for the intensity fluctuations in (A) (blue dotted line) and (B) (red solid line) have a HWHM of 0.25 µm and 0.75 µm, respectively. The black dashed line denotes the measured spatial correlation.

Figure S13 :

 S13 FigureS13: Mutual information between two spatial channels. (A) Mutual information (MI) of experimental bit streams produced by (i) thresholding (red squares), (ii) keeping 3 LSBs from N digit = 6 bit conversion (blue triangles), (iii) calculating XOR of (ii) with a bit stream from a spatially distant channel (black circles). The MI is averaged over all channels. The MI decreases with the channel separation. (B) The residual mutual information of (iii) in (A) at channel separation ∆x ≥ 1 µm is inversely proportional to the bit stream length N , indicating it results from the finite stream length. The red circle indicates the length of the bit stream used in (A).

Figure S14 :

 S14 Figure S14: Concatenating time traces. (A) Three continuous time traces of the emission intensity in a single spatial channel are generated by simulating the many-mode interference in the laser. The shaded areas represent the temporal range measured by the streak camera. (B) A time trace is formed by concatenating the shaded parts of the three traces in (A). (C) The Cohen-Procaccia entropy rate estimate of a single 1.5-µs-long trace as in (A) is equal to that of the concatenated trace as in (B). The sampling period τ is 1.5 ps, and the embedding dimension d is 2.

Fig. S14C shows

  Fig.S14Cshows complete agreement of the two curves, indicating that the process of concatenating windows from different time traces does not affect entropy generation. We also concatenate the same time windows of different traces, and the Cohen-Procaccia entropy rate is the same as well. As the original time trace already reaches the maximal entropy rate (Fig.4A), concatenating separate traces cannot increase the entropy rate any further.

4 .

 4 Random bit evaluation 4.1. NIST tests We evaluate the quality of the generated random bits with the NIST SP 800-22 test suite. Figure S15 shows the results of NIST tests conducted on a single bit stream containing 1000 sequences with 1 Mbit length per sequence (in total 1 Gbit). The pass proportions and the composite P-values are all above the criteria recommended by NIST, indicating that the quality of the experimentally generated random bit stream is acceptable.

Acknowledgments H.C. and K.K. thank Rajarshi Roy and Marc Sciamanna for stimulating discussions. We acknowledge the computational resources provided by the Yale High Performance Computing Cluster (Yale HPC). Funding: This work is supported by the US National Science Foundation under the Grant No. ECCS-1953959 and by the Office of Naval Research under the Grant No. N00014-21-1-2026. For the work at Nanyang Technological University, funding support is acknowledged from the National Research Foundation Competitive Research Program (NRF-CRP-18-2017-02 and NRF-CRP-19-2017-01) and A*Star AME Programmatic Grant under Grant A18A7b0058). Support by the Science Foundation Ireland (SFI) under grant 18/RP/6236 is gratefully acknowledged. Author contributions: K.K. conducted experiment and data analysis, and prepared the manuscript; S.B. built the experimental setup and guided data analysis; Y.Z. and Q.J.W. did sample fabrication; S.G. and O.H. guided the analysis of semiconductor laser dynamics; and H.C. conceived the idea, initiated the project and supervised the research.

Competing interests. The authors declare no competing financial interests. Data and materials availability. All data needed to evaluate the conclusions in the paper are present in the paper or the Supplementary Materials. We test lasing in the broad-area semiconductor lasers with curved facets. Similar LI curves and emission spectra are obtained for multiple lasers of identical geometry but varying size. Fig. S2A shows the measured LI curve for a cavity of length L = 600 µm, width W = 424 µm, and radius of the curved facets R = 345 µm. The lasing threshold current is 570 mA. At a pump current of 1200 mA, the power of lasing emission collected by the objective lens (NA = 0.4) from one end facet of the cavity is 93 mW. The collection efficiency is estimated to be about 20%, with the divergence angle of the lasing emission, the transmission and numerical aperture of the objective lens, and the collection from a single facet of the cavity all taken into account.

Supplementary Materials

Materials and Methods

Supplementary Text

Thus the total emission power is about 470 mW, which corresponds to a quantum efficiency of 0.74 W/A. We record the far-field speckle pattern created by the output beam passing through a diffuser and calculate its intensity contrast (27). The number of transverse lasing modes M estimated from the measured speckle contrast of 0.07 is about 200. Hence the emission power per transverse mode is on the order of 1 mW.

Figure S2B shows the measured emission spectrum at two times of the lasing threshold integrated over a 0.6-µs-long pulse. Individual lasing modes cannot be resolved, as the wavelength separation between adjacent modes is smaller than the resolution of the spectrometer. The fullwidth-at-half-maximum (FWHM) of the emission spectrum is 1.3 nm. The spectrum S(λ) is fitted by a Gaussian function centered at λ 0 = 800 nm with a standard deviation of 0.55 nm. The free spectral range (longitudinal mode spacing) of the 600-µm-long cavity is λ FSR = λ 2 0 /2nL = 0.16 nm, where the effective refractive index n = 3.37 is calculated from optical confinement in the direction perpendicular to the cavity plane for transverse electric (TE) polarization, random bit generators. In Ref. [START_REF] Lihua | Tenth International Conference on Computational Intelligence and Security[END_REF], a pass rate of 59%-71% was reported for some wellknown pseudo-random number generation algorithms including those recommended by NIST.

In Ref. (16), a pass rate between 65% and 75% was obtained for a physical RBG based on a chaotic laser. By accounting for the correlations of the sub-tests included in the NIST test suite, the upper bound of the pass rate was estimated to be 80.99% [START_REF] Lihua | Tenth International Conference on Computational Intelligence and Security[END_REF]. Compared to the pass rates in these prior studies, our pass rate of 75% for all channels (Fig. 3E) is considered acceptable for reliable random bit generators.

To verify the independence of the random bit streams from different spatial locations, we combine the bit streams that are generated in parallel, using the procedure illustrated in Fig. S16A.

We first test pairs from spatially adjacent locations which have potentially stronger correlations than non-neighboring pairs. Among all bit streams (after XOR operation on two spatial channels), the even (or odd) bits from the (2m -1)-th bit stream and odd (or even) bits from the 31 The entire test suite passes with a 95% confidence interval for p-values between 0.0001 and 0.9999 (33). Figure S17A shows the Diehard test results for 100 Mbit from a single channel.

The p-values from all the statistical tests are within the interval [0.0001, 0.9999], thus the entire test suite is passed. the Diehard test, we evaluate the pass rate over all bit streams. Among the 127 bit streams, 118 completely pass the Diehard test, yielding a pass rate of 93%. We repeat the tests for 10 different sets of data, and the pass rates are listed in Table S1. They range from 91% to 95%, and the average pass rate is 93 ± 1%. As a reference, we repeat the tests with random bits generated by one of the widely used pseudo-RBG algorithms -Mersenne-Twister. With the same amount of random bits, the average pass rate over 10 tests is 92 ± 2%. It is very close to the pass rate of all bit streams generated by our laser, thus certifying the randomness of our parallel RBG.


Jakob J Verbeek
email: jverbeek@wins.uva.nl

Using a Sample-dependent Coding Scheme for Two-Part MDL

In this paper we report on our investigation on how using sample-dependent c o d i n g s c hemes can lead to poor results in applying Rissanen's Minimum Description Length (MDL) principle Ris89]. The MDL principle is one of the many k n o wn model selection methods in the eld of 'machine learning', 'statistics' or 'inductive inference'. We analyze the experimental results presented in KMNR97] and provide a method to avoid the over tting. We do so by using a di erent c o d i n g s c heme than in KMNR97].

Introduction

Rissanen's Minimum Description Length (MDL) principle Ris89] is one of the many known model selection methods in the eld of 'machine learning', 'statistics' or 'inductive inference'. The general idea of MDL can be explained as follows : We w ant to model some regularities or dependencies in the data. These regularities in the data can be used to compress the data. Models describe regularities in the data. The more regularity of the data some model captures, the more it can compress the data. Therefore, the model that enables us to compress the data the most is the best model. Two-Part MDL (TP-MDL) is one way to use this idea. What we do is describe the data in two parts, hence the name. First, we describe the model. Second, we describe the data using the model. The model describes the regularities in the data, the description of the data using the model codes the deviation of the data from these regularities.

In the following we w i l l assume that models can be ordered by their complexity. For now think of model complexity a s t h e n umber of parameters in a model.

Suppose we w ant to model some data generating system and we are given more and more examples of the behavior of the system. TP-MDL would typically start selecting simple models when the sample size is still small. As more and more examples become available TP-MDL usually selects more and more complex models, until it converges to some constant model complexity. 1 However, in the experiments described in KMNR97] TP-MDL selected very complex models until some numberm of examples had been reached. The model complexity o f the selected models converged very rapidly if the number of examples was greater than m . The experimental results obtained for TP-MDL in KMNR97] surprised us. We had two main questions: Would TP-MDL be applied correctly in KMNR97]? Can we devise a coding scheme such that the performance of TP-MDL is improved?

Experiments

In this section we discuss the experiments we conducted to answer the second question. Furthermore, we give some arguments why w e think one should not use sample-dependent c o d i n g s c hemes for TP-MDL.

Problem Domain

The interval function selection problem, as de ned in KMNR97], can be de ned as follows. The data source is a function f d : 0 1] ! f 0 1g that divides the interval 0 1] into d + 1 equal intervals, each labeled zero or one. We call such functions interval functions. The examples are drawn according to a uniform distribution on X = 0 1]. The data sequences consist of pairs (x 2 X y2 Y = f0 1g). The value f d (x) is zero in the rst and all other odd intervals, it is one in all other intervals. The task of the inference algorithms is to output some interval function. Ideally the algorithm should output the model that models the data generating function the best. The training data is corrupted with some noise to make the setting more realistic.

Coding Schemes

Now, we g i v e the coding scheme that was used in KMNR97] to code the interval functions, we call this coding scheme K MD L . KMDL In KMNR97] the set of possible parameter values is restricted to the x-data. This is motivated by the observation that the x-data is guaranteed to yield optimal parameter values. So in this coding scheme the data is not only used to select a h ypothesis, the data is also used to describe a h ypothesis. To c o d e a model using this coding scheme, we rst code the number of parameters t. Next, we c o d e t h e i n d e x o f t h e function in an enumeration of all interval functions with t parameters, such that the parameters are restricted to the x-data. That is: we code the index of the function in an enumeration of all subsets with t elements of the 1 See BRY97]. Why should we use a sample-independent c o d i n g s c heme to code the models? We g i v e t wo arguments below, followed by the sample-independent coding scheme we used.

1. When we use a sample-dependent coding scheme to code the models, we do not describe a model that maps some arbitrary x-data to the Y -domain. What we obtain if we decode the codeword generated by K MD L is a function that gives us a function f : X ! Y if we p r o vide a set of x-data points. Thus, K MD Lcodes a function f 0 : X m ! (X ! Y ), where X m is the set of sets with m values from the X-domain. Usually, we use a function f : X ! Y , as these functions give us a dependency between X and Y , the dependency we w ant to model. 2. If we w ant to use the KMDL coding scheme to code a model (dependency between X and Y ), we should provide the function f 0 and the x-data of the sample. So, KMDL provides us with an e cient way to code a model, given some speci c x-data. However, KMDL is quite an ine cient w ay to code a model for some unspeci ed x-data set. So we can use the model to compress the observed data, but in general it will not enable us to compress some other data. We i n tuitively would like t o h a ve a model that allows us to compress any data (given that it is generated by the same data generating function). Now, we g i v e a sample-independent coding scheme to code models.

MDLpss In this coding scheme we do not restrict parameter values to the x-data. First, we c o d e t h e n umber of bits p, that is used to code one parameter value. Second, code the number of parameters t of the hypothesis. Third, code the index number of the function in the standard lexicographical enumeration of all functions with t parameters encoded using p bits precision.

To e v aluate the performance of TP-MDL using these coding schemes we conducted several experiments.

Results and Conclusions

Using M D L pss improves performance in two w ays (see Figure 1): First, the over tting that resulted from using K MD Ldoes not occur, instead under tting occurs. Second, the number of examples needed to output the target number of parameters is considerably smaller when using M D L pss.

However, using M D L pss, the quality o f parameter estimation decreases dramatically. With the quality o f parameter estimation we m e a n h o w good the algorithm estimates the parameter values. Due to the minimization of code length the number of bits used by M D L pss to code the parameters will be relatively small for small samples. The di erence is even so dramatic that K MD Loutperforms M D L pss in terms of generalization error. However, we can compensate for this by using M D L pss only to select the number of parameters for the model and using Maximum Likelihood to select the actual parameter values. Further research should be carried out to nd out whether this is in general a good method. A more elaborate report on this research i s g i v en in Ver98].

Figure 1 :

 1 Figure1: Plot of sample size, m, against hypothesis complexity, using a noise-rate of = 0 :2.