Prof Guangkuo Li

Dr Yuhua Li

Dr Feng Ye

Qian Li

Shao Hao Wang

Dr Benjamin Wetzel

Dr Roy Davidson

Brent E Little

Sai Tak Chu
email: saitchu@cityu.edu.hk

T Sai

Chu

Prof Hao Shao

Wang

Supercontinuum Generation in Dispersion Engineered Highly Doped Silica Glass Waveguides

Keywords: supercontinuum, highly doped silica glass, dispersion engineering, slot waveguide

come   L'archive ouverte pluridisciplinaire

Introduction

Supercontinuum (SC) generation arises from the propagation of a high-power optical pump through a nonlinear medium. This highly complex dynamical process involves both the nonlinearity and dispersion of the medium, and can be extremely sensitive to these characteristics as well as those of the input pump signal. Due to the variety of SC broadening scenario and associated output spectral properties [1], SC has found numerous applications as ultra-broadband light sources in many fields such as frequency metrology [2,3], spectroscopy [4,5], and optical coherence tomography [6,7]. In addition, SC broadening dynamics have also been investigated from a more fundamental viewpoint, being a testbed for the study of exotic phenomena such as rogue waves [8,9] and dark solitons [10,11], or even exploited as an onset nonlinear broadening stage in cascaded architectures towards ultrashort UV or mid-IR pulse formation [1].

The introduction of photonic crystal fibers (PCF) has made these light sources more accessible as the group velocity dispersion (GVD) of PCFs can be engineered to enable strong nonlinear interactions over significant fiber propagation lengths [12]. Such dispersion engineering lowers the energy requirement for the nonlinear spectral broadening processes, allowing the replacement of bulky high power OPO lasers by cost-effective semiconductor diode lasers to serve as the optical pump [13,14]. Besides optical fibers, SC can also be generated in integrated nonlinear optical waveguides, which have been demonstrated in platforms spanning silicon [15][16][17][18], silicon nitride [19][20][21], silicon germanium [22][23][24], chalcogenide [25][26][27], aluminum nitride [28,29], InGaP-on-insulator [30] and highly doped silica glass (HDSG) [31]. These waveguides generally exhibit higher nonlinearity and tighter mode confinement compared to optical fibers. The generation of the SC spectrum in an integrated fashion allows for the direct connection of the broadband light source with other waveguide components on the same chip.

The advantage is two-fold: the integrated approach not only drastically reduces the overall footprint of the device and the overall system, but it can also provide direct access to a host of optical functionalities compatible with advanced and reconfigurable optical pulse processing within a single (stable) on-chip photonic framework [32,33].

Among those integrated architectures, HDSG constitutes a versatile CMOS-compatible platform that possesses very low linear and nonlinear losses. The Kerr nonlinearity of n2 = 1.15 × 10 -19 m 2 W -1 with nonlinear coefficient γ is approximately two hundred times larger than the single mode fibers (SMF). Its high nonlinear figure-of-merit [34] makes the platform highly efficient in nonlinear applications [35][36][37][38]. In 2010, Duchesne et. al., demonstrated SC generation in the HDSG spiral waveguide using a 45 cm long channel waveguide with core dimension of 1.45 μm by 1.50 μm and with a refractive index of n = 1.7 [31]. Using a 100 fs pulse with 1450 W peak power to pump the waveguide anomalous dispersion region at 1550 nm, they obtained a SC spectral bandwidth of 300 nm, where the SC spectrum coverage was limited by the dispersion profile of the waveguide. In fact, the span of the obtained SC spectrum is comparable to only the anomalous dispersion region of the waveguide. To boost SC spectral broadening dynamics, it is generally desirable to design waveguide dispersion profiles that are flat across a wide spectral range, and ideally featured with low anomalous dispersion. The flat and low anomalous GVD region indeed limits the temporal walk-off (and eventual pulse spreading) during SC formation, thus maximizing the phase matching and enhance the nonlinear spectral broadening throughout propagation [39]. Furthermore, such a dispersion engineering approach can also be leveraged to adjust the zero-dispersion wavelength (ZDW) location and therefore tailor the waveguide for an optimized pump wavelength [40].

In this work, we report a systematic study of the role and impact of a thin oxide layer inclusion within the core of HDSG channel waveguides for SC generation. In particular, we analyse how the thickness of the oxide layer can be used to engineer the dispersion profile of the waveguide and the corresponding output SC spectrum. We show that using a 0.1 μm oxide slot extends the SC spectrum coverage to over 1.5 octaves, which is 664 nm wider compared to SC generated from the conventional channel waveguide without the slot [31]. We also perform numerical simulations based on a generalized nonlinear Schrödinger equation (GNLSE) model to study soliton dynamics, spectral broadening processes and SC generation in the fabricated HDSG waveguides. We show that numerical results are in good agreement with experimental measurement, and that the generated supercontinuum exhibits both phase and intensity stability attested by a high coherence degree and a high signal-to-noise ratio (SNR) across its entire spectral bandwidth.

Dispersion Engineering

Dispersion plays a critical role in the spectral broadening of optical pump signals in nonlinear media. For instance, soliton existence and propagation physically require anomalous GVD waveguides to produce a negative chirp able to balance the positive chirp produced by selfphase modulation (SPM). In contrast, in the normal GVD region, the pulse is stretched due to the uncompensated SPM positive chirp, which in turn limits its peak power and associated SPM efficiency. Thus, a wide anomalous GVD region with a flat and close-to-zero dispersion can reduce the walk-off between the high and low frequency components and strengthen the nonlinear effects and further gradually compress the pulse width and broaden the spectral span of SC correspondingly. A popular waveguide design that allows flexible engineering of the waveguide dispersion is the slot waveguide structure. It has a thin layer of lower index oxide, the slot, placed in the waveguide core. This type of structure was shown effective in engineering the dispersion profile of the silicon rich waveguides and greatly expanded their SC spectrum [41][42][43]. The addition of the slot transforms the conventional channel waveguide into a coupled structure with the symmetric and asymmetric super-modes [44,45], where the symmetric mode exhibits a concave dispersion that can be used to help flatten the convex shape of the channel waveguide dispersion. Although the symmetric mode exists in both polarizations, the TM mode with the field normal to the slot interfaces does have a much higher dispersion flattening effect due to the intensity enhancement at the slot interfaces. Hence, dispersion flattening utilizing the TM mode is more common in the silicon, silicon nitride and silicon germanium platforms because of their high refractive indices, where the dispersion of the conventional waveguide is highly anomalous [44,46]. Compared to the Si-rich platforms with typical index contrasts of over 50%, the HDSG platform has an index contrast of only 17% and requires a relatively lower dispersion compensation for the flattening effect than the Si-rich platforms. In addition, the compensation from the TE mode is sufficient to provide a flat and low dispersion profile. The effectiveness of the slot structure will be investigated in the dispersion engineering of the HDSG waveguides in this study. (e), respectively. For TE polarization, the slot splits the power into the upper and lower core regions, where most of the energy is confined. In contrast, the bulk of the TM mode power lies in the vicinity of the slot region due to the discontinuity of the E-field at the slot interfaces. The influences of the slot thickness and the waveguide width on the dispersion of the waveguide structure with total thickness of 1.5 μm are shown in the calculated dispersion profiles in Figure 2. The SEM image is used to extract the waveguide sidewall angle of 4.1º in the calculated profiles. Figures 2(a) and(b) show the effect of the slot thickness on the TE and TM dispersions, respectively. It can be seen that the increase of the slot thickness induces a downward shift to the dispersion curve for both the TE and TM modes and effectively moves the ZDW location of the TE dispersion curve toward the longer wavelength end. Here, the ZDW is shifted from 1060 nm to 1620 nm by the addition of a thin layer of 0.2 μm thick oxide slot. For the TM case, with the dispersion of the waveguide without the slot is already very near zero, the dispersion curve can be moved entirely into the normal dispersion regime with the introduction of a thin ~0.05 μm thick slot. The comparison of the downward shifts introduced by the slot in Figures 2(a) and (b) shows that the amount of dispersion compensation is highly polarization dependent, with the TM curve shifted more than the TE curve caused by the field enhancement of the normal field at the slot interface.

