1. Introduction

The study of Supercontinuum (SC) generation in optical fiber has been a field of intense research over the last decade, leading to a number of important advances in both fundamental and applied science [1–3]. Although the physical mechanisms leading to SC spectral broadening are now well understood [4], recent studies have focussed on the instability properties of the SC. These results have yielded greatly improved insight into the noise-sensitivity of the nonlinear dynamics of SC generation, particularly with regard to establishing intriguing links with the study of extreme events and rogue waves in other systems [5–8].

Another field of optics research that has seen comparable recent interest has concerned the inherent randomness of nonlinear dynamical systems. A particular motivation here has been the generation of random numbers using physical (rather than algorithmic) approaches, as there are many applications of physical random numbers in information theory, cryptography, Monte Carlo simulation and so on. The advantage of optical techniques is that they can exploit a physical random process to generate random numbers at high repetition rate and at optical wavelengths directly compatible with future demands of all-optical integration. Examples of such physical optical random number generators include chaotic lasers and optoelectronic systems, photon counting and homodyne detection of vacuum fluctuations, and spontaneous emission and superluminescent diodes [9–15]. Although the results obtained have been promising, the systems studied have generally operated over a limited wavelength range, whereas future needs of optical random numbers may require physical generators at essentially arbitrary wavelengths.

To this end, we present in this paper a numerical study exploring the potential of the broadband optical SC as a source of physical randomness. Using numerical simulations of SC generation in the incoherent regime, we study the randomness inherent in SC generation from two different perspectives. In particular, although we indeed demonstrate the potential of SC fluctuations for physical random number generation, we begin with a more general discussion of the randomness properties of SC generation, showing that the intensity and phase fluctuations of the SC can be interpreted in terms of the random processes of random walks and Lévy flight-like evolution. We stress that although the intensity and phase fluctuations of the SC spectra that we study are computed numerically, convenient experimental procedures are readily available for practical implementation: intensity measurements are straightforward, and optical phase measurement is enabled through heterodyne detection (with a local oscillator). Since SC generation has been reported over a wide parameter range and using pump sources from the MHz-GHz range, our results suggest that the SC should provide a versatile platform for the study and application of random processes at optical wavelengths. Our aim is to anticipate possible future studies of random walk processes and random number generation using optoelectronic implementations, and to show that the supercontinuum can provide a convenient physical source of random fluctuations that can be applied for this purpose.

2. Supercontinuum intensity and phase fluctuations

We first review the general properties of SC intensity and phase instability under conditions where the dynamics are well known to lead to severe shot-to-shot fluctuations [4]. Specifically, we consider results from noise-seeded generalized nonlinear Schrödinger equation (GNLSE) simulations where SC spectral broadening is seeded by picosecond pulses, and where the dynamics are dominated by noise-driven modulation instability [6, 16].

The general characteristics of these fluctuations can be seen clearly by inspecting the numerical results shown in Fig. 1. The parameters used in the simulations correspond to a realistic photonic crystal fiber with zero dispersion wavelength at 1060 nm. We consider 3 ps pulses of 200 W peak power injected at 1064 nm slightly in the anomalous dispersion regime, and we consider a 20 m fiber length where the nonlinearity coefficient is γ_NL = 11 W⁻¹ km⁻¹ and the dispersion coefficients at the pump wavelength are: β₂ = −4.10 × 10⁻¹ ps² km⁻¹, β₃ = 6.87 × 10⁻² ps³ km⁻¹, β₄ = −9.29 × 10⁻⁵ ps⁴ km⁻¹, β₅ = 2.45 × 10⁻⁷ ps⁵ km⁻¹, β₆ = −9.79 × 10⁻¹⁰ ps⁶ km⁻¹, β₇ = 3.95 × 10⁻¹² ps⁷ km⁻¹, β₈ = −1.12 × 10⁻¹⁴ ps⁸ km⁻¹, β₉ = 1.90 × 10⁻¹⁷ ps⁹ km⁻¹, β₁₀ = −1.51 ×10⁻²⁰ ps¹⁰ km⁻¹. The optical shock timescale used was τ_shock = 0.56 fs [4].

 

Fig. 1 (a) Individual output spectra from each of 1000 individual realizations (gray) and calculated mean (black). The mutual coherence is plotted in the upper subfigure. Histograms from all 200000 events calculated from the energy filtered from a 20 nm bandpass filter with central wavelength (b) λ =880 nm; (c) λ =1100 nm; (d) λ =1270 nm.

