Observation of rogue events in non-Markovian light

The efforts to understand the physics of rogue waves have motivated the study of mechanisms that produce rare, extreme events, often through analogous optical setups. As many studies have reported nonlinear generation mechanisms, recent work has explored whether optical rogue events can be produced in linear systems. Here we report the observation of linear rogue events with tunable height, generated from light imprinted with a non-Markovian wavefront. Moreover, if the non-Markovian wavefront is allowed to propagate through a nonlinear medium, extraordinarily long-tailed intensity distributions are produced, which do not conform to the statistics previously observed in optical rogue wave experiments.


Introduction
Over the past centuries, rare incidents of unexpected, extraordinarily large oceanic waves have been reported, together with their devastating consequences. 1 The effort to understand the mechanisms driving the generation of such events, referred to as rogue waves, has motivated the study of systems with similar behavior, particularly in optics. In a pioneering experiment, Solli et al. observed the generation of optical rogue waves in nonlinear media with modulation instability. 2 This demonstration opened the door to much research on the generation of optical rogue waves in nonlinear media. Optical rogue waves generated with modulation instability were modeled 3,4 and observed in several systems, such as optical fibers, 5 supercontinuum generation, 6,7 Raman fiber amplifiers, 8 photorefractive crystals, 9 nonlinear optical cavities, 10 and femtosecond filamentation. 11,12 Furthermore, the interaction of nonlinearity and turbulence 13,14 was studied with random light sources, 15 Kerr resonators, 16 and ring-cavity lasers. 17 Optical rogue waves have been demonstrated also in dissipative nonlinear systems such as femtosecond lasers 18,19 and fiber lasers. 20 In light of these results, the question arose as to what role does nonlinearity play in generating rare, extreme events and whether it is a fundamental re- †These authors contributed equally to this work. *Corresponding author: hfrostig@bu.edu. quirement for their existence. [21][22][23] Deviations from Rayleigh statistics, a criteria that the optical community often uses for the existence of rogue waves, can be observed when the conditions of Rayleigh scattering are violated. 24 Intensity distributions obeying Rayleigh statistics are generated when a homogeneous wave is Rayleigh scattered by a large number of random scatterers, creating a uniform distribution of phases in the interval [−π, π] in the scattered wave, and the scattered wave is observed at a distance from the scattering layer (fully-developed speckles). 25 The distribution of intensities in the scattered wave is a decaying exponential, conventionally plotted as a linear distribution in semilog scale. By devising systems that defy the conditions of Rayleigh scattering, long-tailed intensity distributions have been theoretically predicted 26,27 and observed in various systems without the use of nonlinearity, including microwave resonators, 28 nanowire mats, 29 partially developed speckle fields, 30 inhomogeneously illuminated scatterers 31 and speckle fields with strong phase modulations. 32,33 Particularly, wavefront shaping provides a means of generating non-uniform phase distributions 34 in a controllable manner, and has been used to generate long-tailed distributions by numerically retrieving the phase mask necessary to create them. 35 Here we study an optical wave model that uses uniform illumination of a large number of random scatterers, yet does not lead to a Rayleigh distri-bution in the far-field. The generation mechanism we study is non-Markovian, 36,37 inspired by the non-Markovian behavior of dynamical turbulent systems in the oceanic environment, such as sea surface winds. [38][39][40] Non-Markovian distributions are distributions with long-range correlations and hence some degree of memory. We show that a non-Markovian light source creates rare, extreme events even in the linear regime, yet their intensities are drastically enhanced by nonlinearity in the medium through which they propagate. Our system allows for tunability of the height of the rogue events, which can be used to study their behavior in different regimes, including a regime with statistics resembling generalized extreme value (GEV) statistics.