The effect of the waveguide width on the TE and TM dispersions is shown in Figures 2(c) and (d), respectively, and it shows a different signature than the influence on the dispersion by the slot thickness. Although decreasing the waveguide width can also lower the dispersion curve, it only applies to the wavelengths after the elbow of the dispersion curve at λ > ~1000 nm, while keeping the dispersion below the elbow almost the same. As can be seen in the calculated integrated dispersion in Section 3, the dispersive wave (DW) is located just before the elbow of the dispersion curve, hence changing the waveguide width has less influence on the location of DW, but it does affect the soliton dynamics which depend on the dispersion characteristics at the longer wavelengths. It is worth mentioning that one can also tailor the dispersion of the channel waveguide by making adjustment to the waveguide dimension. Figures 3(a and ts = 0.1 μm, and the channel waveguides with similar waveguide dimension, with changing only the width or the thickness, respectively. Like the situation in Fig. 2(c), a moderate amount of change in either the width or the thickness of the channel waveguide will only alter the dispersion profiles at the longer wavelength side of the dispersion curve, while leaving the dispersion at below the elbow relatively the same. Hence, the slot waveguide design is a more effective method for the engineering of the specific properties of the dispersion, where the slot dimension can be used to determine the location of the DW, and the waveguide width can be used to fine-tune the soliton dynamic.

In our proposed structure with w = 1.8 μm, t = 1.5 μm and ts = 0.1 μm, the TE mode has a very low anomalous dispersion of below 43 ps km -1 nm -1 from 1260 nm to 2300 nm, and very flat dispersion with variation of less than 6.9 ps km -1 nm -1 from 1550 nm to 2300 nm. Although the slot waveguide structure provides great flexibility in the engineering of the dispersion profile, the added slot layer increases the propagation losses by 0.06 dB cm -1 for the TE and by 0.08 dB cm -1 for the TM modes.

Results and Discussions

Supercontinuum measurements

The pump source used in the experiment was a femtosecond pulse fiber laser (Amonics AFL-1560-FS-300-03-AOCP-SA) centered at 1560 nm with full width at half maximum (FWHM) pulse duration adjustable from 180 fs to 293 fs, and a corresponding to pulse energy Ep ranging from 621 pJ to 280 pJ at repetition rate of 40 MHz. The pulse duration was retrieved from GRENOUILLE measurements (Model 15-100-USB) [START_REF] Trebino | Frequency-Resolved Optical Gating: The Measurement of Ultrashort Laser Pulses[END_REF]. Separate optical spectrum analyzers (OSA) were used to record the SC spectra in the experiment (YOKOGAWA AQ6373B for 500 to 800 nm, AQ6370D for 800 to 1600 nm and AQ6375 for 1600 to 2300 nm). Figures 4(a-c) show the measured spectral evolution within three waveguides having the same length of 3.1 cm but featured with different dispersion profiles. At 1560 nm, the pump wavelength is in the anomalous GVD region for the waveguides with ts = 0 and 0.1 μm, and in the normal GVD region for the ts = 0.2 μm waveguide. When pumped in anomalous GVD region, as seen in Figures 4(a,b), the spectra experience a similar slow broadening process under the effect of SPM when the energy Ep is below 400 pJ. Although both spectra experience a stronger broadening for Ep > 400 pJ, the spectrum for ts = 0.1 μm spreads at a much faster rate than the spectrum for ts = 0 μm due to its lower dispersion. In fact, at Ep = 450 pJ, the spectrum for the waveguide with a slot of 0.1 μm had spread over 1200 nm, close to the maximum bandwidth reached in this experiment. In comparison, for the waveguide without slot, the spectrum had extended to approximately 500 nm and continue to spread to higher energy values of Ep. For these two cases with the pump in the anomalous GVD region, soliton fission occurs when Ep is between 380 pJ to 420 pJ, such high order soliton dynamics play a key role in the expansion of the spectra after the initial SPM stage. When the pulse is compressed to its Fourier limit under the combined action of GVD and SPM, higher-order dispersion and symmetry breaking effects (e.g. nonlinear self-steepening) perturb the pulse, inducing soliton fission where individual solitons are ejected from the initial pulse and are shifted towards longer wavelengths via soliton self-frequency-shift (SSFS) [START_REF] Lee | [END_REF]. Here, each ultrashort fundamental or low-order soliton carries a large number of different frequency components to form a broadband SC spectrum. Conversely, through their spectral shift and associated reshaping [START_REF]Nonlinear Fiber Optics[END_REF], the ejected solitons shed part of their energy to a linear dispersive wave (DW) in the normal GVD region, where the phase matches with the soliton in anomalous GVD region. For the ts = 0 μm waveguide, because of the large distance between the pump and the normal GVD region, the generated SC spectrum is mostly confined within the anomalous region, where the DW only starts to appear when Ep reaches a relatively high value of 500 pJ (and is detached from the main body of the spectrum). In contrast, the addition of the slot flattens the dispersion curve of the ts = 0.1 μm waveguide allowing the SC spectrum to spread across the anomalous region into the normal region to generate a large DW peak. In this case, DW generation begins shortly after soliton fission and the resulting DW peak is connected to the central part of the SC spectrum. As for the broadening of the ts = 0.2 μm waveguide, illustrated in Figure 4(c), it is mainly due to the SPM and is much narrower than the previous two cases. With the pump in the normal GVD region, pulse propagation cannot induce any high order soliton dynamics nor satisfying the required phase matching condition, so that no soliton fission nor DW generation is triggered.

Figures 4(d-f) show the input and output spectra for the three waveguides pumped at Ep = 621 pJ, with 180 fs FWHM pulse duration and 2835 W peak power. In Figure 4(a), the spectrum of the ts = 0 μm waveguide with 30 dB bandwidth spans from 1210 nm to 1912 nm with the DW located at 620 nm, whereas the dispersion engineered ts = 0.1 μm waveguide exhibits a SC spectrum with a 664 nm wider bandwidth, shown in Figure 4(b), spanning from 817 nm to 2183 nm, with the DW at 900 nm. For the ts = 0.2 μm waveguide, shown in Figure 4(c), it produces a classical symmetrical spectrum between 1335 and 1796 nm that is due to SPM broadening.

Supercontinuum simulations

Numerical simulations based on the GNLSE model [12] were performed to investigate the spectro-temporal behavior of the SC generation process in these waveguides. Both the GVD and higher-order dispersion are accounted for in the numerical model. We determined the nonlinear parameters of the waveguides using the rigorous definition from [START_REF] Sato | [END_REF]. It yields γ = 0.170, 0.148 and 0.135 m -1 W -1 for the TE mode of the ts = 0, 0.1 and 0.2 μm waveguides, respectively. This shows that the introduction of the oxide slot has slightly reduced the waveguide nonlinearity. Figures 5(a-c) show the simulated SC spectra using the input pulse shape retrieved from the experimental measurement of Figures 4(a-c). Except for the spectral region above 2000 nm, the simulated results agree well with the measurements and are able to replicate most of the features observed in the experimental spectra, including the locations of the DW for the ts = 0 and 0.1 μm waveguides. The discrepancies in the longer wavelength region are mainly due to HDSG possessing a strong absorption peak at 2030 nm, which leads to the appearance of artifacts ("peaks" at 2100 nm) in both ts = 0 and 0.1 μm waveguides. The measured DW positions of these waveguides coincide very well with the phase mismatch βint = 0 shown in the upper panels of Figures 5(a,b):

where β(ω) is the propagation constant, ωs is the soliton center frequency, and vg is the group velocity [40]. We note that the addition of a lower index slot has shifted the phase matching position toward a longer wavelength, thus easing the DW generation process for ts = 0.1 μm waveguide. The simulated pulse evolution within the three waveguide structures are shown in Figures 6(ac). The simulations were carried out for propagation distances up to 6 cm (approximately twice the actual length of the waveguides in the experiment) to fully illustrate the soliton dynamics and the DW generation process. In general, soliton dynamics are better understood from the perspective of feature lengths, where the input pulse can be regarded as a high order soliton with a soliton order N = (LD/LNL) 0.5 . Here, LD is the dispersion length obtained as LD = T0 2 /|β2| where T0 is the half width at 1/e intensity, and LNL is the nonlinear length where LNL = 1/γP0. The soliton fission length Lfiss = LD/N and the modulation instability (MI) length LMI = (2LDLNL) 0.5 are used to describe the lengths at which the corresponding effect will take place [START_REF]Nonlinear Fiber Optics[END_REF]. Table 1 shows the calculated feature lengths for the ts = 0 and 0.1 μm waveguides. Due to the lower dispersion and nonlinearity, the soliton order N of the ts = 0.1 μm waveguide is 9.4, which is slightly higher than 8.1 for the ts = 0 μm waveguide. For the same input pulses, they both have a very short LNL (0.25 and 0.29 cm) and a long LD (16.3 and 25.5 cm) in ts = 0 and 0.1 μm waveguides, respectively. It indicates that SPM is predominant during the onset of propagation, which can be observed in the range of 0 to ∼1.5 cm in Figures 6(a,b). In the case where the positive and negative chirp compensates each other, the pulse is continuously compressed. Their shorter fission length (2.02 and 2.71 cm) compared to the MI length (2.85 and 3.83 cm) indicates that soliton fission plays a more important role in the broadening process compared to MI in both waveguides. Solitons are deterministically ejected before MI sidebands are fully developed, without the formation of random localized structures in the temporal domain [51]. This avoids the observation of typical SC instabilities associated with broadband noise amplification through spontaneous MI processes [52]. The longer Lfiss of 2.71 cm for the ts = 0.1 μm waveguide compared to 2.02 cm for the ts = 0 μm waveguide means that a longer propagation length is required to trigger soliton ejection, which is observed in Figures 6(a,b).

After soliton fission occurs at 1.6 and 2.2 cm (for the ts = 0 and 0.1 μm waveguides, respectively), some of the ejected solitons have ultrahigh peak power and ultrashort width, which makes LD comparable with (or smaller than) LNL, so that dispersion starts to dominate and causes the individual pulses to break up and spread further apart. Meanwhile, the solitons shed a part of their energy to form DWs between 0.1 and 0.8 ps. The DW generation is stronger in the ts = 0.1 μm waveguide, which is due to the smaller distance between the pump and phase matching point. The temporal overlap of solitons with DW results in their mutual interaction through cross-phase modulation (XPM) and four wave mixing (FWM) [START_REF]Nonlinear Fiber Optics[END_REF]. Such XPM and FWM yield the generation of various peaks merging together to form a broadband supercontinuum. For the case considering normal GVD region pumping, displayed in Figure 6(c), the dispersion and nonlinearity are not mutually compensated. The pulse is progressively deformed, until eventually reaching optical wave breaking (pulse steepening and shock wave formation in Figure 6(f) after ∼5 cm of propagation [53]). In the meantime, SPM becomes progressively weaker and the corresponding SC spectral coverage is significantly limited. (d-f) show the simulated input (blue) and output at 3.1 cm (red) temporal waveforms of the waveguides corresponds to the spectrum in Figures 5(a-c). For the ts = 0 and 0.1 μm waveguides in Figures 6(d,e), soliton fission happens within the 3.1 cm length of the waveguide and, at this distance, the generated solitons features a peak power of 8.3 and 4.65 times higher than the input pulse. All these solitons are much shorter than the input pulse while the shortest one being narrower by a factor of 2N -1 [53]. This results in the shortest soliton after fission is ∼10 fs. When pumping in the normal GVD region, as seen Figure 6(c), the pulse loses its Gaussian shape without chirp compensation, the corresponding peak power being reduced by ∼25.5%, and will continue to collapse during subsequent propagation. To further verify the phase and intensity stability of the SC spectrum generated from the ts = 0.1 μm waveguide, we performed Monte-Carlo simulations including random vacuum noise (one photon per spectral mode with random phase [12]) as well as intensity noise with a relative amplitude set at 10 -4 of the pump pulse. Figure 7(a) shows 256 of those individual simulations with different random noise in gray and their average power in black. From those simulations, we retrieved both the SC spectral coherence and SNR. The first order coherence is calculated with the method in [54] while the SNR is evaluated as SNR(ω) = μ(ω)/σ(ω), where μ(ω) is the average power and σ(ω) is the standard deviation. The calculated coherence and SNR are shown in Figures 7(b,c), respectively, they indicate the excellent stability and coherence degree across the whole SC bandwidth. Here, the coherence remains above 0.8 from 817 nm to 2183 nm, except for some wavelengths with minimal energy. The overall coherence is 0.92 which is comparable with the SC generated from silicon [55] and silicon germanium [56] waveguides.

From the SNR perspective, the highest SNR of ∼17 dB is found at the spectral location of DW formation, and otherwise fluctuates between -2.5 dB and 17 dB across the SC spectrum while the highest SNR is comparable with the SC generated from chalcogenide PCF [57].

Conclusion

In summary, we have demonstrated and studied the dynamics of broadband SC generation in integrated CMOS-compatible and dispersion engineered HDSG waveguides. We have shown that the addition of a lower index oxide slot in the middle of the core of the conventional HDSG channel waveguide can be used to engineer the waveguide dispersion, where one can control the ZDW, phase matching points and the flatness of the dispersion curve. By introducing a 0.1 μm oxide slot, a low and flat anomalous dispersion region can be achieved around the C-band, which leads to SC generation with a spectrum spanning over 1.5 octave, from 817 nm to 2183 nm, when pumped at 1560 nm, with high degree of coherence. We verified our experimental results through numerical simulations based on a GNLSE model, that showed a good match between DW locations and specific spectral features. With the ability to engineering its dispersion, the proposed HDSG slot waveguides thus offers a very promising realization of a fully integrated SC source that can find applications in frequency metrology and frequency synthesis, but also pave the way towards the integration of advanced on-chip functionalities within waveguide architectures with minimal dispersion, compatible with broadband and potentially reconfigurable optical signal processing.

Figure 1 :

 1 Figure 1: (a) Schematic of the slot waveguide structure. (b) SEM image of the fabricated waveguide. (c) Refractive index for thermal oxide (black), PECVD oxide (blue) and HDSG (red), (d) Normalized TE and (e) Normalized TM mode power distributions.

Figure 2 :

 2 Figure 2: Dispersion profiles with fixed waveguide width w = 1.8 μm and changing slot thickness t s = 0.0, 0.1 and 0.2 μm, for the (a) TE and (b) TM modes. Dispersion profiles with fixed slot thickness of t s = 0.1 μm and changing waveguide width w = 1.4, 1.8, and 2.2 μm for the (c) TE and (d) TM modes.

Figure 3 :

 3 Figure 3: Comparison between the TE dispersion profiles of the optimized slot waveguide and (a) the channel waveguide with t = 1.5 μm and changing width w = 1.4, 1.8 and 2.2 μm; (b) the channel waveguide with w = 1.8 μm and changing thickness t = 1.3, 1.4 and 1.5 μm.

  ) and (b) are the comparisons of the TE dispersion profile of the optimized slot waveguide design, with w = 1.8 μm, t = 1.5 μm