Download Full Size | PDF

In the results that follow we have used a noise model of a one photon per mode background with random phase, but near-identical results are obtained using intensity noise in the time domain; in fact, the qualitative nature of the fluctuations observed are largely independent of the nature of the noise model used [17]. Of particular interest in our work, however, is the fact that there is noise on the input field over only a 45 nm bandwidth about the pump wavelength. Outside this bandwidth, the noise level is at the numerical precision of the computation. The use of a finite bandwidth for the initial noise is so that we can can be sure that the randomness we study at any particular wavelength is not simply amplification of initial noise at that wavelength, but rather due to the intrinsic chaotic nature of the nonlinear dynamics [18] that transfers the noise around the pump to a much broader wavelength range across the SC. Similar noise transfer processes are common in other systems in nonlinear optics such as amplifiers and wavelength conversion [19, 20], but to our knowledge have not been explicitly studied in the context of SC generation.

In our simulations, we generate a large ensemble of 200000 different realizations of SC broadening over 20 m propagation, each realization using identical parameters aside from the initial noise seed. A subset of 1000 individual output spectra from this ensemble is shown as the superposed gray curves in Fig. 1(a). The solid line is the calculated mean of these results, and we clearly see the presence of significant fluctuations between realizations. These intensity fluctuations are accompanied by shot-to-shot variation in the spectral phase at each wavelength, which can be readily quantified through the degree of mutual coherence $\left|g_{12}^{(1)}(\omega )\right|=|〈{\tilde{A}}_i^{*}(\omega ){\tilde{A}}_j(\omega )〉/〈{\left|{\tilde{A}}_i(\omega )\right|}^2〉{〈{\left|{\tilde{A}}_j(\omega )\right|}^2〉}^{1/2}|$where angle-brackets indicate ensemble averages over the spectral amplitudes at each wavelength and the indices i, j with i ≠ j indicate different realizations in the ensemble [21]. This regime of unstable SC generation is associated with essentially near zero mutual coherence ( $g_{12}^{(1)}<0.05$) across the spectral bandwidth shown in the upper sub-figure.

The statistical properties of the shot-to-shot fluctuations are shown in Figs 1(b)–1(d) where we filter the SC at different wavelength regions using a 20 nm bandwidth filter as shown, and determine histograms based on the energy of the corresponding temporal pulses [22]. These results illustrate the different nature of the SC fluctuations at different wavelength ranges. Near the spectral edges [(b) and (d)] we see highly skewed and long-tailed distributions, whereas near the center [(c)] we observe a narrower near-symmetric distribution. The insets plot the histograms on log scales to highlight the characteristics in the wings. The qualitative characteristics of these histograms can be complemented by standard statistical tests for distribution fits. We find for example, that the near-symmetric histogram in Fig. 1(c) can be fitted by a Gaussian distribution at the 0.05 significance level (chi-squared statistic), whereas no such meaningful Gaussian fit can be made to the long tailed distributions in Figs 1(b) and 1(d). On the other hand, the histograms in Figs 1(b) and 1(d) are well-fitted by the long tailed Frechet distribution at the 0.05 significance level (chi-squared statistic), and are indeed qualitatively very far from Gaussian distributions.

3. Random walks and Lévy flights

The intensity and phase fluctuations in SC generation are often represented in terms of results such as those in Fig. 1, as the spectral fluctuations and their probability distributions can be directly accessible via experiment. In this section, however, we consider a novel and alternative representation when we consider how the SC instabilities can be represented in a way that makes links with other classes of physical random process more apparent. This approach provides new insights into the nature of the intensity and phase fluctuations in SC generation.

The first process we discuss is the generic case of the random walk (Brownian motion) [23, 24]. In particular, we consider how the loss of spectral coherence in SC generation arising from spectral phase fluctuations can be used to construct an isotropic random walk in two dimensions. Using the simulation results above, the first step in this approach is to consider phase fluctuations at a wavelength where the coherence is near-zero, and verify that the phase is uniformly distributed across the ensemble in the range 0 – 2π. This is straightforward using standard statistical tests. For the results in Fig. 1, the SC phase distribution was found to be fitted by a uniform distribution at the 0.05 confidence level at all wavelengths where the mutual coherence satisfied $g_{12}^{(1)}<0.02$. Note here we consider the phase extracted at a particular wavelength at the resolution of the simulation discretization.

