Characterizing the Statistics of a Bunch of Optical Pulses Using a Nonlinear Optical Loop Mirror

We propose in this work a technique for determining the amplitude distribution of a wave packet containing a large number of short optical pulses with different amplitudes. The technique takes advantage of the fast response of the optical Kerr effect in a fiber nonlinear optical loop mirror (NOLM). Under some assumptions, the statistics of the pulses can be determined from the energy transfer characteristic of the packet through the NOLM, which can be measured with a low-frequency detection setup.The statistical distribution is retrieved numerically by approximating the solution of a system of nonlinear algebraic equations using the least squares method. The technique is demonstrated numerically in the case of a packet of solitons.


Introduction
Considering the ever-growing number of applications involving short and ultrashort optical pulses (sensing, micromachining, medicine, etc.), their precise characterization is an essential task.When the pulses are so short that their direct measurement by optoelectronic means becomes extremely costly or impossible (when the required bandwidth lies beyond that of the fastest photodetectors), all-optical techniques, like the widespread optical autocorrelation [1], are available.The optical autocorrelation technique relies on an ultrafast nonlinear process (typically, second-harmonic generation in a crystal) to measure the autocorrelation function of the optical pulse intensity profile.Although this technique is able to estimate the duration of pulses down to the few fs range, it provides little additional information.The autocorrelation however has been adapted and extended and several techniques were developed that allow complete pulse profile characterization (amplitude and phase of the electric field) [2][3][4][5][6].One noticeable example is the FROG (Frequency Resolved Optical Gating) technique [2], which combines optical autocorrelation with spectral analysis.However, complex data processing is usually required to retrieve the pulse profile from the measured data.Besides, the free-space implementation of such devices and the phasematching requirement (if the nonlinear process is the secondharmonic generation) impose careful beam alignment and crystal orientation and cause enhanced sensitivity to perturbations.Moreover, the mechanically imposed scanning range limits the flexibility of the measurement span.Linear techniques were also developed for the complete characterization of ultrashort pulses, being usually simpler in processing and requiring much lower intensities than their nonlinear counterparts [7][8][9][10][11][12].The drawbacks, however, include the need for complex electronics or high-speed detection, or for the precise characterization of a reference pulse (using a nonlinear technique).
Aside from the precise measurement of single or periodic identical pulses, like the stable trains of periodic solitons produced by mode-locked lasers, or relatively simple multipulse compounds, there is also a need for characterization of by far more complex waveforms consisting of a large number of pulses having different amplitudes, possibly varying over time.One typical example is the "sea of solitons" produced by the breakup of wide, smooth ns pulses in a nonlinear and dispersive fiber [13].In the frame of supercontinuum generation under ns pumping, pulse breakup into a large number of solitons constitutes the initial stage of the process [14].Therefore, an experimental method to estimate the statistics of the formed solitons would be helpful to understand the physics of supercontinuum generation in the ns regime, when numerical calculations tend to be prohibitively long.Another example is offered by the so-called noiselike pulses [15][16][17][18][19][20][21][22][23][24].Noiselike pulses are long (∼ns) packets containing from thousands to millions of sub-ps subpulses with randomly distributed amplitudes.Their potential applications include supercontinuum generation [25,26], nonlinear frequency conversion [27,28], sensing [29], and micromachining [30].Although noiselike pulsing is gaining recognition as a universal mode of operation of passively mode-locked fiber lasers [24], it differs from the more conventional (conservative and dissipative) soliton regimes in that the inner details of the noiselike pulse (i.e., the subpulses) vary widely after each round-trip.Considering the very large number of subpulses and their variability, a statistical approach seems adequate to tackle these complex objects.Such an approach could thus be considered to reinforce the few recent efforts [31] addressing the challenging task of developing novel techniques for noiselike pulse characterization.Once it will become possible to characterize precisely noiselike pulses using such techniques, the path will be open towards a better understanding of the physics underlying their formation and dynamics, a subject that has fueled controversy over the years [15][16][17][18][19].
In this work we propose the use of a fiber nonlinear optical loop mirror (NOLM) [32] to determine the statistical distribution of the amplitudes of pulses in a bunch, from lowfrequency energy measurements.The NOLM is a versatile device that has been widely used for ultrafast switching, signal processing, and laser mode-locking applications, among others [33][34][35].The ultrafast (∼fs) response time of the optical Kerr effect in silica fiber is exploited, which requires no phase matching.An all-fiber solution avoids alignment issues, ensures a high robustness, and reduces costs, whereas the Sagnac architecture (in which both beams travel along the same path) improves the stability, especially under strict polarization control.

