Xiaole Fang
email: xiaole.fang@univ-fcomte.fr

Benjamin Wetzel

Jean-Marc Merolla

John M Dudley

Laurent Larger

Christophe Guyeux

Jacques M Bahi

Noise and chaos contributions in fast random bit sequence generated from broadband optoelectronic entropy sources

Keywords: Random number generation, chaos, optoelectronics, noise, entropy sources, statistical tests

During the last 4 years, chaotic waveforms for random number generation found a deep interest within the community of analogue broadband chaotic optical systems. Earlier investigations on chaos-based RNG were proposed in the 90s and early 2000, however mainly based on piecewise linear (PL) 1D map, with bit rate determined by analog electronic processing capabilities to provide the PL nonlinear function of concern. Optical chaos came with promises for much higher bit rate, and entropy sources based on high complexity (high dimensional) continuous time (differential) dynamics. More specifically in 2009, Reidler et al. published a paper entitled "An optical ultrafast random bit generator", in which they presented a physical system for a random number generator based on a chaotic semiconductor laser. This generator is claimed to reach potentially the extremely high rate of 300 Gb/s. We report on analysis and experiments of their method, which leads to the discussion about the actual origin of the obtained randomness. Through standard signal theory arguments, we show that the actual binary randomness quality obtained from this method, can be interpreted as a complex mixing operated on the initial analogue entropy source. Our analysis suggests an explaination about the already reported issue that this method does not necessarily require any specific deterministic property (i.e. chaos) from the physical signal used as the physical source of entropy. The bit stream randomness quality is found to result from "aliasing" phenomena performed by the post-processing method, both for the sampling and the quantization operations. As an illustration, such random bit sequences extracted from different entropy sources are investigated. Optoelectronic noise is used as a non deterministic entropy source. Electro-optic phase chaotic signal, as well as simulations of its deterministic model, are used as deterministic entropy sources. In all cases, the extracted bit sequence reveals excellent randomness.

I. INTRODUCTION

Random Number Generators (RNGs) are widely used in science and technology. They are a critical component in modern cryptographic systems, communication systems, statistical simulation systems, and any scientific area incorporating Monte Carlo methods and many others [START_REF] Nielsen | Quantum Computation and Quantum Information[END_REF]- [START_REF] Stinson | Cryptography: Theory and Practice[END_REF]. The unpredictability of the bit stream and the speed at which the random bits can be produced are usually reported as very important aspects in the quality of the generated bit sequence. Other factors like system complexity, cost, reliability and so on, are also important for establishing successful RNGs. There are usually two methods for RNGs: One relates to deterministic algorithms implemented in hardware and software, the pseudorandom numbers being generated from a single "seed" (such generators are named pseudorandom number generators, or PRNGs in short); Another one counts on high entropy signals, whether from mainly nondeterministic and stochastic physical phenomena (see [START_REF] Williams | Fast physical random number generator using amplified spontaneous emission[END_REF]- [START_REF] Yamazaki | Performance of random number generators using noise-based superluminescent diode and chaos-based semiconductor lasers[END_REF]), or from deterministic but chaotic dynamical systems [START_REF] Yamazaki | Performance of random number generators using noise-based superluminescent diode and chaos-based semiconductor lasers[END_REF]- [START_REF] Hirano | Characteristics of fast physical random bit generation using chaotic semiconductor lasers[END_REF]. A potential advantage of the latter physical high entropy signal, arises in its deterministic features that might be used to achieve chaos synchronization as it has been already demonstrated [START_REF] Pecora | Synchronization in chaotic systems[END_REF] and widely used for secure optical chaos communications [START_REF] Argyris | Chaos-based communications at high bit rates using commercial fiberoptic links[END_REF]. However, synchronization possibility of the random binary sequence extracted from the chaotic physical signal is still an open problem, which resolution could lead to the efficient and practical use of the one-time pad (a symmetric encryption algorithm derived from the Vernam cipher, which is proven to be impossible to crack if used under appropriate conditions).

The PRNGs based on deterministic algorithms can be implemented in any computational platform, some can even be easily adapted with discrete chaotic iterations to improve output quality of randomness [START_REF] Bahi | An optimization technique on pseudorandom generators based on chaotic iterations[END_REF]. They however suffer from the vulnerability that the future sequence can be deterministically computed if the seed or internal state of the algorithm is discovered. The main advantage of PRNGs is that no hardware cost is added, and the speed is only counted on digital processing hardware. Their algorithms are developed to prevent guessing of the initial conditions, and the rate might be slowed down while increasing the complexity of such algorithms.

Physical RNGs rely on chaotic or stochastic physical processes. Such random number generators are building the random bits from inherently random or chaotic physical process [START_REF]Randomness requirements for security[END_REF], for example, radioactive decay [START_REF] Walker | Hotbits: Genuine random numbers, generated by radioactive decay[END_REF], chaotic electrical and optical circuits [START_REF] Bernstein | Sercure random number generation using chaotic circuits[END_REF], and so on. Up to recently, the implementations of physical random generators have been limited to much slower rates than PRNGs because of limitation of the mechanisms for extracting bits from physical randomness without degrading statistical properties. Typically 10 Mb/s could be achieved by using electronic oscillator jitter [START_REF] Bucci | A high-speed oscillator-based truly random number source for cryptographic applications on a smart card ic[END_REF] and 4 Mb/s using quantum optical noise [START_REF] Dynes | A high speed, postprocessing free, quantum random number generator[END_REF]. One should notice however, that such physical implementation was recently directly developed at the chip level in personal computer processors, finally achieving more than reasonable speed performances (> 3 Gb/s), and actually also very good randomness quality [START_REF] Taylor | Digital randomness[END_REF].

Considerable improvements for the rate of chaotic random bits generation have been however reached by using a semiconductor laser in the presence of external feedback [START_REF] Mukai | New route to optical chaos: Successivesubharmonic-oscillation cascade in a semiconductor laser coupled to an external cavity[END_REF], a well known setup in chaotic optical systems. The dynamical processes involved in optical systems can indeed be very fast. Moreover, high complexity chaotic dynamics can be practically obtained, whether due to intrinsic complex nonlinear coupling between light and matter interactions in lasers, or due to the presence of a large delay feedback cavity enabling dynamics with large number of degrees of freedom. Chaotic optical signal might consist of pulses with a width of few 10ps and with random amplitude and time positions, which provide attractive potentials to easily generate random bits at fast rates. In [START_REF] Uchida | Fast physical random bit generation with chaotic semiconductor lasers[END_REF], a first attempt already reached 1.7 Gb/s RNG, the physical randomness originating from two independent chaotic semiconduct lasers. Each laser intensity signal is practically sampled at an incommensurate rate with respect to the individual optical feedback delay times. Then a threshold value is set for comparison with each signal level and to obtain a Boolean sequence. Lastly the random bit sequence is produced by executing a XOR function between the two Boolean sequences. More recently, Reidler and colleagues [START_REF] Reidler | Ultrahigh-speed random nuber generation based on a chaotic semiconductor laser[END_REF], [START_REF] Kanter | An optical ultrafast random bit generator[END_REF] claim that they successfully demonstrated another method in generating random bit sequence from ultra fast optical chaos, at much faster rate. In such method, the output of a single chaotic laser, with the optical feedback delay time incommensurate with the sampling clock frequency, was digitized by an 8-bit analog-to-digital converter (ADC, practically provided by an ultra-fast digital scope). Then the difference between adjacent, but not nearest, points from the 8bit digitized time series is performed (it is defined as a pseudo-"derivative" operation). At last, a few LSBs only of the values resulting from the subtracted samples are retained to generate the binary sequence. Following this combination of broadband photonic chaos and digital post-processing, generation rate as high as 300 Gb/s are claimed. Many additional papers have then appeared [START_REF] Argyris | Implementation of 140 gb/s true random bit generator based on a chaotic photonic integrated circuit[END_REF]- [START_REF] Li | Direct generation of all-optical random numbers from optical pulse amplitude chaos[END_REF], utilizing similar bit stream extraction method, but on different alternative photonic setups. These reported works have claimed to have achieved comparable high speed and high randomness quality bit stream. None of these papers has however discussed the actual and respective roles played by the photonic chaotic waveform on one side, and by the post-processing method on the other side.

In this paper, the study of using broadband optical signal to generate random binary sequence according to the method proposed by Reidler, is going to be deepen. We propose to apply the same method on the chaotic waveform generated by another class of broadband photonic oscillations, and to analyze the different post-processing steps involved in this method. We will analyze three key factors in the scheme of [START_REF] Reidler | Ultrahigh-speed random nuber generation based on a chaotic semiconductor laser[END_REF] and [START_REF] Kanter | An optical ultrafast random bit generator[END_REF]: the sampling, the difference of distant samples, and LSBs retaining.

The paper is organized in the following way. In Section II, the original architecture that we propose as the physical system to generate fast random bit sequences, is described in details. Analysis of the conversion of the physical time series into a binary random sequence as proposed in [START_REF] Kanter | An optical ultrafast random bit generator[END_REF] is recalled, but also analyzed in terms of signal processing arguments. Then an entropy evaluation for the binary sequence is processed in Section III, thus providing insight for the entropy rate capability for the generated binary sequence. Finally, the randomness quality of the obtained bit sequence is compared, depending on the used physical signal: whether a chaotic laser intensity, or a noisy signal covering a similar Fourier spectrum. The randomness quality is evaluated in Section IV via the results of four statistical testing packages. The paper ends with a discussion of the obtained results, and concludes with possible future work.