Principle
To generate non-Markovian light, we imprinted a coherent wavefront with a two-dimensional spatial phase mask containing random phases with a tunable amount of long-range order. A number set is said to be non-Markovian when the value at a given position within the set is statistically dependent on numbers at other positions, in addition to its nearest neighbors on both sides. Such a set would result, for example, from drawing distinct balls from a sac without replacing them, so that the outcome of each draw depends on all previous draws. In order to produce a two-dimensional phase pattern that would be non-Markovian on the one hand, and uniformly distributed within the range [−π, π] on the other hand, we followed the procedure described by Fischer et al. 36 and illustrated in Fig. 1. In short, the phase values are chosen by solving overlapping Sudoku puzzles. Each puzzle is a 9x9 square in which the numbers 1-9 must appear only once in each row and each column. After the values of a single square have been assigned, the next square is defined such that some of its columns (or rows) overlap with the previous square. We denote the number of overlapping columns (or rows) by r. The larger r is, the more memory there is in the generated pattern and consequently the larger the degree of non-Markovianity. Note that r = 9 -s, where s is the shift parameter defined by Fischer et al.. This procedure is repeated until the entire matrix is filled, resulting in a random correlated number set. 41,42 This set of numbers, N, is then translated into phase values, φ, according to φ = 2π N 9 . The inset of Fig. 1 shows the probability distribution of the phases of the r=7 mask, which has the largest degree of long-range correlation. We see that adding long-range order to the random phase masks indeed does not modify the uniform distribution of their phases in the range [−π, π].
Our experimental setup is outlined in Fig. 2 be- (a) A one-dimensional non-Markovian phase mask is generated by solving overlapping one-dimensional Sudoku puzzles of length 9. In each such puzzle the number 1-9 should appear only once. (b) Two-dimensional non-Markovian phase pattern are similarly generated by solving overlapping twodimensional Sudoku puzzles of length 9. In each such puzzle the number 1-9 should appear only once in each row and each column. The cases shown in a and b are for an overlap of r=2. Inset: the probability distribution of the phases in the r=7 non-Markovian phase mask. The phase values are uniformly distributed in the range [−π, π] despite the non-Markovian properties.
low. The generated non-Markovian phase mask is imprinted onto light from a continuous-wave laser with a two-dimensional spatial light modulator. The shaped wavefront undergoes Fourier transform by a lens and then propagates through a photorefractive crystal (SBN:60, see methods). Since the nonlinear refractive index of the crystal depends on the strength and sign of the electric field applied to it, the crystal is used as either a linear, positive nonlinear, or negative nonlinear medium. The output facet of the crystal is imaged onto a CCD camera and the recorded intensity distribution is analyzed. The experiment is repeated for several overlap amounts, namely r=0,4,7.