Design of the Setup and General Principle of the Method
A NOLM displays a power-dependent transmission (switching) characteristic, which is a sinusoidal function of input power.In most conventional, power-imbalanced schemes, in which an asymmetric coupler is used [32], destructive interference between the counter-propagating beams is not total, and therefore minimal transmission is not zero.Although zero minimal transmission can be achieved in a powerimbalanced scheme (in particular, using an attenuator located asymmetrically in the loop [34]), another option is offered by the polarization-imbalanced NOLM [36].In this scheme, a 50/50 coupler is used, so that the counter-propagating  powers are equal, and only the polarization symmetry is broken through the use of a quarter-wave retarder (QWR) in the loop.NOLM switching is then obtained thanks to the polarization-dependent Kerr effect [37].Because the device relies heavily on polarization, light polarization at the NOLM input should be strictly controlled, as well as the polarization of the counter-propagating beams in the loop; the latter condition can be met in practice if the fiber is twisted (twist makes the fiber optically active, so that light ellipticity is maintained during propagation) [38].On the other hand, polarization control offers enhanced flexibility of the switching characteristic.Flexibility can be further enhanced if a polarization selection is performed at the NOLM output [39].For example, assuming circular (say, right) input polarization, and selecting the orthogonal (left) polarization component at the NOLM output (see Figure 1), the NOLM power transmission is given by [40] where  in ,  out , and   are the input, output, and NOLM switching powers, respectively (  = 6/(), where  is the nonlinear coefficient for linear polarization and  is the loop length), and the phase shift Δ is related to the QWR angle through Δ = 4 (we chose to reference the QWR angle from the orientation for which Δ = 0).Equation (1) was derived using a continuous-wave approach [36], which disregards phenomena like group velocity dispersion and the group velocity mismatch between the two circular polarization eigenmodes in the twisted fiber, whose effects increase as the pulse duration reduces.Some highly nonlinear optical fibers are available, however, which are especially designed to have zero second-order dispersion and a low value of the dispersion slope at a specific wavelength (in particular at 1550 nm).On the other hand, the twist-induced group velocity mismatch is relatively moderate: only ∼0.3 fs/m at 1550 nm in a fiber twisted at a rate of 5 turns/m [41].In these conditions, (1) remains valid for pulses as short as 1 ps and below, even if they are chirped.
Although a transmission function such as (1) can be obtained using other schemes than the one depicted in Figure 1 (including power-imbalanced schemes), the linear dependence of Δ on the single parameter  in the proposed setup strongly facilitates the adjustment procedure of the transmission function.Besides, strict polarization control ensures that the transmission function will not drift over the total measurement time, like it usually happens with conventional schemes under slow (e.g., thermally driven) variations of the fiber residual birefringence (proper operation of the technique will suppose, naturally, that the transmission function is precisely known and does not vary over time).Finally, the requirement of a particular input polarization state may limit the range of applicability of the technique (e.g., solitons in a bunch may be differently polarized); on the other hand, the choice of circular input polarization makes sense considering the recent discovery that, in a twisted fiber, solitons tend to be circularly polarized [42,43].
The method developed here is an extension of a previously proposed method for determining the temporal profile of a short optical pulse through measurements of its energy transmitted through a NOLM with known power transmission characteristic [44,45].Because the pattern of energy transmission through the NOLM is pulse-shape-dependent, it contains information that can be exploited to retrieve the pulse profile.Simply stated, the principle of the method is as follows.We assume a periodic train of identical short pulses whose temporal profile is unknown.The energy of the pulses injected in the NOLM is varied (which can be done using a variable attenuator or amplifier), and, for each value of input energy, the pulse energy at the NOLM output is measured, which can be done using low-frequency electronics.An energy transfer characteristic is thus obtained this way.If now the NOLM transmission in power given by (1) (or equivalently, the instantaneous power transfer characteristic,  out versus  in ) is known, it can be used to calculate the energy transfer function for an arbitrary number of profiles.For proper NOLM adjustments, a different energy transfer function is obtained for each profile.The profile of the pulses under study is then determined by comparing the measured energy transfer characteristic with the calculated ones and choosing the profile whose curve is closest to the measured one.In practice, the unknown temporal pulse profile is approximated as being composed of a few tens of rectangular sections or "slices," whose amplitude approaches the pulse amplitude across the profile.Hence, instead of considering a continuous pulse profile in function of time, the problem is made discrete, and only a finite number of amplitudes have to be found to estimate the profile.The solution can thus be obtained by resolving a system of algebraic equations.Whereas the method is able to determine the amplitude (and even the duration) of the slices, because in nature it is based on the sole measurement of the overall pulse energy transmitted through the NOLM, it gives no information about the relative temporal position of the slices, that is, in which order the slices should be arranged so as to reproduce the searched pulse profile.Therefore, pulse profile retrieval is only possible if some assumption is made on the profile.For example, if a monotonously decreasing profile is assumed, the slices should be arranged from the highest to the lowest.
In the present case, where the signal to be characterized is a bunch of pulses with different amplitudes, each slice will represent one pulse or, more generally, a subcategory of pulses in the bunch having very close values of amplitude.As we are only interested in the statistics of the pulses, the temporal order of the slices is irrelevant and the technique is naturally suited for this purpose.Figure 2 illustrates the correspondence between the problem of determining the statistics of a bunch of pulses and determining a monotonously decreasing amplitude profile, taking as examples the Gaussian and uniform distributions.The figure also illustrates how the statistics of the pulses in a bunch affects the energy of the bunch at the NOLM output; this information can then be exploited to determine the statistics.In summary, the problem of determining the statistics of a bunch of pulses is completely equivalent to that of pulse profile retrieval and can thus be resolved using a similar approach, with the additional advantage that no assumption on the order of slices is needed.