II. METHOD FOR RANDOM BIT SEQUENCE GENERATION FROM AN OPTOELECTRONIC SIGNAL

In this section we describe the physical setup from which we expect to obtain a fast random bit sequence. We also describe the binary sequence extraction method from the continuous time signal generated by the physical setup, as it was formerly proposed in [START_REF] Reidler | Ultrahigh-speed random nuber generation based on a chaotic semiconductor laser[END_REF], [START_REF] Kanter | An optical ultrafast random bit generator[END_REF]. Additionally, theoretical interpretations and discussions of this extraction method are proposed in terms of basic signal processing and sampling theory. These interpretation and discussion are intended to give insight on the possible mechanisms at the origin of the bit stream randomness quality.

A. Setup delivering a broadband optoelectronic signal

In order to additionally support our work with experiments on the generation of optical broadband signals, data recorded from physical chaos generator as well as from optoelectronic noise sources, have been studied. These experiments are different from the ones described in [START_REF] Reidler | Ultrahigh-speed random nuber generation based on a chaotic semiconductor laser[END_REF] and [START_REF] Kanter | An optical ultrafast random bit generator[END_REF], although they are also originating from broadband optoelectronic devices.

A twofold physical source of entropy has been used (see Fig. 1), both having been tested for their randomness quality. One source (referred as "Optoelectronic noise" in Fig. 1) is originating from physical noise sources in the semiconductor laser light generation process (known as RIN: relative intensity noise), in combination with the electronic noise of the photodetector and its integrated electronic amplifier (thermal noise and semiconductor photodiode junction noise, amplified by the noise figure of the electronic amplifier). A comparable (and even cheaper) optoelectronic noise source was also proposed in [START_REF] Li | Scalable parallel physical random number generator based on a superluminescent led[END_REF], [START_REF] Williams | Fast physical random number generator using amplified spontaneous emission[END_REF]. On the contrary, the other physical source of entropy is originating from a strongly deterministic process, which was used for a field experiment demonstrating (analogue) chaotic optical masking of 10 Gb/s data signals, transmitted over an installed fiber optic link [START_REF] Lavrov | Nonlocal nonlinear electro-optic phase dynamics demonstrating 10 gb/s chaos communications[END_REF]. The strong determinism of this entropy source indeed enabled to implement accurate broadband chaos synchronization at the receiver, in order to remove the chaotic masking signal and thus to retrieve the original binary data stream. The dynamics of the electrooptical phase chaos generator is ruled by a nonlinear dual delay differential equation implemented in an optoelectronic and electro-optic feedback loop.

Each of these two signals obtained from noise or chaotic optoelectronic systems, has been processed by using the method proposed in [START_REF] Reidler | Ultrahigh-speed random nuber generation based on a chaotic semiconductor laser[END_REF]. By doing so, the aim is to support our signal processing analysis on the extraction method of the bit sequence, which is inferring that in both cases, the randomness quality must be very similar. As we shall illustrate and discuss later, and as it has been pointed only in the literature [START_REF] Williams | Fast physical random number generator using amplified spontaneous emission[END_REF], chaos is not a necessary condition for a good randomness quality of the extracted bit stream. Global spectral and statistical features in the original analogue source of entropy appear to be enough, such as a broadband Fourier spectrum, and a sufficient spreading of the amplitude probability distribution. Instead of using necessarily noise corrupted signal as it is always the case in experiments, one can also generate time traces numerically, thus resulting in an even stronger determinism in the entropy source. In that case, the noise level is significantly reduced to numerical rounding and numerical integration errors. One however needs a suitable deterministic model in order to generate such numerical time traces. Our deterministic chaotic signal is modeled by the physical solution of a nonlinear dual delay dynamics ruled by the following differential equation:

θ -1 t 0 x(ξ)dξ + τ dx dt (t) + x(t) = β sin 2 [x T -x T +δT + Φ 0 ], (1)
where x T stands for the delayed signal x(t -T ), θ and τ are the characteristic times of the low and high cut-off frequency respectively, which are involved in the bandpass feedback filtering of the RF filter. From a signal processing viewpoint, such a dynamical system can be interpreted as a nonlinear delayed feedback oscillator. This oscillator is ruled by the dynamics of a linear bandpass filter h(ω), which is driven by a nonlinear transformation of two delayed "echos" (delays T and T + δT ) of the filter output, x(t). Chaotic solutions are obtained when the feedback loop gain β is high enough, of the order of 5. This gain is adjusted via the tuning of the laser light intensity. A typically observed chaotic solution is a white noise like signal which is covering the spectral range of the broadband bandpass feedback RF filter, i.e., ca . This results in a fast noise-like large amplitude signal x(t), which is expected to be suitable for high speed RNG based on physical device. It is worth noticing that this chaotic signal generation process can be viewed as a balanced equilibrium between the RF feedback filter (limiting the spectral span of the signal x(t)), and the spectral broadening performed by the nonlinear transformation (sin 2 -function of the difference delayed signals x T -x T +δT ). The offset phase Φ 0 is typically adjusted through the static interference condition, which interference phenomena is physically at the origin of the sin 2 nonlinear transformation. A more accurate description of the generated signal x(t) should also include (small amplitude) noise sources in the equation. The latter noise is actually of the same kind than the one involved in the noisy optoelectronic signal (laser and photodiode noise, also complemented here by the RF electronic amplifier noise). Equation ( 1) can be confidently used solely without the noise contribution, when the large amplitude chaotic motion only is of interest. To support this asumption, one might notice that such numerically obtained waveforms were found to exhibit very good qualitative agreement with the one observed in the experiment [START_REF] Lavrov | Electro-optic delay oscillator with non-local non linearity: optical phase dynamics, chaos, and synchronization[END_REF]. The proposed model has been also successfully used to derive analytically [START_REF] Weicker | Crenelated fast oscillatory outputs of a two-delay electro-optic oscillator[END_REF] most of the bifurcation features indeed observed in the experiment, thus indicating that the large amplitude solutions can be confidently calculated numerically from the noise-free model.

B. Extraction method for the random binary sequence

A schematic view of the algorithm used to extract a random binary sequence from a broadband physical signal, as proposed in [START_REF] Kanter | An optical ultrafast random bit generator[END_REF], is depicted in Fig. 2. On the basis of a physical setup delivering a broadband signal, as the one described in the previous section, a real time oscilloscope is first involved to perform an analogue to digital conversion of the output signal of the setup. This conversion is typically achieved via an 8-bit digitizer at a sampling rate of 40 GHz. In the next subsection, we will discuss from the signal theory viewpoint some particular processing issues that are found to significantly contribute to the actual randomness of the final binary sequence. More precisely, sampling issues will be discussed, quantization issues, and also post-processing operations (such as distant sample difference, and LSB-only retaining). This signal processing is performed before getting the actual final random binary sequence, to be tested for their randomness quality via standard statistical test suites.

1) Sampling issues: aliasing for enhanced entropy

In the following, we assume that samples are originally acquired by a real time digital scope measuring a broadband complex time trace. Such equipment is designed to follow the classical Shannon sampling theorem: the sampling rate f S is matching the instrument analogue input bandwidth, which defines the maximum Fourier frequency f M that can be captured by the instrument. The Shannon sampling theorem indeed states that a limited bandwidth signal can be digitized without loss of information, when the sampling frequency is at least twice the maximum signal frequency. The sequence of the samples {s n = x(nT ), n ∈ Z} can be defined as a function of the continuous time as follows:

s(t) = x(t) • ⊔⊔ TS (t), (2)
where A typical illustration of the proof for the sampling theorem is presented in Fig. 3, as one describes the spectrum of such a sampled signal s(t),

⊔⊔ T (t) = k=+∞ k=-∞ δ(t -k T ).
S(ν) = FT[s(t)] = X(ν)⋆FT[⊔⊔ T (t)] = 1 T X(ν)⋆⊔⊔ 1/T (ν), (3)
where we have used the well known result that the Fourier Transform (FT) of a comb is also a comb. The convolution product in the Fourier domain reveals that the spectrum of the sampled signal is the result of the superposition of an infinite number of regularly spaced replica of the original signal spectrum X(ν) =FT[x(t)], two neighboring replica being separated by the sampling frequency f S = 1/T . Thus, if the maximum frequency f M of the bounded support of X(ν) is less than f S /2, the replicated spectra do not overlap (see Fig. 3(a)). It is then obvious that a suitable window filtering of the sampled signal s(t) allows to recover in the Fourier domain exactly the same spectrum as the one of the original signal x(t). Such a filter typically transmits perfectly all the Fourier components in a frequency band such as [-f S /2, +f S /2], and rejects all the other Fourier components for the other frequency ranges. When undersampling is used, the aliasing phenomenon occurs in the Fourier domain. It consists then of overlaps between the replicated spectra due to the comb convolution. The actual spectrum of the sampled signal s(t), can be viewed as a complex mixing of the frequency components of the original signal x(t), due to the overlapped replica of X(ν) =FT[x(t)]. The procedure of selecting only one sample every n from the original sampled sequence, is thus equivalent to an aliasing operation with an undersampling of order n. The original goal of cementification of the extracted sample sequence, can thus be viewed as an aliasing technique resulting in a complex mixing of the original frequency components. The consequence is an increased entropy of the output sequence. When viewed in the time domain, this procedure results in the vanishing of the short time correlations, since a long time intervalle is then separating two successive samples, compared to the width of the autocorrelation function of the original signal. On the contrary, the short time scales correlations are necessarily present when the conditions of the sampling theorem are fulfilled (two successive samples would thus keep the information of short time scale correlations). Another consequence is that such an operation is unidirectional, in the sense that original information is actually lost after an aliasing process. Because of the complex mixing of the Fourier frequency components, the original spectrum cannot be recovered with a "simple" unmixing procedure.