We show in Fig. 2(a) how these phase fluctuations can be used to construct a two dimensional random walk. Here, the randomly-varying phase from the different realizations k in the ensemble is used to determine the direction of a series of sequential unit length complex vectors. In particular, each SC realization yields a vector r_k = exp(iφ_k), and by combining different realizations we trace an n-step trajectory $\mathbf{r}(n)={\sum }_{k=1}^n{\mathbf{r}}_k$ in the complex plane. Typical results are shown in Fig. 2(a) where we construct a trajectory from 1000 realizations (i.e. a walk of length n = 1000 steps) with the phase properties extracted from the SC at 1100 nm. Note that the figure superposes the results of 20 such 1000-step trajectories to highlight how the uniformity of the SC phase distribution leads to directional isotropy over multiple walks. The evolution of any particular trajectory resembles qualitatively that expected of a random walk, with the inset illustrating clearly the random direction of different steps near the end point of a particular trace. We stress that there was no noise on the input field at the SC wavelength at which these phase fluctuations are extracted; their origin is from the nonlinear transfer of noise over a limited bandwidth near the pump to other wavelengths across the SC spectrum.

 

Fig. 2 (a) Representation of 20 walks of 1000 unit steps in the complex plane based only on phase fluctuations in the SC using a fixed imposed unit step length. The phase is extracted from the SC at 1100 nm. (b) Mean squared displacement (MSD) plotted as a function of step number calculated from an ensemble average of 200 realizations. Results for the phase extracted from three wavelengths are shown: 880 nm (black diamonds), 1100 nm (blue circles), 1270 nm (red squares). The black line shows the expected MSD 〈|r(n)|²〉 = n for an ideal random walk of unit steps.

Download Full Size | PDF

We can readily test these results quantitatively using standard theory. Specifically, we calculate the mean squared displacement (MSD) of a particular trajectory after n steps as MSD(n) = 〈|r(n)|²〉 where r(n) is the displacement from the origin after n steps. The results in Fig. 2(b) show the calculated MSD as a function of the number of steps n up to n = 1000 and where the ensemble average 〈 〉 is evaluated over 200 such trajectories (so that we use the full ensemble). We show results using the phase extracted from three different wavelengths in the SC (880 nm, 1100 nm, and 1270 nm), and we note that the results show essentially identical behavior. We can also readily verify that the MSD scales as expected for an ideal unit step random walk according to 〈|r(n)|²〉 = n; this expected result is shown as the solid black line in the figure. Significantly, the near-identical random walk behavior at different wavelengths highlights the independence of intensity and phase statistics; we see phase isotropy at all wavelengths even though the results in Fig. 1 show very different intensity statistics at the three wavelength ranges considered.

The results above combining the random phase properties of the SC with unit steps can be readily generalized to study a wider class of random walk process. In particular, motivated by important analogies with hydrodynamics, there has recently been much interest in the highly skewed probability distribution of SC amplitude noise at wavelengths near the long-wavelength edge of the spectrum; we show here that we can exploit this also to study more general classes of random walk. In particular, we use the properties of an optical SC in the regime of long-tailed intensity distribution to construct two dimensional trajectories, and we discuss how this leads to characteristics similar to the Lévy flight, an important class of random walk where steps are made in isotropic random directions but with step-lengths governed by a probability distribution that is heavy-tailed. Lévy flights and Lévy statistics are of increasing importance in understanding many classes of system in diverse fields such as physics and biology, and it is interesting to see how the same characteristics can also appear in the spectral intensity characteristics of an optical SC.

The presence of long tailed distributions near the spectral edges is seen clearly in Figs 1(b) and 1(d), arising from physical mechanisms associated with soliton dynamics that are particularly sensitive to input noise [25]. For our purposes, the wavelength ranges where such long tailed distributions are observed are those where we can consider extracting SC characteristics for the construction of Lévy flight-like trajectories. In contrast to the results in Fig. 2, the idea here is to use both the intensity and phase properties of the SC in constructing the random walk. That is, the long-tailed distributions of pulse intensity at particular wavelengths are used to determine the lengths of particular steps in the walk, whereas the corresponding direction at each step is determined by the phase taken from the same wavelength region.