Development of the Method for a Bunch of Rectangular Pulses
In a first approximation, we will assume that the bunch can be viewed as a series of a large number  of rectangles with equal duration and different amplitudes.Similarly to the case of pulse profile characterization, the technique relies on the observation that the NOLM energy transfer characteristic depends on the bunch statistics.As shown in [45], this dependence is stronger when the phase bias Δ in (1) is slightly negative (i.e., when the power transmission characteristic presents a minimum for some nonzero value of input power <   ).This is illustrated in Figure 3, where the transmission in energy of the bunch is plotted in function of the average power of the slices for a few common statistical distributions.The transmission in energy is given by where   is the amplitude of the th slice at the NOLM input and  is the total number of slices.
The technique assumes that we are able to measure   , or equivalently the energy transmitted by the NOLM, for several values of the input energy of the bunch (which can be varied, e.g., using a variable attenuation or gain).Such measurements can be performed using a low-frequency detection setup (low-frequency photodetector and oscilloscope).In most practical implementations, the different energy measurements will be performed sequentially, which supposes that the bunch can be reproduced periodically, conserving at least the same statistics at each period.Let us call   the attenuation/amplification coefficient for the th measurement and   the corresponding measured value of output bunch energy, for  = 1 to  ( being the total number of measurements).We will now model the bunch as consisting of  rectangular "slices," each of which regroups rectangles  c, f) illustration of the basic principle of the method: for the same energy (shaded area) of the input bunch, the bunch energy at the NOLM output is very different for each distribution (in this example, it is higher in the uniform case).
from the bunch which have similar amplitudes (each slice regroups the same number / of rectangles).The amplitude   of each slice is thus an approximation of the amplitudes of the rectangles in that slice.We then have that where Δ =   / is the duration of each slice (  is the total duration of the bunch), which is supposed to be known.Equation (3) is a system of  nonlinear algebraic equations in the  variables   .If  ≥  (i.e., if the number of energy measurements is at least equal to the number of slices), system (3) can be solved numerically to find the values of   .Note that, if   is unknown, the minimum number of equations is increased by 1; that is,  ≥  + 1; see [44,45].System (3) can be conveniently solved numerically by the least squares method, by minimizing the expression In order to keep computational time within reasonable limits, in practice the number  of variables (slices) should not be higher than a few tens, which is usually insufficient to assess with some precision the statistical distribution.There are a few ways, however, to increase the volume of the data.
First, more than one set of  energy measurements can be recorded.For example, if the setup is automated (automatic sweep of the NOLM input energy and automatic output energy measurement), the number  of measured values of energy   can be readily increased, say an integer number  of times, totaling  values that are then organized into  sets of  values.System (3) is then solved  times (once for each set of energy values), yielding  sets of solutions   .Besides, taking advantage of the sensitivity of the result to the initial estimate, to the degree of overdetermination of the system ( − ), and so forth, different sets of solutions   can be obtained even if only one set of measured output energies is available.
The proposed technique will now be tested numerically for different statistics of the amplitudes of the rectangles.The NOLM power transmission is given by (1), with   = 1 and Δ = −/5.First, the output energies (which, in a real situation, would be gathered experimentally) are computed using (2).In each case, the bunch consists in a large number  = 1000 of rectangles, whose amplitudes are random numbers generated according to a specific distribution.The values of   are increased until we get slightly beyond the minimum of the  curves shown in Figure 3 (for the technique to work properly, only a fraction of the switching power   should thus be reached; in any case, increasing power beyond   should be avoided, as it may cause strongly overestimated values of   ).For each distribution, several sets of  = 40-45 values of energy are calculated, by generating several times the random amplitudes of the 1000 rectangles that form the bunch and/or by varying the values of   .Once the different sets of   have been calculated, it is assumed that the distribution of the bunch is unknown.The bunch is then modeled as a series of  = 40 rectangular slices of unknown amplitudes   (hence, each slice actually models a subset of / = 25 rectangles of the bunch with similar amplitudes).For each set of   values, the minimization of ( 4) is performed several times, varying the initial estimates of   (equal values, random values, etc.), for different degrees of overdetermination in (3), and so forth.Several tens of sets of values   are then generated, from which those associated with an anomalously large value of the residual (meaning bad convergence) are discarded.After this selection, the remaining data is used to estimate the statistics of the bunches.
Figure 4 presents the results obtained for 3 common distributions: Gaussian, uniform, and chi squared with 4 degrees of freedom.It appears from the figure that the method yields reasonably good estimates of the actual statistics of the bunches.