2) Further post-processing: difference sequence between distant samples

In [START_REF] Reidler | Ultrahigh-speed random nuber generation based on a chaotic semiconductor laser[END_REF] and [START_REF] Kanter | An optical ultrafast random bit generator[END_REF], computing the difference sequence between two neighbor samples are named as "derivative". However, mathematically speaking, the term "derivative" of x(t) is used for the asymptotic value (x(t+ ∆t)-x(t))/∆t when ∆t → 0. In the physical case of a finite sampling rate, the neighbor samples are obviously not infinitely close in time, hence we prefer not to use "derivative" here. More precisely, we are dealt here with significantly separated samples in time, since strong aliasing is first operated (see Section II-B1), with an undersampling number up to n = 16. The initial 40 GHz sampling rate is respecting the oscilloscope analogue input bandwidth of 12 GHz, but the final series obtained after retaining 1 sample over 16, is corresponding to a 2.5 GHz undersampling rate. The samples obtained after this distant sample difference (which will be called from here DSD) operation can be described as follows:

{d n k } k∈N = {x[k T ] -x[(k -n)T ], k ∈ N}. (4)
If we try again to analyze in the Fourier domain the meaning of this second processing, one obtain the following expression for the Fourier spectrum of n-undersampled difference signal:

D(ν) = 2 i e -iπνnT X(ν) sin(πνnT ) ⋆ ⊔⊔ 1 nT (ν).
(5) This expression reveals a so-called channeled spectrum filter, which applies a periodic sinusoidal modulation of the original spectrum X(ν). One could notice that the maximum transmission of this filter is centered at half the undersampling rate ((nT ) -1 /2) where aliasing is maximally symmetric (thus somehow selecting the frequency components that are most affected by aliasing), and the frequencies of null transmission are centered at zero and ±(nT ) -1 (where the aliasing phenomenon is the less pronounced in the Fourier spectrum, as long as n is not too large). When focusing on the low frequency domain only, another comment about the action of this DSD processing could be made: the very low frequencies are filtered out, which consequence is to asymptotically set to zero the mean value of the corresponding sample set, and thus also improving the symmetry around zero of the amplitude probability distribution. The Fourier analysis of the DSD processing is however not as obvious as for the aliasing issue in terms of randomness enhancement, or entropy amplification. A more meaningful discussion can however be made through the analysis of the statistical sample distribution of the DSD compared to the original one. More precisely, Fig. 4a shows the evolution of the amplitude statistics when the DSD processing is iterated several times. One clearly sees that the statistics is more and more symmetric, resembling closer and closer to a Gaussian distribution. We have checked that introducing noise in the simulation does not change this result. This convergence towards a Gaussian statistics through DSD can be qualitatively explained through the analysis of the DSD principle. Since the difference is performed between the same sample sequence but shifted in time over a quantity large enough compared to the correlation time, one can interprete the DSD as the superposition of nearly independent pseudorandom processes. The central limit theorem can then be used to explain qualitatively the amplitude distribution convergence towards a Gaussian one (limit of the amplitude distribution for the superposition of an asymptotically large number of independent random processes). At the same time, as more and more DSD are operated, the amplitude range is increased along the horizontal axis, whereas the maximum of the statistics along the vertical axis is inversely decreased. Figure 4c shows the numerical evidence of a hyperbolic relation between width and height for the successive iterated processings. This asymptotically Gaussian distribution obtained after a few DSD steps, finally prepares a kind of optimal statistical conditions for the last post-processing operation proposed by [START_REF] Reidler | Ultrahigh-speed random nuber generation based on a chaotic semiconductor laser[END_REF], leading to the final random bit sequence: LSBs retaining.

3) LSB retaining, and quantization noise effect

Surprisingly, the qualitative signal processing analysis of this last post-processing step, is involving very similar theoretical insight compared to the ones related to the aliasing phenomenon underlined in section II-B1. The main idea for the analysis of this last step of LSB retaining, actually originates from an elegant and powerful analogy between the temporal discretization occurring during sampling, and the amplitude discretization occurring as quantization is concerned [START_REF] Widrow | Statistical theory of quantization[END_REF]. Following the results of this statistical theory of quantization, one finds that the actual consequence of the LSBs retaining method is practically to provide a nearly constant (flat) amplitude probability distribution for the quantization noise, for quantization levels as high as the root mean square of the signal to be quantized. In terms of LSBs retaining, this means that flat noise distribution is obtained for the corresponding LSBs amplitude, even if the LSB maximum amplitude reaches the root mean square of the quantized signal. As illustrated in Fig. 9 of Ref. [START_REF] Widrow | Statistical theory of quantization[END_REF], the flat amplitude distribution results in a kind of aliasing from the original symmetric (Gaussian) amplitude distribution, which is segmented and superimposed over the central small amplitude interval limited by the resolution of the quantized amplitude (i.e. the amplitude range encoded by the retained LSBs).

Thus, after a complex mixing in the Fourier frequency domain due to a sampling theory aliasing effect, the LSB retaining process results in a strong flattening of the amplitude probability distribution for the retained LSBs. This flattening of the amplitude probability distribution can be also related to a kind of aliasing, but performed on the initial Gaussian amplitude probability distribution, according to the statistical theory of quantization. The Gaussian profile is implicitely provided through Reidler's method when a sufficiently high number of DSD processing is performed.

4) Discussion of the randomness origin

A straightforward issue can then be raised about the actual source of randomness leading to the final bit sequence, as proposed in [START_REF] Reidler | Ultrahigh-speed random nuber generation based on a chaotic semiconductor laser[END_REF]. The randomness origin has been many times attributed to the chaotic character of the solution generated by the original physical system: A SC laser diode subject to proper optical feedback, which is well known to exhibit chaotic motion. However, from the previously analyzed postprocessing steps, no single argument related to any chaotic property of the original signal was involved. The only necessary requirement was to have a certain broadband character in the original signal, such that first standard aliasing could occur in the Fourier domain. Second the DSD processing leads to a symmetrization of the amplitude statistics, with a convergence towards a Gaussian distribution. And last, LSBs retaining performs naturally a flattening of the final amplitude probability distribution corresponding to the selected bits.

To investigate this issue, we performed a similar analysis as the one done in the previous subsections on the chaotic motion of a nonlinear electro-optic delay dynamics, but with a physical signal a priori originating from physical noise sources only, without any deterministic chaotic compound. This signal is chosen to be the output of the amplified photodiode of the same setup, but without the nonlinear delayed feedback loop at the origin of the chaotic time series: The amplified photodiode signal is issued from the laser intensity noise, it comprises also the photodiode junction noise and the electronic amplifier noise (see "Optoelectronic Noise" output in Fig. 1). Although the electrical signal level is significantly lower, we used the scope magnification to get a time trace of a comparable amplitude with respect to the scope vertical amplitude range, thus resulting in an effective digital scope quantization over a comparable number of bits with respect to the chaotic signal.

We have reported in Fig. 4 and 5 the statistics evolution of the digitally acquired optoelectronic noise signal and its width / height evolution. The figures clearly show very similar features. From this rough analysis of the influence of the two physical signals (the optoelectronic noise and the electrooptic chaos), we realize that the post-processing leads to qualitatively equivalent final bit sequence. This observation, and the previous analysis of the post-processing steps, support the assumption, at least qualitatively, that the randomness of the final bit sequence might be mainly issued from the postprocessing steps. The chaotic feature of a time series appears as an actually sufficient but not required condition, for the generation of a random bit sequence when the discussed postprocessing is used (undersampling, DSD, and LSBs retaining). A simple noisy signal with similar spectral extend is found to lead to similar final output bit stream, as already reported previously [START_REF] Williams | Fast physical random number generator using amplified spontaneous emission[END_REF], but not yet analyzed and interpreted as we have proposed.

III. EFFECT OF NOISE ON THE ENTROPY RATE IN THE BINARY SEQUENCE

In this section, the time evolution of the entropy in the final binary sequence is evaluated under different choices for the method used to build the final random bit stream from the chaotic signal. The aim is to get insight in the origin of the entropy creation mechanism involved in the construction of the final random bit stream. More precisely, we aim at discriminating under which conditions the deterministic feature of the chaotic signal (the determinism coming from the dynamics described by Eq.( 1)) is indeed involved in the entropy of the extracted bit stream. To achieve this goal we reproduce the method proposed in [START_REF] Mikami | Estimation of entropy rate in a fast physical random-bit generator using a chaotic semiconductor laser with intrinsic noise[END_REF], which is intended to measure the sensitivity to initial condition (SIC) of the deterministic chaotic motion in the presence of additional small noise. This measure consists in calculating the temporal entropy evolution for the generated binary random sequence, with respect to several different noise realizations.

A. Introducing noise in the simulated chaotic dynamics