Experimental results
For initial evaluation of the experimental results, we plotted 1D slices of the recorded 2D output intensity patterns for different conditions. Fig. 3 compares 1D Figure 2: Generating rogue events with a non-Markovian light source. A 532nm laser beam is expanded using a telescope to illuminate a two-dimensional SLM. The SLM is used to imprint a non-Markovian spatial phase onto the beam. The beam is then focused into a photorefractive crystal. The output facet of the crystal is imaged onto a CCD camera.
slices taken at the location of the maximum intensity of the output patterns for the cases of propagation through a linear (blue), positive nonlinear (red), and negative nonlinear (green) media, for various overlap values. The r=0 case (Fig. 3a) represents the non-Markovian phase mask with the least amount of order, since the phase of a given pixel is statistically dependent only on the other pixels within its own 9x9 square. In contrast, the r=7 case ( Fig. 3c) represents the non-Markovian phase mask with the most order, since the large overlap between the squares places the largest number of constraints on the value of a given pixel possible, while maintaining random assignment. An established criterion for rogue waves is whether their height exceeds two times the significant wave height (SWH), defined as the mean amplitude of the highest third of the wave events. 22,43 Dashed lines at two times the SWH value of each distribution are plotted in the corresponding colors for reference. We see that the combination of positive nonlinearity and non-Markovianity yields extreme events with an intensity that can be tuned with the long-range order parameter, r. For r=7 and positive nonlinearity, the extreme events produced exceed 20 times the SWH.
The intensity probability distributions of the experimental results are presented in Fig. 4. Figures 4a-c compare the probability distributions of the light intensity at the crystal output for positive nonlinearity, negative nonlinearity, and no nonlinearity, for varying overlap values. All probability distributions are normalized by the maximum probability of their distribution. Dashed vertical lines indicate the SWH of the distributions. For no overlap (r=0), the linear probability distribution (blue data in Fig. 4a) is approximately a linear line in semilog scale, indicating that it is quite similar to a random, non-correlated field. As expected, when this beam propagates under the influence of positive nonlinearity, the probability distribution develops a long tail, indicating the existence of high intensity, rare events (red data in Fig. 4a). [44][45][46] This distribution obeys super-Rayleigh statistics, similar to the statistics generated using modulation instability in a positive nonlinear medium. Conversely, when the r=0 imprinted beam propagates under the influence of negative nonlinearity, the probability for high intensity events is suppressed and the probability for low intensity events rises, obeying sub-Rayleigh statistics (green data in Fig. 4a). Similar relations can be seen between the blue, red and green data in figures 4b and c.
However, as the long-range order in the imprinted phase pattern is increased, super-Rayleigh statistics develops even without nonlinear propagation. This can be seen in the blue data in figures 4a-c, and becomes more apparent when comparing the linear intensity probability distribution for the r=0 (purple) and r=7 (turquoise) patterns directly, as in Fig. 4d. The distributions were fitted with a stretched exponential of the form P pow (I) = exp(−aI b + c), where I is a vector of intensities, P(I) is the normalized probability to detect a specific intensity, and a, b, and c are parameters. While the logarithm of the linear r=0 distribution fits a power-law with exponent 1, the logarithm of the linear r=7 distribution fits a power-law with exponent 0.78, corresponding to super-Rayleigh statistics.
Though non-Markovianity alone serves to generate extreme events, their intensities are not substantially larger than the intensities observed at the end of the tail of the intensity probability distribution of the r=0 imprinted beam. Yet if a medium with positive nonlinearity is 'seeded' with non-Markovian light, the resulting probability distribution has an extremely long tail, longer than that created by either effect on its own. This can be observed in Fig. 4e, which presents a comparison of the distributions obtained from non-Markovian light with no overlap (r=0) and large overlap (r=7) after propagation in a medium with positive nonlinearity. Though both distributions are super-Rayleigh, the distribution obtained for r=7 (turquoise) follows a more extreme trend than that obtained for r=0 (purple). In fact, the distribution obtained for r=7 after nonlinear propagation does not seem to converge to the power-law distribution P pow , typical of rogue waves generated via modulation instability in nonlinear media. 9,31,44 Instead, the distribution is best described by an exponential function of the form P exp (I) = exp[a exp(−bI) − c] which resembles the Gumbel distribution, 47 a particular case of the generalized extreme value (GEV) distributions. Both fits are shown in Fig. 4e.
Overall, the extreme events generated by the com- bination of non-Markovianity and positive nonlinearity are ∼ 2.5 times more intense than those created just by nonlinearity (best seen when comparing the maximum of the red data in Fig. 3c to the maximum of the red data in Fig. 3a), and ∼ 3 times more intense than those created just by non-Markovianity (red data in Fig. 3c vs. blue data in Fig. 3c). Furthermore, the enhancement of nonlinearity and non-Markovianity combined, N N mN l , is about a factor 5 (red data in Fig. 3c vs. the blue data in Fig. 3a), whereas the two individual enhancements are N N m 1.5 for non-Markovianity alone (blue data in Fig. 3c vs. the blue data in Fig. 3a) and N N l 2 for nonlinearity alone (red data in Fig. 3a vs. the blue data in Fig. 3a). Therefore N N mN l > N N m * N N l , showing that the two effects are synergistic and their interplay serves to radically increase the wave amplitudes.