Generalization of the Technique to Any Pulse Profile
Although the technique works well to estimate the statistics of a bunch of rectangles, rectangular pulses are seldom encountered in practice, and it is clear that the above technique will not operate properly if the pulses adopt a different shape, like Gaussian, for example.It is nevertheless quite straightforward to adapt this technique to any particular pulse shape.For the moment, we will continue assuming that all pulses in the bunch have the same duration.It is also assumed that the pulses are sufficiently separated in time so that they do not overlap substantially.First, using (1), we will compute the energy transfer characteristic  1 of a single pulse adopting the profile under consideration, (), as a function of input peak power.Defining () as the profile with normalized peak power,  1 is written as It should be noted that  1 is independent of the duration of the pulse and only depends on its profile and peak power   .Indeed, if  1 () is the profile with unit duration, then the corresponding profile with duration  is given by () =  1 (/); after substituting in ( 5) and proceeding to a change of variable   = /, it appears that  1 remains unchanged when () is replaced by  1 ().Hence, using ( 1) and ( 5),  1 can be readily computed for a given profile as a function of peak power.If now we define  = ∫  1 (), the energy of a pulse with duration  and peak power   is given by   ∫ () =   .Finally, (3) is replaced by where   now represent the peak powers of the pulses.Like in the former case, the bunch is modeled by a series of  "slices," each of which represents a subset of pulses in the bunch with similar amplitudes.Like previously, each slice regroups the same number / of pulses (hence, Δ = /, where  is the duration of individual pulses).Once   have been measured, the problem can be solved using the same procedure as before, replacing (3) by ( 6); the main difference is that the NOLM characteristic in ( 6) is no longer , but  1 , which is profile-dependent.
Once again, we checked the technique numerically.Bunches of  = 1000 Gaussian pulses were generated according to a few common distributions and were used to compute the values of energies   .The energy transfer function of a Gaussian pulse was then computed according to (1) and ( 5), in function of peak power.Equations ( 6) were then solved numerically in the same way as previously (a polynomial fit of  1 was used to improve computational speed).The obtained results are presented in Figure 5.It can be seen that they are comparable to those of Figure 4.
Finally, it has to be mentioned that although the above technique is able to estimate the statistics of a bunch of pulses, the total number  of pulses in the bunch, as well as their duration, remains unknown.This complementary information could be accessed through an autocorrelation measurement, for example.