For the entropy calculation, we first consider a transientfree chaotic solution of Eq.(1) (β = 5). To achieve such a solution, Eq.( 1) is integrated under the proper parameter conditions known to lead to a high complexity chaotic solution. This preliminary numerical integration is performed over a duration long enough compared to the slowest characteristic time scale of the dynamics (θ), so that the asymptotic trajectory is free of any transient. Once this corresponding chaotic attractor is supposed to be reached via the numerical integration, this asymptotic solution can be associated to a single temporal waveform covering only the longest time delay of the dynamics, i.e. T + δT : this is defining the initial condition of the corresponding delay based, and noise-free, chaotic dynamics, from which noise influence will be explored. We then introduce in the right hand side of Eq.( 1) an arbitrary small additive noise term (small perturbation along the chaotic trajectory). The noise amplitude is arbitrarily set so that the Signal-to-Noise Ratio (SNR) is 40 dB. After further integrating Eq.( 1) with the noise term and starting from the initial condition corresponding to the calculated noise free chaotic trajectory, one is able to obtain a continuously noise-perturbed chaotic trajectory. When repeating this calculation with several different noise realizations, one then expects to observe the effect of SIC when comparing the different noise perturbed chaotic trajectories. This property manifests itself through a progressive amplification (as time is running) of the small perturbations materialized by the added noise. Comparing the different calculated waveforms, they consequently looks all the same right after the noise addition, but they split apart (differently for each pair of such time series) after a typical time scale related to the inverse of the largest Lyapunov exponent of the chaotic dynamics (see Ref. [START_REF] Mikami | Estimation of entropy rate in a fast physical random-bit generator using a chaotic semiconductor laser with intrinsic noise[END_REF] for details). Two such simulated waveforms are represented in Fig. 6, after the undersampling procedure, and before the DSD and bit retaining processes for the final extracted binary sequence. The waveforms thus do not appear anymore as continuous in time due to undersampling. For these two realizations and with the chosen SNR for the noise amplitude, one clearly sees that the two time series separate one from each other after a typical time scale of ca. 300 ns. This time scale is of the order a few tens of the largest time delay T + δT , which is corresponding to a few tens round trips of the chaotic signal in the nonlinear delayed feedback loop. This is fully consistent with the typical order of magnitude of the inverse largest Lyapunov exponent, i.e. it is of the order of the largest time delay in the dynamics.

B. Entropy estimation for each bit cell

Many different noisy chaotic time series (N = 10 3 ) are simulated to generate as many random bit sequences. Each realization is calculated from the same initial conditions (the noise free chaotic waveform over one largest delay time interval), but with different added small noise perturbations. From each obtained time series, one can explore various bit stream extraction methods, e.g. with or without DSD, or even with several successive DSD processing, with the LSB retaining or with the MSB, . . . For a fixed bit extraction method, the N obtained bit sequences can be used to calculate, at each time t k of a new extracted bit, the probabilities P 0 (t k ) and Fig. 6: Two simulated waveforms of chaotic laser intensity starting from the same initial conditions with different noise sequences added at time 0. The noise amplitude is set so that SNR is 40 dB (the signal energy being calculated on the noisefree chaotic trajectory) P 1 (t k ) for obtaining a bit 0 or 1 respectively. This time varying probability distribution is then used to calculate how the statistical bit entropy evolves in time,

H(t) = - 1 i=0 P i (t) • log 2 P i (t). (6)
As described in [START_REF] Mikami | Estimation of entropy rate in a fast physical random-bit generator using a chaotic semiconductor laser with intrinsic noise[END_REF] if SIC of a chaotic dynamics is indeed involved in the final bit sequence, the N extracted bit sequences are initially strongly correlated. This is because the bits are originating from the same initial chaotic waveform. Consequently, the influence of the small added noise is negligible at the initial times of the deterministic chaotic dynamics: the entropy at small t is expected to be close to zero if indeed the dominating phenomena is deterministic. However, as time is evolving, SIC is amplifying the influence of the small noise amplitude on the large amplitude chaotic motion, and the N bit streams realizations are more and more decorrelated, leading progressively to a maximum binary entropy of 1. This unity entropy means an equal probability for obtaining bit zero and one, independently of any deterministic motion. The influence of the small added noise term is then dominating the output random bit stream. One should notice here that the bit extraction method used in Ref. [START_REF] Mikami | Estimation of entropy rate in a fast physical random-bit generator using a chaotic semiconductor laser with intrinsic noise[END_REF] concerns MSB only. In our case, we apply the complex bit extraction method with multiple steps (aliasing, DSD, LSBs retaining) proposed in [START_REF] Reidler | Ultrahigh-speed random nuber generation based on a chaotic semiconductor laser[END_REF], and we are interested in the resulting entropy growth rate under these particular conditions, with actually expected significant differences.

Figure 7 represents the obtained binary entropy calculated with different bit extraction methods. In the first column, the binary entropies obtained from a direct 8-bits ADC for different bits from LSB to MSB, are plotted as a function of time. According to the description of [START_REF] Mikami | Estimation of entropy rate in a fast physical random-bit generator using a chaotic semiconductor laser with intrinsic noise[END_REF], memory time is defined as the time required for the entropy to reach a value close to one. One clearly sees that as the bits chosen for the random bit stream moves from LSB to MSB, the memory time of the related bit cell is increasing. This illustrates that MSBs are most reflecting the deterministic features of the signal, whereas this determinism is earlier lost as smaller LSBs are concerned: The SIC amplification of the small initial differences (noise realizations in that case) occurs naturally earlier for LSBs than MSBs. Differently speaking, this illustrates the fact that deterministic properties are more prononced with MSBs, or equivalently, LSBs are containing less deterministic features than MSBs. Columns 2 and 3 in Fig. 7 are showing similarly this determinism loss from MSBs to LSBs, when one and two DSD steps are processed respectively. One clearly sees that the smooth entropy transition, which is a signature of the chaos determinism, even completely disappears for LSBs (a zoom-in over the first ns, would show that the entropy already starts at values very close to unity).

The conclusion on Fig. 7 is that fastest entropy rate (down to the actual sampling) can be achieved when LSBs are used and when several DSD processing are performed. This is however achieved at the cost a full lost of determinism (zero memory time, without any smooth entropy growth). This corresponds to the plots on the upper right positions, for which unit entropy is already achieved very close to the time origin. This support the fact that deterministic chaos does not exist anymore in the obtained final bit stream randomness. On the contrary, the MSBs are showing a non-zero memory time, and thus a signature of a maintained determinism. The most pronounced determinism is revealed in the lower left plots, for MSB and without any DSD. One could notice that this actually corresponds to the 1-bit ADC used in [START_REF] Uchida | Fast physical random bit generation with chaotic semiconductor lasers[END_REF], where deterministic chaotic origin does contribute to the randomness of the final bit stream. In this case, the good randmoness quality is obtained only by carefully combining two uncorrelated chaotic signals (two chaotic lasers, with incommensurate characteristic time scales).

Figure 8 shows the plots of every bit cell entropy averaged over 10 3 trajectories. Each plot of entropy is obtained as a function of time for an ensemble of time series starting with exactly the same initial condition at time t = 0. Eight plots are shown in Fig. 8 corresponding to eight different position bit cell of the value. These curves are the smoothed versions, due to averaging, of the curves represented in the first column of Fig. 7. Again, it can be seen that more time is required to converge to a unity entropy when using MSB compared to the use of LSB. Differently speaking, the memory time depends on bit cell selecting, MSB and LSB appearing as the slowest and fastest entropy increasing rate, respectively. IV. STATISTICAL TESTS Additionally to the previous signal theory analysis of the processing steps used in the bit extraction method, this section is intended to qualify the final bit stream in terms of their benchmarking from several standard randomness test suites. We thus verify in this section that the analyzed and used method proposed in [START_REF] Reidler | Ultrahigh-speed random nuber generation based on a chaotic semiconductor laser[END_REF] and [START_REF] Kanter | An optical ultrafast random bit generator[END_REF] have led also in our experiment to quasi-equivalent random bit stream quality, whether from the deterministic EO phase chaos generator or with the optoelectronic noise source.

A. The tested streams

First of all, we give here a brief description of the tested methods that have been formerly proposed in [START_REF] Reidler | Ultrahigh-speed random nuber generation based on a chaotic semiconductor laser[END_REF] and [START_REF] Kanter | An optical ultrafast random bit generator[END_REF]. On the one hand, in [START_REF] Reidler | Ultrahigh-speed random nuber generation based on a chaotic semiconductor laser[END_REF], authors have used a DSD method to generate a random bits stream. The chaotic laser signal is sampled by using a 2.5 GHz ADC, and then 5 LSBs of every first DSD value are joined together to generate the final random sequence. On the other hand, in [START_REF] Kanter | An optical ultrafast random bit generator[END_REF], the chaotic laser signal is sampled thanks to a 20 GHz ADC. Then the DSD operation is processed 4 times and 8 LSBs of each value are joined to produce the pseudorandom bit stream.

These two schemes are both adapted to optoelectronic noisy signal. The generated streams sourced from chaotic laser and noise are compared by standard statistical tests in the next subsections.

B. Statistical tests

Considering the properties of binary random sequences, various statistical tests can be designed to evaluate the assertion that the sequence is generated by a perfectly random source. We have performed some statistical tests for the optoelectronic noise and electro-optic chaos generators proposed here. These tests include NIST suite [START_REF] Rukhin | A statistical test suite for the validation of random number generators and pseudo random number generators for cryptographic applications[END_REF], DieHARD battery of tests [START_REF]Diehard test suite[END_REF], ENT program [START_REF] Walker | Ent: a pseudorandom number sequence test program[END_REF], and and some comparative test parameters. A brief description of each of the aforementioned tests is given in the following paragraphs.

1) NIST statistical test suite Among the numerous standard tests for pseudorandomness, a convincing way to show the randomness of the produced sequences is to confront them to the NIST (National Institute of Standards and Technology) Statistical Test, because it is an up-to-date test suite proposed by the Information Technology Laboratory (ITL). A new version of the Statistical Test Suite has been released in August 11, 2010.