Simulation results
In order to verify the experimental results, we performed numerical simulations of non-Markovian light propagation through nonlinear media. When the nonlinear medium is a photorefractive crystal, the propagation is described by a modified version of the nonlinear Schrödinger equation: 48 (1) Where k 0 = 2π λ and λ is the wavelength in vacuum, n 0 is the linear refractive index of the medium, and γ is the voltage dependent nonlinear coefficient of the crystal. I bkg represents the intensity of a second beam incident upon the crystal, of homogeneous incoherent light, which serves as background illumination (see methods). While I bkg was kept constant in the experiment and simulations, the voltage on the crystal was varied such that γ > 0, γ < 0 or γ = 0. To simplify Eq. 1, the last term can be approximated as −γE 1 − |E| 2 /I bkg + |E| 4 /I 2 bkg for |E| 2 I bkg . For positive nonlinearity (γ > 0), this expression includes the well-known Kerr selffocusing contribution as well as a higher-order defocusing contribution which prevents the collapse of small features in E(x, y) to sizes below the resolution of the simulation. For negative nonlinearity (γ < 0), the term γ|E| 2 E/I bkg causes self-defocusing and the physics is well described with only the third-order term, so the higher order term, −γ|E| 4 E/I 2 bkg , was neglected. The intensity distributions after propagation through the medium were analyzed for non-Markovian phase patterns with r=0,4,7 and the corresponding probability distributions are shown in Fig. 5. Again, Fig. 5a-c show a comparison between the distributions resulting from propagation in a lin-ear (blue), positive nonlinear (red) and negative nonlinear (green) medium, for r=0,4,7. Dashed vertical lines in the corresponding colors are plotted at two times the SWH of each distribution. Fig. 5d-e show a comparison between the intensity probability distributions for the r=0 and r=7 imprinted beams, after linear and positive nonlinear propagation, respectively. We see that the simulation results qualitatively agree with the experimental results, showing distributions with increasingly long tails for increasing r values. Moreover, the combination of non-Markovian light and positive nonlinearity results in an extraordinary-long tailed distribution, which is best fit with an exponential of an exponential, as in the experimental results (the fit to a power law is also shown). We note that the effects of nonlinearity on the distributions are stronger in the simulated data than in the experimental data, since we set γ/I bkg to a slightly higher value in the simulation in order to make the trends more clear.