Figure 4 :

 4 Figure 4: Comparison of measured SC generation in t s = 0, 0.1 and 0.2 μm waveguides: (a-c) SC evolution vs input pulse energy; (d-f) output spectrum at E p = 621 pJ, corresponds to pulse width of 180 fs and peak power of 2835 W propagating through 3.1 cm long waveguides.

Figure 5 :

 5 Figure 5: Top: computed integrated dispersion. Bottom: comparison between the experimental (blue) and simulated (red) output spectra of the t s = 0, 0.1 and 0.2 μm waveguides.

Figure 6 :

 6 Figure 6: Temporal evolution (a-c) of optical pulses in 6 cm long HDSG waveguides and corresponding temporal profiles at 3.1 cm (d-f). Numerical simulations considering waveguides with t s = 0 (a,d), 0.1 (b,e) and 0.2 μm (c,f).

Figures 6

 6 Figures 6(d-f) show the simulated input (blue) and output at 3.1 cm (red) temporal waveforms

Figure 7 :

 7 Figure 7: Numerical assessment of SC stability: (a) Spectra from 256 Monte-Carlo simulations with random noise (gray) and corresponding average spectrum (black); (b) Calculated first-order coherence degree; (c) Calculated SNR.

Table 1 :

 1 Feature lengths for t s = 0 and 0.1 μm waveguides.

Acknowledgements

This work was funded by the National Natural Science Foundation of China (62105291). This work was funded by the Shenzhen Fundamental Research Program (GXWD20201231165807007-20200827130534001). This work was funded by the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme under grant agreement No. 950618 (STREAMLINE project). Dr. Benjamin Wetzel further acknowledges the support of Région Nouvelle-Aquitaine (SCIR & SPINAL projects).


Michele Alessandro Bucci

Onofrio Semeraro
email: onofrio.semeraro@universite-paris-saclay.fr

Alexandre Allauzen
email: alexandre.allauzen@dauphine.psl.eu

Sergio Chibbaro
email: sergio.chibbaro@universite-paris-saclay.fr

Lionel Mathelin
email: lionel.mathelin@lisn.upsaclay.fr

Preprint

CURRICULUM LEARNING FOR DATA-DRIVEN MODELING OF DYNAMICAL SYSTEMS A PREPRINT

The reliable prediction of the temporal behavior of complex systems is key in numerous scientific fields. This strong interest is however hindered by modeling issues: often, the governing equations describing the physics of the system under consideration are not accessible or, when known, their solution might require a computational time incompatible with the prediction time constraints. Not surprisingly, approximating complex systems in a generic functional format and informing it ex-nihilo from available observations has become common practice in the age of machine learning, as illustrated by the numerous successful examples based on deep neural networks. However, generalizability of the models, margins of guarantee and the impact of data are often overlooked or examined mainly by relying on prior knowledge of the physics. We tackle these issues from a different viewpoint, by adopting a curriculum learning strategy. In curriculum learning, the dataset is structured such that the training process starts from simple samples towards more complex ones in order to favor convergence and generalization. The concept has been developed and successfully applied in robotics and control of systems. Here, we apply this concept for the learning of complex dynamical systems in a systematic way. First, leveraging insights from the ergodic theory, we assess the amount of data sufficient for a-priori guaranteeing a faithful model of the physical system and thoroughly investigate the impact of the training set and its structure on the quality of long-term predictions. Based on that, we consider entropy as a metric of complexity of the dataset; we show how an informed design of the training set based on the analysis of the entropy significantly improves the resulting models in terms of generalizability, and provide insights on the amount and the choice of data required for an effective data-driven modeling.

Introduction