The NIST test suite SP 800-22 is a statistical package consisting of 15 tests. They were developed to test the randomness of binary sequences produced by hardware or software based cryptographic pseudorandom number generators. These tests focus on a variety of different types of non-randomness that could exist in a sequence.

For each statistical test, a set of P -values (corresponding to the set of sequences, P T ) is produced. The interpretation of empirical results can be conducted in various ways. In this paper, the examination of the distribution of P T to check for uniformity is used. The distribution of P-values is examined to ensure uniformity. If P T 0.0001, then the sequences can be considered to be uniformly distributed.

In our experiments, 100 sequences (s = 100), each with 1,000,000-bit long, are generated and tested. If the P T of any test is smaller than 0.0001, the sequences are considered to be not good enough and the generating algorithm is not suitable for usage. In Table I, the random streams generated by the chaotic laser and by the noisy signal have both obtained a 100% passing rate when considering the NIST battery of tests, thus it is impossible to found a difference between the two streams using the NIST suite.

2

) DieHARD battery of tests

The DieHARD battery of tests has been the most sophisticated standard for over a decade. Because of the stringent requirements in the DieHARD test suite, a generator passing DieHARD battery of tests can be considered good as a rule of thumb.

The DieHARD battery of tests consists of 18 different independent statistical tests. This collection of tests is based on assessing the randomness of bits comprising 32-bit integers obtained from a random number generator. Each test requires 2 23 32-bit integers in order to run the full set of tests. Most of the tests in DieHARD return a P -value, which should be uniform on [0, 1) if the input file contains truly independent random bits. Occasional P -values near 0 or 1, such as 0.0012 or 0.9983 can occur. However, an individual test is considered to be failed if its P -value approaches 1 closely, for instance, P -value > 0.9999. Results derived from applying the DieHARD battery of tests to the two random streams computed from experimental time series, reveals that both sequences can pass sucessfully all the tests. This confirms that the quantitative randomness of the bit stream taken from the chaotic laser intensity indicates similar features compared to the one obtained by the optoelectronic noise source (thus without any chaotic origine).

3) ENT test program ENT test program applies various tests to sequences of bytes stored in files and reports the results of those tests. The program is useful for evaluating random number generators for encryption and statistical sampling applications, compression algorithms, and other applications where the information density of a file is of interest [START_REF] Walker | Ent: a pseudorandom number sequence test program[END_REF].

There are 5 tests contained in the program:

1) Entropy test: Entropy in bits per character (or byte), which corresponds to the incompressibility of the sequence (as a perfectly random sequence cannot be compressed, since no part of it can be expressed in terms of other parts). Hence entropy of 8 bits/byte means perfect randomness in the sense of incompressibility. 2) χ 2 test: χ 2 testing is very common for goodness-of-fit of sample distributions of random numbers. It is known to be very sensitive to deficiencies in random number generators (when it is located between 5% to 95%, data are treated as random). 3) Sample test: Sample test means can be tested for bias in random number generation. In binary mode, the expected mean is 0.5 while for bytes, the expected mean is 127.5. 4) Monte Carlo test: a Monte Carlo approximation of π, which is simply the evaluation the area of the unit circle using the N generated random numbers (X i ,X i-1 ), i = 2, ..., N . 5) Serial Correlation test: Serial correlation coefficient evaluated from < X i , X i-1 > / < X i , X i >, for i = 2, ..N . The intended value for perfect random sequences is 0. In Table II, it is shown that the results for each pair of random streams, considering these five tests detailed above, are very closed one to each other. They all achieved to pass the threshold of the Chi-squared test, and the results are very similar for the other tests. To sum up, all these streams satisfy the same random-like behavior according to the ENT battery.

4) Comparative test parameters

Five well-known statistical tests [START_REF] Menezes | Handbook of applied cryptography[END_REF] are used too as simple comparison tools. They encompass frequency and autocorrelation tests. In what follows, s = s 0 , s 1 , s 2 , . . . , s n-1 denotes a binary sequence of length n. The question is to determine whether this sequence possesses some specific characteristics that a truly random sequence would be likely to exhibit. Standard tests intended to answer this question are Frequency test (monobit test), Serial test (2-bit test), Poker test, Runs test and Autocorrelation test (we refer the reader to [START_REF] Menezes | Handbook of applied cryptography[END_REF] or [START_REF] Bahi | Randomness quality of ci chaotic generators: Applications to internet security[END_REF] for detailed definitions).

We show in Table III a comparison between two random bits streams sourced respectively by the chaotic laser intensity and by the noisy optoelectronic signal. The results confirm that the proposed random streams present very closed statistical qualities. This finding implies that to have a chaos-like deterministic origin is not a required condition for high randomness quality in the proposed method.

Finally a comparison of the overall stability from 5 × 10 3 to 2 × 10 5 for these generators is given in Fig. 9. It can be seen that the trends for the amplitude movements of values are more or less in the same scale, which again indicates that all these random sequences share closed random properties.

V. DISCUSSION, CONCLUSION, AND PERSPECTIVES

Random number generation via photonic broadband signal generation does provide nowadays a novel and interesting approach allowing for unprecedented high bit rate of random bit streams. These photonics to digital world conversion are designed so that randomness can be certified according to most of the usual randomness tests such as NIST and DieHARD suites. Among the recently proposed physical systems and related processing intended to extract bit streams from photonic analogue waveforms, two rather different approaches can be identified: when the source of randomness explicitly stems from photonic noise [START_REF] Williams | Fast physical random number generator using amplified spontaneous emission[END_REF]- [START_REF] Wetzel | Random walks and random numbers from supercontinuum generation[END_REF], and when deterministic chaos is claimed to be at the origin of the random bit stream [START_REF] Uchida | Fast physical random bit generation with chaotic semiconductor lasers[END_REF], [START_REF] Reidler | Ultrahigh-speed random nuber generation based on a chaotic semiconductor laser[END_REF]. Whereas the first approach provides obviously, and by definition, a non deterministic random bit stream, the second one has an implicit potential source of determinism, similarly to the algorithmic and fully digital pseudorandom bit sequence RNGs. On the contrary to the photonic noise based RNGs, this indeed can be expected from the chaos-based photonic RNGs, since they also originates, at least partially, from deterministic dynamics, similarly to the algorithmic PRNGs.

In that particular context, we have proposed to address related issues, through the analysis of a particular post-processing method [START_REF] Reidler | Ultrahigh-speed random nuber generation based on a chaotic semiconductor laser[END_REF] applied in many chaos-based photonic RNGs. The actual role of deterministic chaos in this photonic RNG approach was known to be controversial [START_REF] Williams | Fast physical random number generator using amplified spontaneous emission[END_REF], [START_REF] Hirano | Fast random bit generation with bandwidthenhanced chaos in semiconductor lasers[END_REF] (but not analyzed), since purely non-deterministic (noisy) photonic signals were found to lead to similar randomness quality when using the same method.

In this article, we have checked on additional experiments that the proposed method can succesfully lead to high speed and high randomness quality, whether when used on strongly deterministic chaos provided by an electro-optical phase dynamics, or when used on a photonic noise having comparable spectral and statistical features. It is thus confirmed that chaos is not a necessary condition for the method to be successful. We have also analyzed and interpreted this digital postprocessing method. Standard signal theory arguments revealed that the method is actually acting as an entropy enhancement through aliasing phenomena, both with the time discretization Finally, we have also investigated how a typical signature of deterministic chaos, sensitivity to initial condition, can or cannot survive in the final bit stream. It is found that this determinism signature is actually lost precisely due to the random bit stream extraction method. We also conclude that the analyzed bit stream extraction method does not appear as suitable, when further use of deterministic feature is expected. For example, this happens when one wants to achieve synchronization between distant random bit stream (e.g. when cryptographic applications are concerned, such as one time pad cypher). Synchronizable random bit stream generated from deterministic chaos, is thus still an open problem, which would require other bit stream extraction methods able to provide both a strong enough determinism together with a high randomness quality.

Christophe Guyeux has taught mathematics and computer science in the Belfort-Montbéliard University Institute of Technologies (IUT-BM) this last decade. He has defended a computer science thesis dealing with security, chaos, and dynamical systems in 2010 under Jacques Bahi's leadership, and is now an associated professor in the computer science department of complex system (DISC), FEMTO-ST Institute, University of Franche-Comté. Since 2010, he has published two books, 9 articles in international journals, and 25 articles in peer reviewed international conferences dealing with security or chaos. Jacques M. Bahi was born in July, 25th 1961, after a Master of Science in Applied Mathematics, he received a Ph.D in Applied Mathematics from the university of Franche-Comté in 1991. From 1992 to 1999, he was a associate professor of applied mathematics at the Mathematical Laboratory of Besanc ¸on, his research interests were focused on parallel synchronous and asynchronous algorithms for differential-algebraic equations and singular systems. Since september 1999, he became a full professor of computer science at the University of Franche-Comté. He published about 150 articles in peer reviewed journal and international conferences and 2 scientific books. Jacques Bahi is a IEEE senior member since 2009. He is the head of the Distributed Numerical Algorithms team of the Computer Science Laboratory of Besanc ¸on, he supervised 21 PhD students. He is a member of the editorial board of 2 international journals and is regularly a member of the scientific commitees of many international conferences. Currently, he is interested in 1) high performance computing, 2) distributed numerical algorithms for ad-hoc and sensor networks and 3) dynamical systems with application to data hiding and privacy.

Fig. 1 :

 1 Fig. 1: Experimental setup of the optical system used to generate signal as physical random sources for the derivation of random bit sequences.

  Photodetector

Fig. 2 :

 2 Fig. 2: Scheme of the RNG using optical signal