1D simulation
To gain some insight on the reason that non-Markovianity generates long-tailed probability distributions, we take a closer look at the output intensity patterns. In order to make the effect visually clear we use a 1D simulation, which allows for a significantly larger number of statistics compared with a 1D slice of a 2D simulation. The 1D simulation follows the same physics described in section 3.2, yet in one dimension (masks are generated as depicted in Fig. 1a). The results for a completely random phase mask, as well as for r = 0,4,7,8 after linear (blue), positive nonlinear (red), and negative nonlinear (green) propagation are presented in Fig. 6 below. For clarity we have replaced the x-axis, that represented CCD pixles in our previous plots, with the frequencies (denoted as k) that represent the Fourier-conjugate variable of the pixels of the SLM. Fig. 6a presents the case of a completely random phase mask, where the numbers 1-9 are randomly assigned to pixels on the SLM. The far-field pattern is a Rayleigh-distributed speckle pattern for linear propagation, super-Rayleigh distributed speckle pattern for positive nonlinear propagation, and sub-Rayleigh distributed speckles for negative nonlinear propagation, as has been demonstrated previously 44,49 (for plots of the corresponding probability distributions see Fig. S1).
As in the 2D case, the r=0 case represents the non-Markovian phase mask with the least amount of order, since the phase of a given pixel is statistically dependent only on the other pixels within its own 1x9 puzzle. While the far-field is still close to a random speckle pattern (see Fig. 6b), we can gain some understanding of its nature by considering how of- ten numbers can repeat in the near-field phase mask. The constraint that the numbers 1-9 appear once in each puzzle suppresses the probability of very low frequencies, since numbers must repeat at least every 17 entries. This manifests in a dip in the probability of high intensity speckles near k = 0 (the fall off is slow, rather than sharp, due to the random pattern superposed). This attenuation drives some energy into the rest of the spectrum, causing some higher intensity speckles to appear and slightly elongating the probability distributions. Nevertheless the distributions do not significantly deviate from those of a completely random phase mask (see Fig. S1a and. S1b).
As the value of r increases, more long-range order appears in the applied phase masks, with r=4 being an intermediate case (Fig. 6c). The dip in the low frequencies widens, pushing more energy into the speckles in the intermediate frequency region. While the effect is less pronounced for the linear case, it is significant for the positive nonlinear case. As a result, the probability distributions, in particular the positive nonlinear distribution, become appreciably elongated.
The most interesting case is r=7, shown in Fig. 6d. This mask is the non-Markovian phase mask with the most order, yet random assignment is maintained. Since numbers must repeat every 8-10 entries, the most likely frequencies are k = 1/9 and adjacent frequencies, as well as the second and third harmonic of these frequencies. As a result, the output pattern has three regions in which the likelihood of high-intensity speckles is higher, centered around k = 1/9, k = 2/9 and k = 3/9. The region centered at the fundamental frequency, k = 1/9, is the most narrow and sharpest, while the harmonic regions become gradually wider and less sharp. In between the highprobability regions there are dips in the pattern, or regions of low-probability of high-intensity speckles. This pattern of low-and high-probability regions is most pronounced for the positive nonlinearity case, where it contributes to the emergence of an extraordinarily long probability distribution (see Fig. S1d), but is present also for the linear and negative nonlinear cases.
Finally, the r=8 case represents a completely ordered phase mask. Order emerges because when assigning the ninth pixel of a new puzzle whose previous eight pixels overlap with an existing puzzle, only one 'free' number out of the numbers 1-9 is left. Therefore, the first puzzle is randomly assigned and it is repeated periodically thereafter until the whole phase mask is assigned. The near-field pattern can thus be described as a convolution of a Dirac comb and a random phase pattern of a length of 9 entries.
As a result, the far-field pattern after linear propagation, shown in Fig. 6e, is a Dirac comb multiplied by a random intensity pattern, resulting in a comb with random teeth height. The width of the teeth is then broadened by diffraction during linear propagation (the patterns after nonlinear propagation are not presented for this case for clarity). We note that for r=8 the probability distribution is again a Rayleigh distribution (see Fig. S1e). Therefore the more longrange order exists in the non-Markovian mask, the more elongated the intensity distribution. Yet if the random assignment disappears completely, so does the long-tailed intensity distribution. (a) For the r=8 phase mask. As order is increased, regions of high probability for high intensity events emerge, causing the elongation of the intensity distribution (see Fig. S1). Nevertheless when random assignment disappears completely, as in the r=8 phase mask, the Rayleigh intensity distribution reemerges. The nonlinear cases are not shown for r=8 for clarity.

Discussion
In this work we have observed, for the first time, the generation of extreme, rare events from a non-Markovian light source. Moreover, we have shown that when a nonlinear medium is seeded with non-Markovian light, the interplay between these two effects leads to the generation of exceptionally large rogue events. The resulting intensity distribution follows a more extreme trend than distributions conventionally generated in rogue wave experiments, such as those obtained with nonlinearity alone.
The unique statistical properties of our system, as well as its tunability and experimental simplicity, make it a convenient test-bed for studying the physics of extreme events. Furthermore, as the probability distributions created by our system resemble generalized extreme value distributions, it may be a suitable model for extreme-value systems that could not be modeled previously.