We first show in Figs 3(a)–3(c) the general form of such trajectories obtained by constructing random walks from the SC characteristics at the three wavelengths shown in Fig. 1(a) 880 nm, 1(b) 1100 nm, and 1(c) 1270 nm. In each case, each realization is used to yield: (i) a step length determined from the energy calculated over a 20 nm bandwidth about the specified center wavelength, and (ii) a step direction determined from the phase at the wavelength in the center of the 20 nm band.

 

Fig. 3 For wavelengths of (a) 880 nm, (b) 1100 nm, and (c) 1270 nm, the upper panels show results of 20 walks of 1000 steps. The lower panels show the corresponding mean square displacement (MSD) plotted as a function of step number calculated from an ensemble average of 200 realizations.

Download Full Size | PDF

Note that the phase distributions were well-fitted by uniform distributions as described above, but the different intensity distributions in each wavelength region lead to very different characteristics for the associated random walks. Near the spectral edges, where the distribution is highly skewed [see Figs 1(b) and 1(d)], the trajectories are associated with rare long flight segments between turning points about clusters of shorter flight segments. This qualitative behavior, shown in Figs 3(a) and 3(c), is that which is typical of Lévy flights [26–30]. On the other hand, near the pump wavelength, Fig. 3(b) shows trajectories where the intensity distribution is approximately symmetric and Gaussian-distributed [see Fig. 1(c)]. In this case, the calculated MSD exhibits a linear dependence with step number with 〈|r(n)|²〉 = Dn (where we define D as a fitting constant). Note that the value of the constant D = 3.46 × 10⁸ agrees well with the calculated mean square step length 〈|r_k|²〉 = 3.68 × 10⁸ ∼ D.

A characteristic feature of the Lévy-flight like trajectories in Figs 3(a) and 3(c) is the observation of “superdiffusion” where the MSD varies with step number according to a more general power law of the form 〈|r(n)|²〉 = Dn^α with α > 1. The lower subfigures in Figs 3(a) and 3(c) show the evolution of MSD with n for these trajectories. Note that because of our limited size ensemble, fitting the observed MSD evolution with step length is not meaningful, but we show the slope of the strictly linear evolution with α = 1 for comparison.

An additional interesting feature of the trajectories constructed with step lengths associated with long-tailed distributions is the fact that they exhibit clustering in self-similar patterns characteristic of fractals. This is shown by taking particular trajectories from Figs 3(a) and 3(c) and examining their characteristics under successive zooms about turning and clustering points in the trajectory path. Figure 4 shows this behavior for trajectories constructed from SC characteristics at (a) 1270 nm and (b) 880 nm, clearly revealing the qualitative features of fractal scale-invariance.

 

Fig. 4 Plotting successive zooms about turning and clustering points as shown for trajectories at (a) 1270 nm and (b) 880 nm reveals the qualitative features of fractal scale invariance.

Download Full Size | PDF

4. Random number generation

The results above show that the intensity and phase fluctuations of supercontinuum generation indeed exhibit characteristics of random processes. An application that immediately suggests itself is in the generation of random numbers, a field that has received much attention recently in photonics as advances in information technology have highlighted the need for physical random number generators at optical wavelengths [10, 11, 14, 15]. In this section we use simulations to assess the ability of the optical SC to be applied as a physical random number generator for this purpose.

We consider simulation parameters similar to those used above: 3 ps pulses of 300 W peak power injected at 1064 nm and propagating in 10 m of fiber. Note here that we use a higher peak power and shorter fiber length allowing us to generate broad SC bandwidth at a significantly reduced computational cost; this is necessary in order to generate an ensemble of 10⁶ realizations in order to perform meaningful statistical tests as we describe below.

The underlying principle of random number generation from a noisy SC is to convert the intensity at a particular wavelength in the SC into either a 0 or a 1 depending on its value relative to a threshold, exploiting shot-to-shot fluctuations between different realizations in the ensemble to generate a random binary sequence. In practice, the different SC in the ensemble would be generated from an incident pulse train from a mode-locked laser which would determine the sequence repetition rate. For our proof-of-principle study, we use the same approach as above, seeding the SC with noise only over a finite range of wavelengths around the pump, and sampling the SC fluctuations at wavelengths outside this range. This is to distinguish the random characteristics of the sequence due to nonlinear propagation dynamics from the properties of the noise seed.