Fig. 3 :

 3 Fig. 3: Illustration of the properly fulfilled sampling theorem conditions (a), and the incorrect sampling condition leading to aliasing (b).

Fig. 4 :

 4 Fig. 4: Statistics of recorded chaotic time series (experimental). (a) 1st-4th times DSD for 2.5GHz undersampling; (b) Zooming of the centering area of (a); (c) 1st-10th times DSD, log-log plot between weight and height, showing the hyperbolic relationship between the two as DSD is repetitively processed.

Fig. 5 :

 5 Fig. 5: Statistics of recorded noise time series (experimental). (a) 1st-4th times DSD for a 2.5GHz undersampling; (b) Zooming of the centering area of (a); (c) 1st-10th times DSD, the log-log plot between weight and height, showing the hyperbolic relationship between the two as DSD is repetitively processed.

Fig. 7 :

 7 Fig. 7: LSBs entropy evolution as a funtion of time, for an ensemble of simulated chaotic time series with the same initial conditions, but with different small noise added to it. The noise strength is -40 dB, with sampling rate 2.5 GHz 8-bit ADC. (a) For sampling value; (b) For the 1st DSD of the sampling value. (c) For the 2nd DSD of sampling value.

Fig. 8 :

 8 Fig. 8: Simulated averadge growth of bit entropy and its dependence on bit cell selecting (from MSB to LSB), without any DSD

  20GHz, 4th DSD, and 8 LSB method

Fig. 9 :

 9 Fig. 9: Overall Sequence Stability Comparison

TABLE I :

 I NIST SP 800-22 test results (P T )

	Method	2.5GHz, 1st DSD, 5LSB	20GHz, 4th DSD, 8LSB
	Source	Chaotic laser	Noise	Chaotic laser	Noise
	Frequency:	0.935716	0.798139	0.171867	0.834308
	BlockFrequency:	0.040108	0.350485	0.289667	0.867692
	CumulativeSums:	0.334152	0.575225	0.228927	0.688782
	Runs:	0.595549	0.834308	0.851383	0.637119
	LongestRun:	0.191687	0.964295	0.162606	0.304126
	Rank:	0.534146	0.037566	0.637119	0.719747
	FFT:	0.236810	0.514124	0.202268	0.249284
	NonOverlappingTemplate:	0.502510	0.491449	0.521769	0.501830
	OverlappingTemplate:	0.851383	0.964295	0.090936	0.574903
	Universal:	0.798139	0.739918	0.102526	0.319084
	ApproximateEntropy:	0.224821	0.236810	0.435436	0.419021
	RandomExcursions:	0.347389	0.229729	0.471174	0.104312
	RandomExcursionsVariant:	0.217344	0.209317	0.461569	0.350467
	Serial:	0.300289	0.366918	0.237996	0.606177
	LinearComplexity:	0.350485	0.262249	0.224821	0.935716

TABLE II :

 II ENT battery using 10 8 bits for each stream

	Method	Using source	Entropy	Chi-square	Sample	π error	Correlation
	2.5GHz, 1st	Chaotic laser	7.999984	67.18%	127.4988	0.03%	-0.000771
	DSD, 5 LSBs	Noisy signal	7.999986	7.13%	127.5034	0.03%	-0.000392
	20GHz, 4th	Chaotic laser	7.999986	73.48%	127.4973	0.03%	0.000481
	DSD, 8 LSBs	Noisy signal	7.999985	12.37%	127.5011	0.02%	-0.000411

TABLE III :

 III Comparison between the presented sources for a 2 × 10 7 bits sequence

	Subjects	Monobit	Serial	Poker	Runs	Autocorrelation
	Method		2.5GHz, 1st DSD, 5LSBs
	Chaotic laser	0.2509	1.9200	16.6650	16.6215	1.5739
	Noise	0.6019	0.7144	8.5606	17.5156	1.5247
	Method		20GHz, 4st DSD, 8LSBs
	Chaotic laser	1.4580	0.5199	13.1430	28.9460	1.1583
	Noise	0.2554	0.7835	14.0035	22.9136	1.6739
	(PRBS) generators.
	A major interest of the digital PRBS resides in their capability
	of generating a distant and synchronized random bit stream,
	which allows one to apply them in symmetric cryptography.
	A major advantage of PRNG is precisely their perfect de-
	terminism, and perfect control, due to their digital program-
	based generation process. This feature is also at the origin of
	their main drawback: the same absolute digital determinism
	can be used in principle for cryptanalysis, trying to guess the
	seed which can then deterministically and totally allow for the
	random bit sequence reproduction, even by an eavesdropper.
	Also from a more technical viewpoint, the processor based
	architecture of PRNGs defines some speed limitations related
	to the processor clock, and the number of elementary opera-
	tions needed to implement the PRNG algorithm. Noise based,
	and chaos based, photonic RNGs provide at least a technical

answer to the limited bit rate generation provided by purely algorithmic solutions. It was also reported in many attempts on photonics based RNGs, that high quality randomness is possible, since they can pass successfully all the standard NIST and DieHARD test suites. A strong open problem however still remains concerning the capability to control, and reproduce, the random bit stream provided by photonic chaos-based


Daniel Delahaye
email: delahaye@recherche.enac.fr

Stéphane Puechmorel
email: puechmor@recherche.enac.fr

Air Traffic Complexity based on Dynamical Systems

This paper presents a new air traffic complexity metric based on dynamical systems. Based on a set of radar observations (position and speed) a vector field interpolating these data is constructed. Once the field has been obtained, the Lyapunov spectrum of the associated dynamical system is computed on points evenly spaced on a spatial grid. The results of the computations are summarized on complexity maps, with high values indicating areas to avoid or to carefully monitor. A first approach based on linear dynamical system enable to compute an aggregate complexity metric. In order to produce complexity maps, two extensions of the previous approach have been developed (one in space and another in space and time). Finally, an approximation is proposed in order to localize the computation of the vector field by the mean of Local Linear Models.

I. INTRODUCTION

The ATM system has to cope with an increasing number of flights, pushing the capacity to its limits. As an example, the average daily traffic above Europe was 26286 flights/day, with a peak traffic demand in excess of 31000 flights [START_REF]Atfcm and capacity report[END_REF]. Although delays are kept low, it is expected from the same reference that capacity has to be extended in the future. Basically, two strategies can be devised : adapt the demand to capacity (slot-route allocation, collaborative decision making, ...) or adapt the capacity to the demand (Airspace design, 4D trajectory planning, autonomous aircraft, ...). The first approach can be used in the context of current ATM system, while innovative future designs will mainly follow the second strategy. Currently, complexity of the traffic is measured only as an operational capacity : the maximum number of aircraft that ATC controllers are willing to accept is fixed on a persector basis and complexity is assessed by comparing the real number of aircraft with the sector capacity. It must be noted that under some circumstances, controllers will accept aircraft beyond the capacity threshold while rejecting traffic at other times although the number of aircraft is well below the maximum capacity. This simple fact clearly show that capacity as a crude complexity metric is not enough by itself to fully account for the controller's workload. In order to better quantify the complexity, geometric features of the traffic have to be included. As previously stated, depending on the traffic structure, ATC controllers will perceive differently situations, even if the number of aircraft present in the sector is the same. Furthermore, exogenous parameters like the workload history can be influential on the perceived complexity at a given time (a long period of heavy load will tend to reduce the efficiency of a controller). Some reviews of complexity in ATC have been completed, mainly from the controller's workload point of view [START_REF] Hilburn | Cognitive complexity and air traffic control: A literature review[END_REF], [START_REF] Mogford | The complexity construct in air traffic control : A review and synthesis of the literature[END_REF], and have recognized that complexity is related to both the structure of the traffic and the geometry of the airspace. This tends to prove that controller's workload has two facets :

• An intrinsic complexity related to traffic structure.

• A human factor aspect related to the controller itself.

While most complexity metrics tend to capture those effects within a single aggregate indicator, the purpose of this work is to design a measure of intrinsic complexity only since it is the most relevant metric for an highly automated ATC system (no human factors). The first complexity indicator incorporating structural considerations along with the simple number of aircraft is the "Dynamic Density" of Laudeman et al. from NASA [START_REF] Laudeman | Dynamic density : an air traffic management metric[END_REF]. The "Dynamic Density" is a weighted sum of the traffic density (number of aircraft), the number of heading changes (> 15 degrees), the number of speed changes (>0.02 Mach), the number of altitude changes (>750 ft), the number of aircraft with 3-D Euclidean distance between 0-25 nautical miles, the number of conflicts predicted in 25-40 nautical miles. These factors are summed together using weighting factors that were determined by showing different traffic scenarios to several controllers. B.Sridhar from NASA [START_REF] Sridhar | Airspace complexity and its application in air traffic management[END_REF], has developed a model to predict the evolution of such a metric in the near future. Efforts to define "Dynamic Density" have identified the importance of a wide range of potential complexity factors, including structural considerations. However, the instantaneous position and speeds of the traffic itself does not appear to be enough to describe the total complexity associated with an airspace. A few previous studies have attempted to include structural consideration in complexity metrics, but have done so only to a restricted degree. For example, the Wyndemere Corporation proposed a metric that included a term based on the relationship between aircraft headings and dominant geometric axis in a sector [START_REF] Inc | An evaluation of air traffic control complexity[END_REF]. The importance of including structural consideration has been explicitly identified in work at Eurocontrol. In a study to identify complexity factors using judgment analysis, Airspace Design was identified as the second most important factor behind traffic volume [START_REF] Kirwan | Investigating complexity factors in u.k. air traffic management[END_REF]. Histon,et,al. [8], [START_REF] Histon | Structural consideration and cognitive complexity in air traffic control[END_REF] investigated how this structure can be used to support structure-based abstractions that controllers appear to use to simplify traffic situations. The previous models do not take into account the intrinsic traffic disorder which is related to the complexity. The first efforts related with disorder can be found in [START_REF] Delahaye | Air traffic complexity : Towards intrinsic metrics[END_REF]. This paper introduces two classes of metrics which measure the disorder of a traffic pattern.G.Aigoin has extended and refined the geometrical class by using a cluster based analysis [START_REF] Aigoin | Air traffic complexity modeling[END_REF]. All the previous metrics capture only one feature of the complexity and are not able to produce an aggregate metric which can capture all the possible situations (high-low density, how-low convergence, translation organization, round about organization etc ...). The work presented in this paper is based on dynamical systems modeling of the air traffic. A dynamical system describes the evolution of a given state vector. If such a vector is given by the position of aircraft X =[ x, y, z] T , a dynamical system associates a speed vector ˙ X =[ v x ,v y ,v z ] T to each point in the airspace. The key idea is to find a dynamical system which models the observed aircraft trajectories. Based on this dynamical system modeling, a trajectory disorder metric can be computed.