Experimental implementation
Light from a 532nm CW laser (Coherent Verdi) was magnified with a 4x telescope, and reflected off a two-dimensional liquid-crystal spatial-light modulator (Hamamatsu LCOS-SLM x10468). The SLM was used to imprint the beam with a non-Markovian phase mask.
The outgoing field was Fouriertransformed using a 500mm focal length lens in order to obtain a non-Markovian intensity pattern (fullydeveloped speckles), and inserted into a photorefractive crystal (SBN:60) with a 5mm x 5mm crosssection and 10mm length. The output facet of the crystal was imaged onto a CCD camera (UI-1250 IDS) and the recorded intensity patterns were analyzed. The background illumination for the crystal (represented by I bkg in Eq. 1) was generated by splitting the main beam before the telescope, and passing it through a rotating diffuser in order to destroy its coherence. The intensity of the background illumination was ∼4 times the intensity of the non-Markovian light. The SLM pixels were binned into 5-by-5 macro pixels. Each phase mask was displayed on the SLM and the intensity probability distribution was computed from the histogram of the pixel values in the corresponding CCD image, with 50 equallyspaced bins. In distributions where there were less than 50 discrete pixel values (such as the negative nonlinearity case for low r values), the maximum possible number of bins was used. Each distribution was normalized to its maximum value. Pixels in the image with a value of zero were eliminated in order to use a log plot. For each r value, different realizations of the non-Markovian phase mask were generated by computing different solutions of the Sudoku puzzles. The measurements with linear propagation were performed without applying voltage to the crystal, and were averaged over 20 real-izations of the non-Markovian phase mask for each r value. The voltage applied for positive nonlinearity measurements was ∼0.7kV, and the results were averaged over four realizations, with 50 measurements for each realization. The voltage applied for negative nonlinearity measurements was -0.5kV, and the results were averaged over two realizations. In all three cases, the number of realizations was chosen to maintain reasonable experimental run times after verifying numerically and experimentally that the amount of statistics produced created smooth intensity distributions.

2D simulation implementation
The non-Markovian phase masks used were either the same as those used in the experiment, or generated in the same way. Nonlinear propagation was simulated using the split-step method. Time-domain was not addressed specifically in the simulation, as this would result in impractical run times. Defocusing nonlinearity was modeled as Kerr-type, with the nonlinear operatorN = −γ/I bkg |E| 2 . Focusing nonlinearity was modeled with a higher-order defocusing term,N = γ/I bkg (|E| 2 − 1/I bkg |E| 4 ). The higher order term introduces saturation of the self-focusing process, as occurs in practice in photorefractive crystals, 50 and prevents the collapse of the speckles to sizes below the resolution of the simulation. The results were averaged over 20 realizations for linear, positive nonlinear, and negative nonlinear propagation. All histograms were computed using 50 equally spaced bins and were normalized to their maximum values. Zeros values were removed due the log plot.

1D simulation implementation
The simulation implementation followed that of the 2D simulation, but in a single spatial dimension ((1+1)D geometry). The results were averaged over 100 realizations for the completely random case and for r=0,4,7. For r=8 the results were averaged over 1000 realizations since the r=8 output pattern contains significantly less statistics than other r values. All histograms were computed using 25 equally spaced bins and were normalized to their maximum values. Zeros values were removed due the log plot.

February 28, 2020
This document provides supplementary information to "Observation of rogue events in non-Markovian light". The supplementary material includes two sections. The first section presents the intensity probability distributions of the one-dimensional simulated data presented in Fig. 6 of the main text. The second section presents the propagation traces of the one-dimensional simulated data through the linear and nonlinear media.
1 Intensity probability distributions of the onedimensional simulated data Below we present the intensity probability distributions computed from the 1D simulation, corresponding to Fig. 6 of the main text. Figures S1a-e show the distributions generated by the totally random phase mask, the r=0 mask, the r=4 mask, the r=7 mask, and the r=8 mask, respectively. Each sub figure displays the distribution after propagation through a linear (blue), positive nonlinear (red) and negative nonlinear (green) medium. As in the 2D simulation, the distribution becomes more elongated as r increases, and the combination of positive nonlinear propagation and large r values serves to create extremely long-tailed distributions. Yet when the order is increased further and the phase mask becomes completely ordered, as occurs for r=8 (Fig S1e), the distribution becomes linear just like in the completely random case (Fig S1a). †These authors contributed equally to this work. *Corresponding author: hfrostig@bu.edu. Examples of waves exceeding two times the SWH are marked with white arrows. We note that the generated rogue events for all three mediums are transient.