The manner in which random numbers are extracted from the SC ensemble is illustrated in Fig. 5. For each realization, we extract the intensity at a particular wavelength at the resolution of the simulation discretization, and the first step is to process a sub-sequence of realizations in order to determine a median to define a threshold value. A plot of the spectra generated under these conditions is shown in Fig. 5(a), and Fig. 5(b) plots a time series (as the gray points) of filtered intensities at 1140 nm obtained by plotting the intensity at this wavelength over 50000 realizations. The median of the time series (black line) is computed dynamically as a function of the number of realizations processed, and stabilizes rapidly as seen in the figure. Note that for a total sequence length of 10⁶, we found that determining the threshold over the first 50000 realizations yielded good results. After this initial step, intensity values from following realizations are converted to 1 or 0 respectively depending on whether they are above or below threshold. This is illustrated in Fig. 5(c) and allows us to generate a longer binary sequence (10⁶) with the required symmetric distribution of 0 and 1’s [31].

 

Fig. 5 Schematic showing random number number generation from spectral instabilities. (a) Results from 50000 simulations (gray) and calculated mean (black). The red line shows a particular sample wavelength of 1140 nm. (b) For this wavelength, we construct a time series from individual realizations in the ensemble and calculate a rolling median to determine a threshold over the first 50000 realizations. (c) Subsequent intensity values at this wavelength are compared to this threshold to yield a binary sequence.

Download Full Size | PDF

The degree of randomness of computed binary sequences can be readily tested using the standard statistical benchmark provided by the National Institute of Standards (NIST) [32], but a sequence length of only 10⁶ bits is insufficient to pass the NIST benchmarks with appropriate statistical significance. Nonetheless, we note that the results of the NIST test for the raw sequence of 10⁶ bits were consistent with the results obtained using the same sequence length of pseudo-random numbers generated from well-known algorithms.

Generating the required sequence length of > 10⁸ bits numerically for a meaningful NIST test is computationally prohibitive, but we note that it would actually be trivial in an experiment by using a high-repetition rate MHz or GHz source. On the other hand, we have found that we can generate a suitable longer sequence suitable for NIST testing from our simulation results by a spectral multiplexing step, concatenating sequences of 10⁶ bits generated at 200 different wavelengths across the SC extracted uniformly from the wavelength ranges shown as the shaded regions on both short and long wavelength ranges of the SC in Fig. 5(a).

We remove the possibility of any residual bias associated with binary conversion and wavelength correlation using a time-delay exclusive-or (XOR) operation with a delay of Δk = 50 realizations, and in this case the resulting sequence consistently passed all of the NIST statistical tests for randomness, as shown in Table 1 which lists the results of the NIST tests applied to 200 samples (i.e. wavelengths) of 10⁶ bit records obtained from the XORed sequence x[n] ⊕ x[n+50]. In order to pass each of the statistical tests, the composite P-value must exceed 10⁻⁴, and there may be no more than 7 failures out of 200 trials. (The random excursion variant test may have no more than 5 failures out of 120 trials.) The XORed data set passes all of the NIST statistical tests, and we also note that we obtain similar successful results for various delays in the XORed sequences (e.g. Δk = 10, 20, 100) and when using different number of realizations for establishing a comparative threshold (e.g. 10000).

Table 1. NIST benchmark tests results for 200 sequences of 10⁶ bits. For 200 sequences and a significance level α = 0.01, the P-value (uniformity of p-values) should be larger than 0.001 and a proportion of 193 test success for the whole benchmark (115 for random excursion (variant) tests) is required to succeed statistical tests. Note: In case of a test producing multiple result outputs, the worst case is shown.

View Table

5. Conclusion

Supercontinuum generation in optical fibres is a complex nonlinear dynamical process, and any noise present on the input field is well-known to generate significant fluctuations in the output spectra. In this paper, we have shown how these fluctuations can be interpreted in a novel fundamental way in terms of the characteristics of random walks, and we have shown an important applications potential of the supercontinuum as a physical random number generator. Our results also suggest new links with broader areas of optics and physics. In particular, for regimes of supercontinuum generation where long tailed intensity distributions are observed near the spectral edges, our results have shown how the statistics can be used to construct Lévy flight like processes. Since these long-tailed statistics correspond to the regime of rogue wave like behavior, our results suggest an important link between the dynamics underlying the generation of rogue waves and Lévy flights. Although we have studied this link in the specific case of an optical nonlinearity, we anticipate that it will hold generally for rogue waves in other systems.