The first part of this paper will present a linear dynamical system modeling for which the complexity metric can be represented into a complex coordinate system. In this system, it is very easy to identify any speed vector organization pattern. The second part introduces a non linear extension of the previous dynamical system modeling that may model any observed traffic. This extension is done in space only. Such a non linear modeling can be used to produce maps of traffic complexity by identifying areas with high(low) complexity. The third part presents the time extension of the previous model for which works straightly on the trajectory segments instead of speed vectors. Finally, the fourth part, propose an approximation approach based on Local Linear Models for computation improvement. This last model is also able to take into account some uncertainties on aircraft positions in order to improve the robustness of complexity maps.

II. LINEAR DYNAMICAL SYSTEM MODELING

The key idea of this approach is to model a set of aircraft trajectories by a linear dynamical system which is defined by the following equation :

˙ X = A. X + B ( 1
)
where X is the state vector of the system :

X = xyz T (2)
This equation associates a vector speed ˙ X to a position in the space coordinate X. The eigenvalues of the matrix A control the evolution of the system. The real part of those eigenvalues is related with the convergence or the divergence property of the system in the direction of the eigenvector. When such a eigenvalue has positive real part, the system is in expansion mode and when it is negative the system is in contraction mode. On the other end, the imaginary part of the eigenvalues are related with the level of curl organization of the system. Depending of those eigenvalues, a dynamical system can evolve in contraction, expansion, rotation or a combination of those three modes. Eigenvalues loci for several typical situations.The small squares are the initial positions of aircraft at a given time(this represents the observation given by a radar for instance with the associated speed vector). associated to the matrices A are given on the bottom part of figure 1. When the relative distances between aircraft remain unchanged with time (situation 1 and 4), it can be noticed that the real parts of the eigenvalues of matrix A are zero. When the norms of the relative distances between aircraft diminish with time (situation 2) the real part of the eigenvalues are negatives; finally, when those relative distances increase with time (situation 3) the real part of the eigenvalue are positive. The largest those real part are in absolute value, the more those evolutions are fast. When more aircraft are involved, a Least Mean Square regression is applied in order to extract matrix A and vector B. Based on the observations of aircraft (positions and speed vectors), the dynamical system has to be adjusted with the minimum error. This fitting has been done with a Least Square Minimization method. For each considered aircraft i, it is supposed that position X i = [x i ,y i ,z i ] T and speed vector V i =[vx i ,vy i ,vz i ] T are given. An error criterion between the dynamical system model and the observations is computed :

E = i=N i=1 V i -A. X i + B 2
In order to use matrix forms the following matrices are introduced.

X = ⎡ ⎢ ⎢ ⎣ x 1 ... x N y 1 ... y N z 1 ... z N 1 ... 1 ⎤ ⎥ ⎥ ⎦ V = ⎡ ⎣ vx 1 ... vx N vy 1 ... vy N vz 1 ... vz N ⎤ ⎦ C = ⎡ ⎣ a 11 a 12 a 13 b 1 a 21 a 22 a 23 b 2 a 31 a 32 a 33 b 3 ⎤ ⎦
With such matrices, the error criterion E can be written as :

E = V -C.X .
To minimize E is the same as to minimize

E 2 = V -C.X 2 .
The problem is to find the matrix C which minimize such a criterion. By using the matrix derivation properties, the gradient of E 2 is given by :

∇ C E 2 = -2.(V -C.X).X T .
The gradient is then canceled :

∇ C E 2 =0⇔ C.X.X T = VX T .
The optimum matrix C opt is then given by :C opt = V.X T .(X.X T ) -1 . The expression X T .(X.X T ) -1 is the pseudo inverse of the matrix X T and can be written as : X T .(X.X T ) -1 = L T .S -1 .R. where S is the diagonal matrix of the singular values of the matrix X T . This singular value decomposition is very helpful in order to avoid condition troubles. The matrix C is finally given by :C = V.L T .S -1 .R. The matrix A is then extracted from the matrix C, and the associated eigenvalue decomposition is given by A = L.D.U T . This first model enables to produce an aggregate metric associated to any traffic situation and can recognize any global organization pattern. It gives an aggregate measure of a given traffic situation but is not able to identify high (low) complexity areas in the airspace. For such purpose, we have developed a complexity metric based on non-linear dynamical system.

III. NON LINEAR EXTENSION IN SPACE

A. Non linear Dynamical System

A non linear dynamical system is summarized by the following equation :

˙ X(t)= f ( X) ( 3
)
where X is the state vector of the system ( X =[ x, y, z] T ) and f : C 2 vector field, describe systems which integral curves may fit the observed trajectories. This equation associates a vector speed ˙ X to a position in the space coordinate X and then synthesis a particular vector field. Based on the observations of the aircraft (positions, speed vectors), the dynamical system has to be adjusted with the minimum error.This fitting is done with a Least Square Minimization (LMS) method for which the following criterion is used :

E 1 = i=N i=1 V i -f ( X i ) 2 (4)
where N is number of observations. If we consider criterion E 1 only, it can be shown that there is an infinite number of vector fields f which can be adjusted to the observations. In order to keep the smoothest one, another criterion is added which has to be minimized, the so-called "div-curl" criterion :

E 2 = R 3 α ∇div f ( X) 2 + β ∇curl f ( X) 2 d X (5)
with α, β positive weights controlling the smoothness of the approximation by focusing on constant divergence or constant curl. In the following, we will consider α = β = 0.5; in such case :

E 2 = ∆ f ( X) 2
where ∆ f ( X) is the Laplacian of the vector field f . The joint minimization of E 1 and E 2 induces a unique optimum vector field [START_REF] Amodei | A vector spline approximation[END_REF], [START_REF] Suter | Fast evaluation of vector splines in three dimensions[END_REF] :

f ( X)= N i=1 φ( X -X i ). a i + A. X + B 0 V γ γ 0 x x (t,x ) (t,x) γ D(t, )
Fig. 2. Time evolution of a reference trajectory and a perturbed trajectory with a i parameter vectors (one for each observation),

φ( X -X i )=Q( X -X i 3 )
and

Q = ⎡ ⎣ γ 00 0 γ 0 00γ ⎤ ⎦ with γ = ∂ 2 xx + ∂ 2 yy + ∂ 2 zz
The resulting adjustment is done without error (⇒ minE 1 =0).

When α = β =0 .5, the vector spline function φ has the following structure :

φ( X -X i )=12. X -X i (r i = X i -X ).
It must be noticed, that farthest observations has more weight in such calculation. Then, in order to compute the smoothest vector field which fit exactly the measures, all observations have to be taken into account in the computation.

IV. LYAPUNOV EXPONENTS

The metric chosen for complexity computation relies on a measure of sensitivity to initial conditions of the underlying dynamical system called Lyapunov exponents. In order to figure out what Lyapunov exponents are, let consider a point and look at its evolution when transported by the dynamical system. Let x 0 be fixed (initial point) and let γ be a point trajectory of the dynamical system associated to the vector field f given by :

γ(t, x 0 )= x 0 + t 0 f (u, γ(u, x))du (6)
Assume now that trajectory is disturbed by a small perturbation ǫ,w eh a v e:

γ(t, x 0 + ǫ)=γ(t, x 0 )+∇ x f (γ(t, x)). ǫ + o( ǫ )
where ∇ x f (t, γ(t, x)) is the differential of the vector field f at x. Divergence to nominal trajectory with respect to time is thus γ(t, x 0 )γ(t, x) = D(t, s) (see figure 2). γ(t, x) being defined as a flow :

∂γ(t, x) ∂t = f (t, γ(t, x)) γ(0, x)= x
2071 with f a smooth vector field, it is possible to show that D(t) satisfies a differential equation also. Given a nominal trajectory γ(t, x 0 ), then divergence of nearby trajectories can be found up to order one in xx 0 by solving :

∂D(t, x) ∂t = ∇ x f (t, γ(t, x)).D(t, x) D(0, x)= x -x 0
If the three space dimensions are considered (x, y, z), and since the previous equation is linear, it can be extended to the matrix form :

dM (t) dt = ∇ x f (t, γ(t, x)).M (t) M (0) = Id
Where each column of the M matrix corresponds to the divergence associated to the principal coordinate axis. This equation is called the variational equation of the flow.