In the present era of mathematization of Nature, the need for modeling is ubiquitous. Focusing on engineering as an example, models are invaluable for predicting quantities of interest, designing and optimizing physical systems, understanding their behavior in an environment, allowing to derive control strategies. Physical models have historically been developed from a combination of observational data and first principles, depending on the amount and type of prior knowledge. This approach is still followed today where the mathematical structure of a model is typically provided and observational data are used to inform it, i.e. learn the coefficients associated with the different terms of the retained structure. Countless variants exist but they generally follow this philosophy. Nonetheless, in many relevant cases, the physics of the phenomena is very complex and the interplay of observation with physical intuition is not sufficient to sort out a formalized model, let alone more rigorous general laws. In these cases, it is tempting to resort to a pure data-driven approach, in which some model is built up entirely from data. The data-driven approach received some attention after the work of [START_REF] Takens | Detecting strange attractors in turbulence[END_REF], [START_REF] Eckmann | Ergodic theory of chaos and strange attractors[END_REF], [START_REF] Sauer | Embedology[END_REF], that provided the theoretical backbone to predictive models based on time-series analysis of low-dimensional systems [START_REF] Kantz | Nonlinear time series analysis[END_REF]. These approaches found new motivation following the significant progresses in the field of Machine Learning (ML), fostered by the massive augmentation of computer-based activities [START_REF] Goodfellow | Deep learning[END_REF]. Even if a clear framework is still missing, the use of such algorithms to build models from data has been introduced in the last decades [see e.g. Voss et al., 1998, Schmidt and[START_REF] Schmidt | Distilling free-form natural laws from experimental data[END_REF] and, thanks to the recent blooming of data-driven techniques, has been considered in the modeling a variety of different physical phenomena [see e.g. [START_REF] Quade | Sparse identification of nonlinear dynamics for rapid model recovery[END_REF][START_REF] Maarten De Hoop | Unsupervised learning for identification of seismic signals[END_REF][START_REF] Steven L Brunton | Machine learning for fluid mechanics[END_REF]. However, it is important to stress that the introduction of deep learning approaches has not necessarily changed the global picture of modeling. The potential versatility of highly expressive neural networks is often associated with a large amount of necessary observational data to suitably inform the associated coefficients of the model. Yet, not every situation enjoys a massive amount of observational data and can afford to inform a widely generic model. In practice, it is more likely that data are scarce due the cost and/or technical difficulty in obtaining them and complex high-dimensional systems may need a quantity of data significantly exceeding the observing capability, [START_REF] Boffetta | Predictability: a way to characterize complexity[END_REF]. In that sense, the lack of data for inferring physical models might have direct consequences on the relevance of these highly expressive models. A strategy to alleviate this limitation is to simplify the postulated structure of the model so that it is less data hungry, with simpler models here understood as involving fewer parameters to be informed from the available data, hence a weaker constraint on the amount of training samples. Another strategy is to leverage prior knowledge one has on the physical system under consideration to guide the structure and the learning process, hence saving on data. For instance, known symmetries can be enforced by design in the structure or in the representation format of the input data, effectively reducing the need for data. A similar idea applies to invariances one may know the system obeys (conserved quantities such as energy, translational or rotational invariances, etc.) or properties such as stability in a given sense. This physics-aware machine learning approach has now become popular and is the motivation of numerous recent efforts [see [START_REF] Raissi | Physics-informed neural networks: A deep learning framework for solving forward and inverse problems involving nonlinear partial differential equations[END_REF][START_REF] Zhu | Physics-constrained deep learning for high-dimensional surrogate modeling and uncertainty quantification without labeled data[END_REF][START_REF] Kashinath | Physics-informed machine learning: case studies for weather and climate modelling[END_REF], von Rueden et al., 2021, among many]. In that sense, the modeling of chaotic dynamical systems is a natural play-field [START_REF] Pathak | Using machine learning to replicate chaotic attractors and calculate lyapunov exponents from data[END_REF], [START_REF] Borra | Effective models and predictability of chaotic multiscale systems via machine learning[END_REF]. Remarkably, previous works seem to indicate an intriguing capability to give accurate forecasting of chaotic systems, even relatively high-dimensional [START_REF] Pathak | Model-free prediction of large spatiotemporally chaotic systems from data: A reservoir computing approach[END_REF], apparently improving over deterministic embedding techniques [START_REF] Kantz | Nonlinear time series analysis[END_REF]. These results therefore deserve to be thoroughly analyzed. In the present work, we aim to shed some light on the challenges in data-driven learning dynamical systems and to explore the idea of leveraging prior information, here in the form of known organization of the dynamics, to inform the process of collecting relevant observations. Specifically, we adopt a curriculum learning strategy for dynamical systems from first principles, and assess it with respect to common sampling strategies. Curriculum learning relies on the idea that machines, like humans, learn best if first exposed to "easy" situations prior to dealing with "more complex" ones when well into the learning process. A possible rationale is that curriculum learning can be interpreted as an adaptive regularization method for optimization of a non-convex objective function. By progressively introducing more complex samples, the associated cost landscape evolves from near-convex to a more complex function exhibiting many more local extrema. Sorting training examples in this way may hence favor the optimization process toward a better extremum. This concept has been mostly developed in robotics and control of systems, [START_REF] Sanger | Neural network learning control of robot manipulators using gradually increasing task difficulty[END_REF], and was formalized more recently in the context of machine learning, [START_REF] Bengio | Curriculum learning[END_REF]. While it still lacks a firm theoretical ground despite some efforts [START_REF] Weinshall | Curriculum Learning by Transfer Learning: Theory and Experiments with Deep Networks[END_REF], curriculum learning has been widely illustrated and successfully employed in a variety of applications, ranging from robotics to video games, from computer vision to natural language processing and reinforcement learning, [START_REF] Narvekar | Curriculum learning for reinforcement learning domains: A framework and survey[END_REF], [START_REF] Soviany | Curriculum learning: A survey[END_REF]. What makes a situation, or training sample in the context of machine learning, "easy" or "complex" is case-specific. For instance, in classification problems, the complexity of a particular example presented to the machine may be associated with the distance to the decision surface (e.g., SVM) or with the amount of noise in the data. Here, we illustrate the impact of curriculum learning based sampling strategies, by considering the celebrated Lorenz'63 system. While deterministic and low dimensional, the Lorenz system features many challenges encountered in more complex configurations, such as sensitivity to initial conditions and an attractor dimension d A A A significantly lower than the ambient dimension d, here d A A A 2.06 < d = 3. It is thus somehow a minimally complex unit and justifies its relevance for illustrating some challenges in learning dynamical systems. We associate the complexity of a training element with its information entropy, in the sense of Shannon, and recognize that the solution of dynamical systems in the vicinity of an unstable fixed point is well described by a linear model and exhibits a low entropy. This observation is key and allows a principled sampling strategy in collecting observations of increasing complexity from the system at hand as it evolves from an unstable fixed point. We explore and compare different strategies for sampling a dynamical system in the aim of learning its behavior and predicting its state over time. These strategies only differ by the region in the phase space where observational samples are taken from, every other aspects being identical (sampling frequency, number of samples, observed quantity, etc.). These different configurations are chosen to be typical of situations ranging from no prior information, besides an estimated upper bound of the dimension of the dynamical system up to approximate knowledge of some invariants of the system. We carry out a comprehensive statistical analysis in order to systematically test the impact that different sampling strategies have on the quality of the resulting models in terms of prediction abilities and generalizability. From an algorithmic point of view, recurrent structures such as auto-regressive models or recurrent neural networks are popular choices in such situations. In the present study, we focus on the LSTM (Long-Short Term Memory) architecture [START_REF] Hochreiter | Long short-term memory[END_REF] for its widespread use, including applications in prediction of dynamical systems behavior Gers et al. [2002]. Moreover, it has been recently shown that the performance given by different architectures including reservoir computing is similar [START_REF] Pantelis R Vlachas | Backpropagation algorithms and reservoir computing in recurrent neural networks for the forecasting of complex spatiotemporal dynamics[END_REF], although for low dimensional systems like Lorenz'63, LSTMs appear to give the best results. Given this choice, we also provide best-practice guidance concerning the possible impact of initial conditions and the memory effects on the results of supervised learning models, beside analyzing the possibility of using less data while preserving generalizability and quality based on curriculum learning strategies. The paper is organized as follows. General principles in learning a system from data are briefly discussed in Sec. 2 to set up the stage and define notations. The general class of models we consider in this work is presented in Sec. 3. The specific configuration we retain for the present discussion is the Lorenz'63 introduced in Sec. 4, together with the different sampling strategies. The results are gathered and presented in Sec. 5 before the discussion in Sec. 6 closes the paper.

Theoretical background 2.1 Learning dynamical systems

The present work lies in the framework of learning a model to approximate the time evolution of a state of the system at hand from observations s s s(t) ∈ R R R n s . To do so, one needs a representation for the state and a dynamical model for describing its evolution in time. A representation of the state is given by a function g mapping observations s s s Γ (t) over a measurement time horizon t ∈ Γ, to a set of coordinates v v v ∈ R R R d v . Several mappings have been proposed in the literature, including delays by time intervals ∆t, v v v i = s s s(t -i∆t) for i = 0, 1, . . ., and time derivatives of varying orders i, v v v i = s s s (i) , [START_REF] James | Equation of motion from a data series[END_REF], [START_REF] Gibson | An analytic approach to practical state space reconstruction[END_REF]. To obtain the size of the coordinate space, we assume that the state of the system u u u lies on an attractor A A A of dimension d A A A < d that is a boundaryless bounded, smooth sub-manifold of R R R d . In this case, the Whitney embedding theorem [START_REF] Whitney | Differentiable manifolds[END_REF] shows that if the number of coordinates d v is such that d v > 2 d A A A then the dynamics of the state can be entirely captured by the new coordinate system. The Takens embedding theorem [START_REF] Takens | Detecting strange attractors in turbulence[END_REF], [START_REF] Sauer | Embedology[END_REF] provides a way to make this result practical by asserting that one can take this number of time delays. Recent results have extended this theorem by providing guaranteed recovery properties based on the measurement operator, type of attractor, and other parameters of the problem, [START_REF] Eftekhari | Stabilizing embedology: Geometrypreserving delay-coordinate maps[END_REF]. In this work, we assume no prior knowledge of the governing equations of the model and rely on a purely datadriven approach. The dynamical model for describing the time evolution of v v v(t) ∈ R R R d v can be formulated in various ways, including feed-forward neural networks, owing to their high expressiveness. As an alternative, for instance in case an upper bound of the attractor dimension is unknown due to process noise or measurement noise, a recurrent model F θ θ θ can be employed for describing the dynamics of v v v. The model is parameterized by

θ θ θ ∈ R R R n θ θ θ and writes v v v(t) = F θ θ θ (v v v) , (1)
where a dotted quantity denotes its time-derivative. Learning a model for the system under consideration takes the form of a classical supervised learning problem based on the misfit between observations s s s Γ ≡ s s s (u u u (t ∈ Γ)) collected through a measurement operator h : u u u → s s s over a time domain Γ = [t 0 ,t T ] and predictions, in the sense of some norm • Γ . It typically results in an optimization problem in terms of the model parameters, possibly including a regularization term r:

θ θ θ ∈ arg inf θ θ θ ∈R R R n θ θ θ s s s(t) -h v v v t; θ θ θ 2 Γ + r θ θ θ , s.t. v v v t; θ θ θ = F θ θ θ (v v v) .
(2)

Figure 1: Sketch of the retained architecture. Each hidden layer is made of LSTM cells, while the last layer is fully connected.

Collecting data

The main issue impacting the predictions based on observations is the quantity of data actually needed to approximate the dynamics in a given sense. The issue is related to the ergodic theory of dynamical systems, whose founding idea is that the long-time statistical properties of a deterministic system can be equivalently described in terms of the invariant (time-independent) probability, µ, such that µ(S) is the probability of finding the system in any specified region S of its phase space. If the trajectories of a d-dimensional ergodic system evolve in a bounded phase space Ω ⊂ R R R d , the Poincaré recurrence theorem Poincaré [1899] ensures that analogues exist as it proves that the trajectories exiting from a generic set S ∈ Ω will return to such set S infinitely many times. The theorem holds for any points in S, almost surely. It was originally formulated for Hamiltonian systems, but it can be straightforwardly extended to dissipative ergodic systems provided initial conditions are chosen on the attractor. In d-dimensional dissipative systems, the attractor

A A A typically has a dimension d A A A < d. The probability µ B d
x (ε) of finding points on the attractor which are in the d-dimensional ball of radius ε ( with resolution ε

1) around x ∈ A A A ⊂ R R R d : µ B d x (ε) = B d x (ε) dµ( x) ∼ ε d A A A . (3)
When d A A A is not an integer, the attractor and its probability measure are fractal. The attractor is typically multifractal [START_REF] Paladin | Anomalous scaling laws in multifractal objects[END_REF], but we can assume without lack of generality that it is homogeneous. Given that the Poincaré theorem proves the existence of good analogues, the issue is shifted to finding out how much time is needed for a Poincaré cycle to end. This key problem was considered by Smolouchoski and resolved by [START_REF] Kac | Probability and related topics in physical sciences[END_REF], who showed that for an ergodic system, given a set S ≡ B d x (ε), the mean recurrence time τ S (x ) relative to S is inversely proportional to the measure of the set S, i.e.,

τ S (x ) ∝ 1 µ(S) ∼ ε -d A A A , (4)
where the average • is computed over all the points x ∈ S according to the invariant measure. Therefore, if we require ε to be small and if d A A A is large, the average recurrence time becomes huge, a symptom of the curse of dimensionality in data-driven methods.

3 Machine-learning approach

Architecture

As discussed above, a set of d-dimensional observations retained to express a state of the system under consideration might not be a state vector, in the sense that it does not necessarily contain enough information to uniquely define its future evolution as an autonomous deterministic dynamical system. Even though the system were Markovian, it might not be so in the retained representation and may implicitly involve latent variables. This is typically the case in situations of partial observability or when measurements are affected by noise. As discussed in Sec. 2, under conditions rigorously studied in the embedology line of research, [START_REF] Eckmann | Ergodic theory of chaos and strange attractors[END_REF], [START_REF] Sauer | Embedology[END_REF], latent state variables can be substituted with observations over a recent past, that is, temporal information is traded for unavailable high-dimensional phase space information. In the present work, we consider the LSTM (Long-Short Term Memory) architecture. Compared with standard multi-layer perceptrons, LSTMs additionally feature internal variables ζ ζ ζ carried-over through the recurrence and effectively acting as a "memory". Training an LSTM then allows to derive a discrete-in-time model predicting the future "state" v v v k+1 ≡ v v v (t k + ∆t) from the current state and memory:

v v v k+1 = F θ θ θ (v v v k ; ζ ζ ζ k ) . (5)
The retained architecture comprises 2 LSTM layers of 50 neurons each and one fully connected layer, with tanh activation functions in the LSTM cells and sigmoids for the state output (see the sketch in Fig. 1). The present architecture was deemed as a good compromise between accuracy and computational expenses, after an analysis with comparable architectures with one and three layers has been performed. A 1-step ahead model is learned based on a whole trajectory, or a set of trajectories, in the L 2 -sense. The training is performed using the ADAM optimizer and an adaptive learning rate from 0.01 to 10 -5 , with a decay factor of 0.07. The maximum number of epochs is limited to 1000.

Strategies for sampling

Learning a model here consists in estimating the best set of parameters θ θ θ in the sense of the optimization problem defined in Eq. 2. In the context of neural networks, this acts as the loss function used for training. The regularization term may promote several properties of the optimal solution θ θ θ , such as the popular choice of low magnitude (ridge regression), r (θ θ θ ) ∼ θ θ θ 2 , or sparsity, r (θ θ θ ) ∼ θ θ θ 1 , possibly mimicked with a suitable dropout strategy. We adopt a discrete-in-time viewpoint, in a formulation consistent with the discrete model in Eq. 5. The loss function L associated with the training problem then writes

L (θ θ θ ) = K ∑ k=1 s s s k -h (v v v k (θ θ θ )) 2 2 + r (θ θ θ ) , (6)
with K such that T = K ∆t and the 1-step ahead dynamical model

v v v k+1 = F θ θ θ (v v v k ; ζ ζ ζ k ).
The loss is here simply defined in terms of the Euclidean distance but many alternative definitions could be employed. For instance, optimal transportbased metrics, such as Dynamic Time Warping (DTW), can prove useful when observations come at a poorly known or controlled time pace. While the way the prediction from the model is compared with observational data is pivotal, the focus in this work is instead on the relevance of the training data. The present findings are believed to hold, disregarding the definition of the loss. A standard L 2 distance, implicitly assuming regularly sampled measurements, is hence retained. The measurement operator is the identity, h ≡ I, so that the measurements are associated with the state of the system, s s s k ≡ u u u k , ∀ k. No noise affects the observations, so that one is in a full observability situation. In this context, the observed system is Markovian in the representation space of the observations. However, a recurrent architecture is here considered as a general structure for learning, relevant for a wide class of situations. It further has the benefit of improving upon the robustness of the learning by virtue of some degree of redundant representation. Without loss of generality, we disregard the regularization term r (θ θ θ ) in the definition of the loss function which finally writes

L (θ θ θ ) = K ∑ k=1 s s s k -v v v k (θ θ θ ) 2 2 , s.t. v v v k+1 = F θ θ θ (v v v k ) , v v v 0 = s s s 0 . (7)
While this context is favorable for learning (low-dimensional deterministic system, noiseless observations, full observability, etc.), it will be seen below that the measurement strategy still significantly affects the learning performance.

Curriculum learning, metrics and complexity of the datasets

As mentioned in the Introduction, we aim at verifying the extent to which a measurement strategy inspired by the curriculum learning principle can be applied to data-driven learning of dynamical systems and how this impacts the generalizability and quality of the resulting models. The rationale is the adaptive regularization of the objective function when moving from "simpler" cost landscapes associated with simpler data towards more complex ones. From this perspective, any indicator of complexity can be used; here, we consider entropy as a measure and compute it based on singular value decomposition (SVD) [START_REF] Varshavsky | Novel unsupervised feature filtering of biological data[END_REF]. Given an embedding for the measurements s s s based on time delaying

s s s = [s(k), s(k + τ), . . . , s(k + (n e -1)τ)] , (8)
with k an integer and τ the time delays, the matrix S S S is defined as

S S S = [s s s(0), s s s(τ), s s s(2τ), . . . , s s s(N -(n e -1)τ)] , (9)
with n e indicating the number of τ. The SVD is performed on the matrix and M singular values σ SVD i i are obtained such that the entropy H is defined as

H SV D = - M ∑ i=1 σ SVD i log 2 (σ SVD i ), (10)
with

σ SVD i = σ SVD i / ∑ M j=1 σ SVD j .
The SVD entropy provides an estimation of the number of eigenvectors that are needed for an adequate representation of the dataset. In that sense, it measures the dimensionality of the data. With this metric, we will assess in the next section the complexity characterizing each of the training datasets that will be compared.

The second metric introduced in the article for the evaluation of the models is the correlation dimension, computed using the algorithm by [START_REF] Grassberger | Characterization of strange attractors[END_REF]. Considering a ball of radius ε, the probability of finding a point scales as C(ε) ∼ ε d 2 , where d 2 is an indication of the dimensionality. Using the Euclidian norm, the probability C (ε) can be efficiently computed as the correlation sum

C (ε) = lim N→∞ 1 N (N -1) N ∑ i, j=1 i = j Φ (ε -s i -s j ) , (11)
with Φ the Heaviside step function, s s s the components of the embedded variable ranging i, j = 1, 2 . . . N, with a fixed time increment; N is the number of pairs. Due to the divergence of trajectories, most pairs (s i , s j ) with i = j are usually uncorrelated. Here the correlation dimension d 2 is considered for the information provided on the fractal dimension of the attractor as well as the connections with the intrinsic dimension, a metric often found in machine learning [START_REF] Camastra | Intrinsic dimension estimation: Advances and open problems[END_REF]. Alternatively, the Lyapunov exponents can also be adopted as criteria.

4 Learning the Lorenz'63 system

Lorenz'63 system

To illustrate the different learning strategies, we consider the celebrated Lorenz'63 system Lorenz [1963], originally introduced to mimic the thermal convection in the atmosphere. It consists of three coupled nonlinear ordinary differential equations (ODEs):

   ẋ = σ (y -x) , ẏ = ρx -y -xz, ż = xy -β z, with {σ , ρ, β } ∈ R R R + .
The system of ODEs is integrated in time with a 4 th order Runge-Kutta scheme and a time step ∆t = 0.01. Constants σ , ρ, β are positive parameters, related to dimensionless scales such as the Prandtl and the Rayleigh numbers. For ρ > σ (σ + β + 3) / (σβ -1) and σ > (β + 1), the solution u u u := (x, y, z) exhibits a chaotic behavior. It is important for the following to note that the system features three unstable fixed points u u u FP , associated with a vanishing time-derivative:

u u u FP + = β (ρ -1), β (ρ -1), ρ -1 , u u u FP -= -β (ρ -1), -β (ρ -1), ρ -1 ,
u u u FP 0 = (0, 0, 0) . In the following we shall consider the parameters fixed to β = 8/3, σ = 10, ρ = 28 as in the original work, such that the solution is in the chaotic regime. In this regime, the Kaplan-Yorke dimension of the attractor is d A A A 2.06.

Set-up of the learning method

We turn our attention to the training set T = {s s s k } K k=1 for learning the Lorenz'63 system. It is a crucial ingredient of the learning process and largely conditions the quality of the resulting model. We consider four different training sets (T ergo , T split , T ens , and T FP ) of identical settings (sampling frequency, number of samples, measurement operator, etc.) but following different sampling strategies.

Ergodicity-compliant training sets, T ergo and T split

The T ergo dataset is a long trajectory on the attractor. We consider it the reference case as it provides sufficient statistics for learning a reliable and accurate model. With ε = 10 -2 and the Kaplan-Yorke dimension d A A A 2.06 of the Lo-renz'63 for the retained settings {σ , ρ, β }, [START_REF] Kuznetsov | The Lorenz system: hidden boundary of practical stability and the Lyapunov dimension[END_REF], the required number of samples resulting from a constant sampling period of ∆t = 0.01 is about 27, 000 strobes in time. A second dataset, T split , is derived from T ergo : it consists of one, long trajectory on the attractor, yet split in 9 chunks of equal size and reshuffled, such that each chunk is thus composed by 3000 samples. As it will be clear when discussing the results (Sec. 5), this dataset will allow us to investigate some effects of memory. The training set T FP exploits prior knowledge on the structure of the system at hand, namely the fixed points. As mentioned in the Introduction, it is well known dynamical systems are organized around critical features, such as fixed points, and that invariant sets provide information on the behavior of the solution, [START_REF] Gilmore | Connecting curves for dynamical systems[END_REF]. Besides fixed points, these features include higher-dimensional invariant sets such as slow invariant manifolds. The sampling strategy here consists in acquiring observations from the system originating from its (unstable) fixed points, in contrast with other training sets obtained from the attractor. Specifically, three 3000-sample trajectories are acquired from each of the three fixed points of the Lorenz'63 system. These trajectories originate from a small deviation from each fixed points along each of the d = 3 directions given by the eigenvectors of the local Jacobian matrix. Note that if the Jacobian is not available, a set of orthonormalized deviations can also be used. The resulting training set then consists of 3 × 3 × 3000 = 27, 000 samples, consistent in size with the other training sets illustrated in Fig. 2.

Random sampling, T ens

Finally, we consider the sampling of short trajectories on the attractor. Nine trajectories are considered, each originating from a state randomly chosen on the attractor. Each trajectory is 3000 sample long as for the other strategies introduced above so that the total size of the training set is again 9 × 3000 = 27, 000 samples. This training set is representative of the widespread situations where information about a system come in several independent glimpses, as resulting from several measurement campaigns. It is illustrated in Fig. 2.

Entropy characterization

Using the definition of SVD entropy provided Eq. ( 10) with order n e = 3 and time-delay τ = 1, we characterize in the following each of the training sets. The SVD-entropy algorithm has been chosen after a comparison among different strategies for the entropy estimation. First, we considered convergence of the estimates over long trajectories emanating from different initial conditions, in order to evaluate the consistency of the results in ergodic datasets; second, a validation based on entropy estimators requiring a reduced number of data, namely the Approximate Entropy algorithm by [START_REF] Steven | Approximate entropy as a measure of system complexity[END_REF], has been performed. Good agreement was found between the different entropy estimators.

From the physical viewpoint, the entropy for a time series is defined on the probability transition between the embedded states (unfolded through SVD). In Fig. 3, the entropy distributions are compared for each of the training sets: quite neatly, the dataset T FP is associated with the lowest entropy while the remaining datasets are all characterized by the same levels. From the physical viewpoint, this behavior is explainable by considering the higher predictability of the dynamics in the vicinity of the fixed points, where the behavior of the states close to them is well approximated by the linear limit. High predictability corresponds to less required information. As we will show in the following, this leads to an apparent contradiction: despite less information is contained in the structured datasets -in a Shannon sensethe possibility of training LSTM with well structured datasets characterized by increasing entropy in time allows to deploy a curriculum learning strategy. when the same model is evaluated with an unknown initial condition (test set), the performance deteriorates and the accuracy is lost after less than 4 time units. When a unique ergodic-compliant trajectory T ergo is used to train a model (bottom row), the prediction from an unknown initial condition improves significantly.

Results

We now assess the quality of the models learned from various observation datasets of the Lorenz'63 system.

How much training information?

The influence of the amount of information to train a model is first considered. Two models with the same architecture (cf. Sec. 3) are trained from two different datasets. These datasets consist of observations of the solution of the Lorenz system along the same trajectory on the attractor, but of different length. One dataset, T ergo , comprises 27,000 samples while the second one T short is its restriction to the first 3000 samples, so that T short ⊂ T ergo . Once trained, the resulting models, hereafter respectively termed F ergo and F short , are evaluated in terms of their long-term prediction capability.

Results are gathered on Fig. 4. The performance of the models is here evaluated when predicting from a known (i.e., contained in their training set) and unknown (not contained) initial condition, to assess the generalization capability of the models. For all cases, memory is initialized to zero.

A known initial condition is first considered. As can be seen from the top row (left column) where it is plotted against the truth (red), the F short model allows a reasonable quality of the prediction (plotted in blue) over a significant time horizon of about 6 Lyapunov times (7 time units). Since training the models is achieved via a stochastic procedure, e.g. stochastic gradient descent and random initialization of the parameters, and the landscape of the loss function is not globally convex, the underlying optimization problem ends up in a local optimum and bears an aleatoric contribution.

To strengthen the above conclusions to account for the statistical uncertainty of neural network modeling, we hence consider an ensemble of 100 identical independent models and evaluate their performance [START_REF] Agarwal | Deep reinforcement learning at the edge of the statistical precipice[END_REF]; the final evaluation is thus based on an aggregate metric that allows to provide a robust assessment of the relative performance across different trainings based on the different datasets. The right column of Fig. 4 shows the mean and standard deviation of the prediction of the ensemble of models for the same conditions as in the left column. From the top row, one may deduce that the prediction of the ensemble of models is very consistent and that the learning is robust. With the same short training set T short , the trained model is now evaluated on an initial condition not contained in the training data, see Fig. 4 middle row. Compared with the previous situation (top row), the performance is seen to deteriorate, both for a single model (left column) and in average (right). In particular, the time horizon before the prediction significantly differs from the truth is much shorter, indicating a poor generalizability properties. When more training data are used, T ergo , the resulting model remains accurate over a long time horizon even when predicting from an unknown initial condition, as seen from Fig. 4 (bottom row). Again, this is a reliable finding as evidenced from the ensemble performance seen on the right column.

In Fig. 5 a complementary analysis is shown, considering the average cumulative 2-norm error over a short time period of 1 time-unit. The averaging is performed considering 100 models for each of the 4 sampling strategies. The models are trained using random initialization of the memory. In the inset (a), the prediction emanating from trajectories included in the training set is shown. It is possible to see that the error is bounded and constant, at least within the short time window we consider here; the ergodic case is characterized by higher error due to a larger number of strobes back-propagated during the training in each batch. In (b), the error is shown for the testing sets. First -coherently to Fig. 4 -the error increases while time increases; secondly, larger deviations are observed for the random dataset, while the cumulative error of the other cases is more comparable. 7a, while the distribution of the models is reported on the vertical axis using the same number of bins. We order from light cyan to darker blue models performing progressively worse in terms of resulting d 2 accuracy; note that the level set of the blue gradient is normalized for each subplot, such that only a qualitative information on the distribution along r is included in the graphs.

Figure 7: Illustration of the model parameters via a t-SNE representation. The LSTM memory is initialized to zero.

0 1 2 Ergodic 8% 0 1 2 Random 20% 0 1 2 Structured 2% 0 1 2 Ergodic split 10%
Figure 8: The figure depicts the empirical distribution of the d 2 dimension associated with the predictions obtained for models trained by initializing the LSTM memory at random, using the four different learning strategies. For each training set, the proportion of models associated with an estimated dimension differing by more than 25% from the true dimension d A A A = 2.06 is indicated.

These results illustrate the fact that the training data, relied on to learn a model, have to be informative enough of the system under consideration to allow for generalizable models to be inferred. In particular, the superior performances shown in the top row of Fig. 4 and the low error shown in Fig. 5(a), when the initial conditions are contained in the learning data-set appear to be due to overfitting.

Relevance of the data

Besides the mere size of the training set, a critical aspect is its relevance for informing a model in view of a given objective function. To shed light on this pivotal aspect, the four different training sets introduced in Sec. 3.2, namely the ergodic (T ergo ), the ergodic-split (T split ), the random-sampled (T ens ) and the structured based on the fixed points (T FP ), all of the same size of 27, 000 samples, and analyzed using the SVD entropy in Sec. 4.2.4, are now evaluated in terms of the quality of the associated resulting model. To avoid any bias due to the initialization of the parameters, an ensemble of 100 models is trained for each situation and the comparison is made in terms of the ensemble performance [START_REF] Agarwal | Deep reinforcement learning at the edge of the statistical precipice[END_REF]. A first metric to quantify the quality of the resulting model is the d 2 -dimension of the predicted model, computed by estimating the correlation sum in Eq. 11. This criterion is not associated with a given initial condition and is affected by the way the trained model consistently reproduces the fractal nature of the Lorenz'63 system. Results are gathered in Fig. 6 in terms of a histogram of the estimated d 2 dimension for the ensemble of 100 models for each of the four training sets. For a more quantitative appreciation, a threshold is defined when the dimension resulting from the trained model differs by more than 25 % from the truth d A A A = 2.06. As expected, the ergodicity-compliant training set T ergo leads to a good performance, with most models associated with a d 2 dimension close to the true value and a few degenerate models such that d 2 ≈ 0. Only 14 % of the models are beyond the threshold. Models trained from T FP are also seen to achieve a similar level of performance. Again, almost all these models are associated with a rather good d 2 dimension, with very few degenerate models. Only 8% of the models are beyond the threshold, an even better performance than models trained from T ergo . In contrast, the T ens (random) and T split (split) dataset leads to a rather poor performance model, with a significant part of the models associated with d 2 ≈ 0 and respectively 34 and 26 % of the models beyond the threshold. An observation is in order in this regard.

The models with d 2 ≈ 0 are all characterized by being stationary: the dynamics of the system stabilizes towards the fixed points after a short transient. On the other hand, only a limited number of models is characterized by d 2 ≈ 1; in this case the temporal dynamics settles on periodic orbits. Interestingly, in most of the stable models, the fixed points are correctly predicted, also for the random datasets where they are not included explicitly; however, their stability features are not correctly predicted, as already stressed.

It may be not quite surprising that models trained with T ens and T split are associated with a similar performance. Since the only difference between the two datasets is that the samples in T ens originate from a collection of 9 independent segments of trajectory, randomly located on the attractor, while T split is made of 9 contiguous segments. Yet, such behavior for the T split dataset appears to be not trivial, because it contains exactly the same information of the T ergo one about the phase space. Further assessment of these observations can be made using a different technical tool, namely t-distributed stochastic neighbor embedding (t-SNE) of the parameters, shown in Fig. 7a. It is essentially a dimensionality reduction technique convenient for assessing high-dimensional quantities [START_REF] Van Der Maaten | Visualizing data using t-sne[END_REF]. The model parameters vectors are here mapped onto a 2-dimensional space where each model can hence be represented in terms of its corresponding coefficients γ 1 and γ 2 . F ergo and F FP are mostly clustered in the center of the points cloud, indicating that, while the models are different, they are however similar to each other in terms of parameter distribution. Instead, the models obtained from the random sampling strategy (T ens ) and the ones obtained stacking chunks from an ergodicity-compliant trajectory (T split ) are distributed in the outer region. LSTM parameters are initialized as independent realizations of a Gaussian random variable, whose t-SNE representation is typically visually close to a Gaussian as well since it penalizes deviation from a t-Student distribution. Recent works have shown that neural networks associated with parameters weakly evolving during training while yet able to significantly improve the loss function, are often associated with good generalization properties, Jacot et al. [2018]. In the present study, the t-SNE representation of models trained from the ergodicity-compliant and the structured datasets both appear close to a Gaussian distribution, a further indication that they might enjoy good generalization properties. In contrast, the t-SNE distribution of the models trained from T ens or T split appear to deviate significantly from their initial distribution, which corroborates well with their poor performance. This distribution indicates that the best models are robust, in the sense of having similar parameters. Furthermore, training using information about fixed points appears equivalent to that based on a long trajectory. This observation is also quantified in terms of a radial distribution, see Fig. 7b. Both the models trained from T ergo and T FP present a distribution with highest density around zero, while T ens and T split lead to a bias in the parameter distribution, hence being associated with poor models. In view of the entropy levels associated with each of the datasets (see Sec. 4.2.4), we conclude that organizing the datasets around the fixed points, and leveraging the simpler, linear dynamics characterizing the vicinity of those points corresponds to deploying a curriculum learning for the training of fully data-driven models. We postulate that a possible rationale behind the superior performance can indeed be understood by a simpler landscape of the cost associated with the data obtained where the dynamics is well approximated by linear behaviour. On the other hand, the curriculum learning assumption does not fully explain the discrepancies arising between models that -in principle -are qualitatively the same (namely, T ergo and T split ). For this reason, we turn our attention to a crucial element: the role of the memory in the recurrent models.

Handling the memory

Memory is considered in this section. Instead of the null memory (initialized to zero) considered so far and classically found in many applications, an initial memory chosen at random, drawn from a multivariate Gaussian distribution N (0 0 0, I) is used. Compared to the null initial memory, the performances are seen to significantly increase for all models, as evidenced in Fig. 8. This is consistent with the fact that a non-zero hidden state improves the numerical conditioning of the training steps and prevents degenerate gradients, [START_REF] Mostafa Mehdipour Ghazi | On the initialization of long short-term memory networks[END_REF]. The contrast in performance between models learned from T ergo and T split , which rely on the same information, leads to identify a key point.

The LSTM model trained with T ergo learns to best fit one long trajectory. In the case of T split , it learns to best fit 9 short trajectories, which however, when stitched together, match the T ergo long trajectory. Yet, their performances are consistently very different in the null initialization situation, emphasizing that the initial memory of the retained recurrent network (LSTM) strongly affects the resulting performance. In the case of T ergo , the impact of an incorrect initial memory is diluted within a long trajectory and weakly affects the learning quality. When short trajectories are considered instead, the relative impact of an incorrect initialization is stronger and can introduce a significant bias, hence a poor resulting model. This observation is further supported by the model learned from T FP . While it also relies on samples originating from 9 distinct, non sequential, segments, its performance is good, as evidenced by Fig. 6. What is unique about this particular training set is that the samples are from segments of trajectories originating from (unstable) fixed points, thus associated with a locally linear dynamics. This feature allows the incorrect initial memory to be progressively wiped out in time without significantly affecting the predicted dynamics.

Can we use less data?

The results shown so far lead to slightly counter-intuitive conclusions: from one hand, low entropy indicates that the dynamics described by the datasets T FP lead to a simpler training; on the other hand, simpler dynamics means high predictability and thus less information. In that sense, one might wonder whether is possible to leverage these considerations for reducing the amount of data used during the training: in fact, the linear regions are less informative as they are redundant, in particular when we consider "slow" trajectories emanating from the fixed points located in the two lobes of the Lorenz'63 attractor. In order to provide more effective datasets, while preserving the curriculum learning idea, we consider here shortened trajectories where the initial samples are removed. The resulting characterization of the shortened datasets, based on the SVD entropy, is contained in Fig. 9: the entropy increases progressively as the dataset contains (relatively) more data stack in the nonlinear region. By considering these shortened datasets, we train LSTM models and assess their robustness using the d 2 metric applied in the previous sections, and study the associated models. As shown in Fig. 10, it is found that we can reduce the number of samples to a total of 9 × 2, 000 = 18, 000, as the behavior is essentially the same as for the original T FP models, i.e. not prone to overfitting and with low probability to produce biased results. When further reducing the data, results start to worsen until they definitively deteriorate for the last case 9 × 1, 000 = 9, 000. Nonetheless, we can conclude that a reduction of 33% of the samples does not lead to loss of performance for the analyzed case. While the focus of the analysis is on data selection and their impact on the quality of performance, it is worth noting that any constraint applied to the cost function can further reduce the amount of data while preserving the performance of the system, for instance by explicitly introducing information about the stability of the system (Lyapunov exponents or the Jacobian matrix) or introducing multi-objective cost functions [START_REF] Quade | Prediction of dynamical systems by symbolic regression[END_REF][START_REF] William | A probabilistic and multi-objective analysis of lexicase selection and ε-lexicase selection[END_REF]. In that sense, a combination of physics constraints, data selection/distillation and architecture optimization can deeply improve the prediction abilities of data-driven models with a significant reduction of the data necessary for training.

Figure 10: Using the same postprocessing adopted in Fig. 7, the radial distribution of the t-SNE projected model parameters is shown for the model based on the datasets obtained starting from T FP , and shortened of the initial states.

We order from light cyan to darker blue models performing progressively worse in terms of resulting d 2 accuracy; note that the level set of the blue gradient is normalized for each subplot. As a reference, we consider the case with 3000 (e). A rather compact clustering of the models, i.e. r < 0.5, is found up to t = 2000 (insets c and d), suggesting good behavior of the trained models. Shorter trajectories lead to less robustness in terms of model weights as well as larger errors when considering d 2 .

(a)

(b) (c) (d) (e)

Concluding discussion

In this work, we apply the concept of curriculum learning for data-driven modeling of a complex dynamical system: this concept is inspired from human behavior and is based on the idea that better learning is achieved when machines are trained using data of increasing complexity. Among the possible explanations for motivating the effectiveness of this learning strategy, a possible one is related to the shape of the cost landscape, that evolves from being near-convex to more complex as the complexity of the data increase. By sorting the data in this way, the optimization process is favored towards better extrema. Motivated by this methodology, we have performed a thorough analysis of some aspects inherent to the process of learning a complex physical system from observational data using the classical Lorenz'63, owing to its simplicity and hence the possibility to carry out an exhaustive campaign of simulations. As a learning method, we have considered LSTM neural networks, an established recurrent structure often used in such a context. We addressed the following issues:

1. Is there a minimum amount of data needed to obtain a robust model able to generalize? 2. Is it possible to go beyond this limitation thanks to some "data-based" knowledge of the system using a metric? 3. What is the impact of the initial state of the retained model on its resulting performance?

The first two questions are rather fundamental. In particular, the first question concerns the amount of data required for training and it is related to what is known for classical embedding approaches, which is bound to fail when time-series data are not enough to reconstruct the entire phase-space. Generally speaking, in absence of prior information on the system at hand besides an estimated upper bound of its dimension, the ergodic theory, and notably the Kac lemma, provides a reliable criterion on the minimal amount of data necessary to train a model. This amount of data grows exponentially with the dimension of the attractor. Models trained on a smaller amount of training data have been shown to consistently lead to poor performance, as illustrated both via the time-series prediction and the resulting model estimated dimension assessment criteria. Our analysis hence suggests that recent results on dynamical systems showing a larger predictability horizon appear more as a result of overfitting than an intrinsic achievement of neural networks. The second question is addressed by comparing different training sets. In particular, a training set acquired from a system initially close to one of its unstable fixed points exhibits a low entropy, low complexity, behavior owing to the essentially linear local regime, before reaching a higher complexity, fully nonlinear, regime when getting closer to the attractor. This essentially allows to disentangle the contribution of the different regimes of the underlying system and to sequentially focus the learning effort on these different regimes from specific parts of the training set, effectively resulting in a curriculum. Interestingly, shortening the trajectories emanating from the fixed points and following insights based on the entropy contents lead to results consistent with the ones obtained with a larger amount of data. The last issue is related to the objective of introducing best-practice guidelines for practitioners. The numerical test provides evidences on the effectiveness of curriculum learning, as this choice leads to consistent and accurate models even with a limited amount of data with respect to what is prescribed by ergodic theory. Moreover, the analyzed experiments have neatly shown that starting from the fixed points of the dynamical system constitutes a workaround in learning the dynamics of the system from situations where the memory is essentially harmless to the learning process. Indeed, in these memory-based models, flushing an initially incorrect memory is crucial for obtaining a relevant and faithful model able to generalize its prediction beyond the sole situations encountered during the training step. Our numerical experiments clearly show that there is an impact of the initialization of the memory on the quality of the predictions provided by the LSTM models and it is particularly severe if one learns from a set of short-termed observations. Superior performance was consistently achieved with a memory initialized at random to avoid bias. While we here focused on the Lorenz'63 system and LSTM models, we believe these findings to be widely applicable and to provide some evidence-based good practice for data-driven modeling of physical systems. Future work will focus on a finer analysis of the interplay between the structure of the dynamical system and the learning process, as well as on the choice of appropriate metrics complementing entropy levels, toward a principled learning strategy. In that sense, curriculum learning strategies can be envisioned as a valuable tool for the analysis of complex dynamical systems and open up avenues for research within the context of active learning.

  (a) Ergodic datasets: T ergo and T split . (b) Random sampled dataset: T ens . (c) Structured dataset: T FP .

Figure 2 :

 2 Figure 2: The figure qualitatively illustrates the different training sets used in this investigation. Each ergodic dataset (T ergo ) is composed by one long trajectory (a). The datasets with trajectories emanating from the fixed point, exploiting the structure of the attractor (T FP ) and the datasets randomly sapled (T ens ) are composed by an ensemble of 9 trajectories, of which one example for each is shown in (b) and (c), respectively (see also Sec. 4.2).

Figure 3 :

 3 Figure 3: Distribution of the SVD-entropy (Eq. 10) associated with the four different datasets described in Sec. 4.2: ergodic (T ergo ), ergodic-split (T split ), random sampled (T ens ) and structured (T FP ). Only the structured dataset shows substantially lower entropy, while the other datasets share the same values and distribution.

Figure 4 :

 4 Figure 4: Left column: Time-evolution of the y-component of the Lorenz'63 system, both from the true (true, solid red) and the learned model (prediction, dashed blue) from different training sets and initial conditions. Right column: zoom on the [0, 2] time units interval, with the mean prediction (dashed blue line) from an ensemble of 100 trained models, and the associated 1 standard deviation. Top row: The prediction of the best (out of the ensemble set) model trained from T short is shown. The initial condition is drawn from the training set. The prediction is seen to remain reliable up to about 4.3 Lyapunov times (5 time units) before drifting significantly from the truth. Middle row:when the same model is evaluated with an unknown initial condition (test set), the performance deteriorates and the accuracy is lost after less than 4 time units. When a unique ergodic-compliant trajectory T ergo is used to train a model (bottom row), the prediction from an unknown initial condition improves significantly.

Figure 5 :Figure 6 :

 56 Figure 5: Time evolution of the prediction analysed as cumulative 2-norm error, using an average of 100 models. A short windows of 1 time unit is considered. In the inset (a), the average predictions emanating from trajectories included in the training set are shown. In (b), the error is shown when analysing the test sets: larger deviations are observed for the random sampled datasets.

  (a) t-distributed stochastic neighbor embedding (t-SNE) of the parameters of the trained models analyzed in Fig. 6. Each dot corresponds to a different model and its size visually indicates the accuracy: the larger the dot, the larger the error with respect to d A A A . (b) Radial distribution of the t-SNE projected model parameters. r indicates the distance from the barycenter of the cluster in Fig.

Figure 9 :

 9 Figure9: Distribution of the SVD-entropy (Eq. 10) associated with the structured dataset T FP (blue curve) as compared with the entropy computed for the same dataset shortened by progressively removing samples in the vicinity of the fixed points. The shortest trajectories are characterized by only 1000 samples. The entropy progressively increases as the redundancy of the dataset reduces. The highest levels of entropy found for the latter case are comparable to the estimate computed for the reference dataset T ergo , as well as the cases T ens and T split .

Acknowledgement

This work was funded by the French Agence Nationale de la Recherche via the Flowcon project (ANR-17-ASTR-0022) and the Speed project (ANR-20-CE23-0025-01). L.M. gratefully acknowledges stimulating discussions with Alex Gorodetsky (University of Michigan, US). O.S. thanks Luca De Cicco (Politecnico di Bari, Italy) for exchanges on the role of entropy metrics in curriculum learning. S.C. gratefully acknowledge fruitful discussions with Angelo Vulpiani (University La Sapienza, Italy).