Adaptation of the Technique to the Case of Solitons
An important particular case is that of a bunch of fundamental solitons, which can be the result, for example, of the breaking of a long, ns-long pulse in a dispersive nonlinear fiber.Applying the technique developed in Section 4 directly (assuming a sech 2 profile) will yield incorrect results, however, if we disregard the well-known property of first-order solitons that their duration is inversely proportional to the squared root of their peak power.In a fiber with anomalous dispersion  2 (in ps 2 /km) < 0 and nonlinear coefficient  (in Therefore, the higher the peak power, the shorter the duration of the pulse (it has to be clarified here that the pulses are solitons only in the fiber that originated them, and not, of course, in the NOLM used for the characterization, which is still supposed to be dispersion-free).In this section, we will continue regarding the bunch as a series of  slices with equal duration Δ, which regroup pulses with similar amplitudes.However, considering (7), the slices will no longer contain an equal number of pulses: the higher-amplitude slices will regroup a higher number of solitons than those corresponding to a lower value of amplitude.Hence each slice will now contain a number of pulses that is proportional to the squared root of its amplitude: where   are the durations of the solitons included in the th slice, calculated using (7).Hence, to estimate the amplitude distribution of a bunch of solitons on the basis of measurements of energy transmission through a NOLM, the same procedure as in Section 4 can be applied, based on (6), to determine the amplitudes   of the slices.Then, using (8), the number   of pulses in each category can be determined, which finally allows determining the statistics of the bunch.In addition, the total number of solitons in the bunch,  = ∑   , can also be estimated.Hence, thanks to the particular relation between amplitude and duration that prevails in the case of solitons, information such as the total number of pulses in the bunch and their duration is accessible to the measurement.
Finally, we will test numerically the procedure described above to estimate the statistics of a bunch of solitons.Again, we will assume that the pulses are sufficiently separated in the bunch, so that they do not overlap substantially.In the case of solitons resulting from the process of pulse breaking in a fiber, this implies that the bunch should be measured after a sufficiently long distance of propagation along the fiber, in order to ensure that the process of soliton formation and separation is completed.Here we will assume that the "duration of the bunch"   = Δ is not known.Indeed, although in practice the duration of the bunch (ns or longer) is usually accessible using affordable measurement optoelectronics, this value encompasses the pulses as well as the empty intervals separating them, which are not accounted for in   (strictly speaking,   = ∑     is the cumulated duration of all the pulses in the bunch).With Δ unknown, the minimum number of (6) to solve the problem is increased by 1 (i.e.,  ≥  + 1).It is then convenient to use one of these equations, say the th, to eliminate Δ from the remaining equations, which then are rewritten as Once (9) are solved in   , the value of Δ is finally given by where the fact that  = 2 for sech 2 pulses was taken into account.
Again, bunches of  = 1000 sech 2 pulses are generated according to a few common distributions and are used to compute the values of energies   (which, in a real measurement, would be experimental).In addition, we also consider the distribution defined by Equations (2) and (3) in [14].This distribution corresponds to the approximate analytical solution of the pulse breaking problem obtained in the frame of the work of Zakharov and Shabat [13], when the long original pulse is nearly rectangular.For each distribution, the average peak powers of the pulses in the bunch range as a few tens of W. For the NOLM parameters, we use   = 125 W and Δ = −/7 in (1).Such value of critical power is obtained, for example, using the scheme of Figure 1 with 15 m of high nonlinearity fiber ( = 10 rad −1 km −1 ) in the loop.With these values, only adjustable attenuation is required at the NOLM input (amplification is not needed; i.e.,   < 1).The energy transfer function  1 is then computed assuming a sech 2 profile and fitted using a polynomial.Equations ( 9) are then solved numerically using the least squares method, yielding   .We consider  = 40 sections.After calculating Δ through (10), the number of pulses in each slice is determined using (8) (the dispersion and nonlinear coefficients are those of a standard telecommunications fiber at 1550 nm; we used for the latter 2/3, corresponding to circularly polarized solitons).For each distribution, the procedure is repeated several times, varying the values   , the initial estimates of   , the degree of overdetermination of the equations, and so forth, in order to increase the volume of data.Sets of solutions associated with significantly higher-than-average values of residuals are discarded.Finally, the obtained results are presented in Figure 6.Overall a good agreement is observed between the estimated and targeted statistical distributions.Besides, in each case, the estimated total number of solitons in the bunch is found to be within ∼5% of its actual value,  = 1000.
Let us finally provide a few additional details in the case of Figure 6(d).The bunch is the one that would be formed, according to [14], from the breaking of a nearly rectangular pulse with 30 W peak power and a duration of 106.8 ns, in a standard telecommunications fiber ( 2 = −21.7 ps 2 /km at 1550 nm and 2/3 = 1 rad −1 km −1 ).With these realistic values of the parameters,  = 1000 solitons are formed, with amplitudes up to 120 W, corresponding to 4 times the amplitude of the original pulse, and durations down to 0.75 ps (full width at half maximum).Although the values of amplitudes given by Equations (2) and (3) in [14] are deterministic and do not define a probability density function, they yield the distribution shown by the transparent histograms in Figure 6(d).The reconstructed data (solid histograms) present a good overall agreement with the targeted distribution, although a noticeable discrepancy appears at the upper end of the distribution: the large number of pulses in the last category of amplitudes (117-120 W) is strongly underestimated.Instead, the technique finds a certain number of solitons with amplitudes beyond 120 W (up to ∼150 W).It has to be noticed however that if we sum the heights of all the solid rectangles beyond 117 W and calculate the corresponding average peak power of these solitons (123.3W), we get values that reproduce reasonably well both the height and abscissa of the last rectangle in the targeted distribution (black dot in Figure 6(d)).

Conclusions
In this work we propose and develop a technique to characterize the statistics of a bunch containing a large number of optical pulses with random amplitudes.The technique exploits the nonlinear switching characteristic of a fiber NOLM, whose ultrafast response time is provided by the optical Kerr effect.For proper adjustments of the NOLM parameters, the energy of the bunch transmitted by the device as a function of input power is specific of the statistics that obey the amplitudes of the pulses in the bunch.Therefore, on the basis of the knowledge of the energies of the bunch transmitted by the NOLM (which are readily measured in practice using a low-frequency measurement setup), we are able to estimate through numerical calculations the statistics of the pulses in the bunch.The method involves estimating the solution of a system of nonlinear algebraic equations using the least squares method.The technique is demonstrated numerically in the case of rectangular pulses, Gaussian pulses, and in the case of a bunch of solitons.The technique is attractive in particular in the frame of supercontinuum generation using ns pulses, in which pulse breaking plays an essential role, and for the characterization and understanding of the noiselike pulsing regime in mode-locked fiber lasers.

Figure 1 :
Figure 1: Polarization-imbalanced NOLM with output polarization selection.The use of a QWR at the NOLM output in combination with a polarizer oriented at 45 ∘ in the proper direction allows selecting the circular polarization component orthogonal to the input polarization.

Figure 2 :
Figure 2: (a, d) Series of rectangles with random amplitudes obeying (a) a Gaussian distribution with mean 1 and standard deviation 0.4 and (d) a uniform distribution on the interval [0, 1]; (b, e) when arranged in decreasing order of amplitude, each series tends to a specific intensity profile (materialized by the dashed curves);(c, f) illustration of the basic principle of the method: for the same energy (shaded area) of the input bunch, the bunch energy at the NOLM output is very different for each distribution (in this example, it is higher in the uniform case).

Figure 4 :
Figure 4: Histograms presenting, in 40 ranges of amplitude, the number of slices determined by the technique for three distributions: (a) Gaussian, (b) uniform, and (c) chi squared with 4 degrees of freedom.The actual distributions that were used to generate the bunches and compute the   values are superposed for comparison.

Figure 5 :
Figure 5: Histograms presenting, in 40 ranges of amplitude, the number of "slices" (subsets of Gaussian-shaped pulses) determined by the technique for three distributions: (a) Gaussian, (b) uniform, and (c) chi squared with 4 degrees of freedom.The actual distributions that were used to generate the bunches and compute the   values are superposed for comparison.The phase bias Δ = −0.2and   = 1.

Figure 6 :
Figure 6: Histograms presenting, in 40 ranges of amplitude, the number of pulses determined by the technique, for four distributions: (a) Gaussian, (b) uniform, (c) chi squared with 4 degrees of freedom, and (d) the statistics derived from the work of Zakharov and Shabat (see main text for details).The actual distributions that were used to generate the bunches and compute the   values are superposed for comparison.The black dot in (d) shows the combined contribution of the rightmost sections of the retrieved histogram (bracket).