The variational equation describes in some sense a linear dynamical system "tangent" to the original one. Let U t (t)Σ(t)V (t)=M (t) be the singular value decomposition of M (t). The Lyapunov exponents are mean values of the logarithms of the diagonal elements of Σ(t):

κ(s)=- 1 T T 0 log(Σ ii (t))dt ∀Σ ii (t) ≤ 1
When Lyapunov exponents are high, the trajectory of a point under the action of the dynamical system is very sensitive to initial conditions (or to the parameters on which the vector field may depend), so that, situation in the future is unpredictable. On the other hand, small values of the Lyapunov exponents mean that the future is highly predictable (expected to be comfortable for a controller).

So, the Lyapunov exponent map determines the area where the underlying dynamical system is organized. It identifies the places where the relative distances between aircraft do not change with time (low real value) and the ones where such distance change a lot (high real value).

Let us now describe the practical algorithm for computing complexity maps. 1) Dynamical regression based on radar observations.

2) The vector field and its associated gradient is then computed on each point of a cube of airspace. 3) Lyapunov exponents are computed at each point of the cube by the mean of a Runge-Kutta integration. The figure 3 shows an example of Lyapunov exponents map for which full organized miles in trail trajectories (from south west to north east) cross two random traffic situations. This figure shows clearly a complexity valley on the miles in trail direction. This organization may have been detected even if the miles in trail trajectories would have been structured on a curve trajectory. That is the strong point of this metric: Lyapunov exponents are able to identify any kind of trajectory organization (when aircraft follow the same trajectory at the same speed).

This approach is able to produce a picture of the complexity at a given time. If we want to be able to predict the Fig. 3. Miles in trail traffic between disordered areas complexity in a near future, one must be able to take into account time dimension in the model. This is the topic of the following section.

V. N ON LINEAR EXTENSION IN SPACE AND TIME

The previous extension, is able to model air traffic at a given time but is not valid for observation with large time extension.

In such case, the vector field f has to be time dependent :

˙ X = f ( X,t)
As before we are looking for a vector field which both minimize the error criterion E 1

minE 1 = i=N i=1 k=K k=1 V i (t k ) -f ( X i ,t k ) 2
and the smoothest criterion E 2 :

minE 2 = R 3 t ∆ f ( x) 2 + ∂ f ∂t 2 d xdt
The differential operator associated with this variational problem is given by : P = ∂ 2 ∂t 2 -μ∆ 2 An elementary solution can be found in S ′ as :

fp 1 √ μ ξ 2 e -|t| √ μ|ξ| 2
where the previous expression is the Fourier transform in S ′ , the space of rapidly decreasing smooth functions, with respect to the spatial coordinates of the elementary solution (fp stand for "Finite part" in the Hadamard sense). It turns out that this distribution is indeed obtained from a L 1 loc (integrable on compact sets) mapping so that the inverse Fourier transform can be obtained readily by integration as :

p(t, x)= 1 8π 3 R 3 1 √ μ ξ 2 e (-|t| √ μ ξ
2 ) e (i x,ξ ) dξ by Fubini's theorem and polar change of variables it becomes :

p(t, x)= 1 8π 3 √ μ R e (-|t| √ μr 2 S 2
e (i x,rs ) dσ(s)dr

2072

with dσ the solid angle measure. Using a polar parametrization of the unit sphere :

p(t, x)= 1 8π 3 √ μ R e (-|t| √ μr 2 )
2π 0 π 0 e (i x r cos θ) sin θdθdφdr and finally :

p(t, x)= 1 2π 2 √ μ R
e (-|t| √ μr 2 ) sin x r x r dr by Parseval equality :

p(t, x)= 1 4 x π 2 π |t| √ μ x 2π - x 2π e " -π 2 ω 2 |t| √ μ "
dω so in terms of error function erf :

p(t, x)= 1 4π 2 x erf x 2 |t| √ μ
The exact solution for such joint minimization problem is given by the following equation :

f ( X,t)= N i=1 K k=1 ψ( X(t)-X i (t k ) , |t-t k |). a i,k +A. X+ B with ψ(r, t)=diag σ r. √ π .erf r σ . 1
β = ψ(r, t) r 2 r = x 2 + y 2 + z 2
The algorithm used to produce complexity map is the same as in the space extension case unless, here, the field and the associated gradient are time-dependent. As, before, one need the contribution of all aircraft in order to compute the vector field and its associated gradient at a given point. When few observations are involved (≃ 2000) the previous algorithms are tractable. When more observations are considered, one can improve the regression phase by the mean of Krylov spaces trick [START_REF] Gutknecht | A brief introduction to krylov space methods for solving linear systems[END_REF], [START_REF] Bjorck | Numerical methods for least squares problems[END_REF] but the vector field reconstruction and gradient computation become very slow.

Another limitation of the previous approaches is linked to the robustness of the produced complexity maps. When such tools are applied to a given traffic situation, the produced maps measure the current complexity. When uncertainties has to be taken into account (uncertainties in the future position of aircraft or when we want to produce complexity map for the entire airspace in order to compare different countries (in such case flight plans are simulated with relevant noise due to the wind, flight delay, etc...)), one must use an approximation of the previous exact models for which local linear models are used.

VI. APPROXIMATION BY LOCAL LINEAR MODELS

In order to speed up the computation of the vector field, we have developed a local approximation for which only local observations are needed for the regression. Before applying such approximation, the observations are first normalized according to the global linear model (A. X + B) by subtracting the linear part to the observations :

v i = V i -(A. X i + B)
Then, such approximation will locally update the linear regression model in order to fit observations. The first order space-time linear expansion of a vector field f (t, X) is given by the following expression :

f (t 0 , X 0 )= f (t, X)+ ∂ f (t, X) ∂t (t 0 -t)+ + ∂ f (t, X) X ( X 0 -X)+O(|t 0 -t| + X 0 -X )
where ∂ f (t, X) ∂t is the time derivative vector of the vector field f and + ∂ f (t, X) X is the associated space derivatives matrix. This equation represents a local linear model of the vector field f in the vicinity of (t, X). This approximation is only valid near a given point (t, X). We propose to use this model in order to compute the vector field base on a given set of local observations. First, a point (t, X) in the state space is considered (a grid point) for which observations are searched in its vicinity. The field at this point ( f (t, X))i s then computed in order to minimize the errors between local observations v i (t i , X i ) and the linear expansion of the vector field f at the grid point (t, X) v i (t i , X i ) ≃ f (t, X)+ ∂ f (t, X) ∂t (t it)+ ∂ f (t, X) X ( X i -X)

= a + b.(t it)+C.( X i -X)

where a = f (t, X), b = ∂ f (t, X) ∂t and [C]= ∂ f (t, X)

X

We are seeking vectors a, b and matrix [C] which minimize the error criterion J : min a, b,C J = N i=1 v i ( X i ,t i )a + b(t it)+C.( X i -X) 2 .ψ i (7) where ψ i = ψ(t i -t, X i -X) is a space-time weight function which select the closest observations. If we note

X = ⎡ ⎢ ⎢ ⎢ ⎢ ⎣ 1 ...

  Fig.1. Eigenvalues loci for several typical situations.The small squares are the initial positions of aircraft at a given time(this represents the observation given by a radar for instance with the associated speed vector).

1 (t 1 - 1

 111 t) ... (t Nt) (x 1x) ... (x Nx) (y 1y) ... (y Ny) (z 1z) ... (z Nz) ... vx N vy 1 ... vy N vz 1 ... vz N ⎤ ⎦ Ψ= ⎡ ⎣ ψ 1 ... 0 ... ... ... 0 ... ψ N ⎤ ⎦ M = ⎡ ⎣ a x b x C xx C xy C xz a y b y C yx C yy C yz a z b z C zx C zy C zz ⎤ ⎦

49th IEEE Conference on Decision and Control December 15-17, 2010 Hilton Atlanta Hotel, Atlanta, GA, USA 978-1-4244-7746-3/10/$26.00 ©2010 IEEE

The criterion J can be written as :

By differentiating J with respect to M , we find : M opt = V.Ψ.X T (X.ΨX T ) -1 In order to avoid condition trouble, a singular value decomposition will be used to compute such matrix. We then propose the following algorithm to compute complexity maps :

A. Stochastic extension

When countries have to be compared in term of airspace, one must consider many days of traffic with regular flight plans for which observations are sample of a random process. In this framework, positions of aircraft on their trajectories undergo random deviations due to wind, delay, etc ... Even when one want to predict complexity in a given area, one must be able to take into account uncertainties of future positions. For both applications, time shift random process is the most relevant model, for which the time of observation t i is replaced by t i -τi , where τi is a random variables (all random variable will be note with a upper tilde .).

Then t i is now a random variable : ti = t i -τi We have

, where γ(t) is the aircraft trajectory. The criterion J becomes a random variable and we are seeking a, b, C in order to minimize the expected value of J : andE[ Ψ] the expected value of the random matrices Ṽ , X, Ψ the criterion E[ J] can be written as :

By differentiating E[ J] with respect to M , we find :

The same algorithm as the one propose in the deterministic case will be used for computing complexity maps.

VII. CONCLUSION

We have presented in this paper a new air traffic complexity metric based on vector field model of air traffic. Most of the work has been devoted to improvements in interpolating splines used to fit a vector field to the observation. Unlike linear models which produce mean complexity indicators, the non-linear one may give local information, thus providing a way of displaying maps of complexity. A local approximation has also been proposed in order to parallelize the computation and to improve the robustness of the map produced when time shift noise is included in observations.

In a future work, such a tool will be applied to comparison of US and Europe airspace by producing complexity map of USA and Europe for a full day of traffic.