To myself I seem to have been only like a child, playing on the seashore, and diverting myself in now and then, finding a smoother pebble or a prettier shell than ordinary, whilst the great ocean of truth lay all undiscovered before me."

Isaac Newton

 and soliton or breather occurrence [

1

Examples of extreme events in nature and society . . . . . . . . . 2

Schematization of the experimental setup used in [Solli 2007] . . . 3

Histogram of the pulse heights as in [Solli 2007 Vortices and extreme event formation as in [Gibson 2016] . . . . . 7

Numerical results as in [Chaté 1999

Abstract

Extreme events are phenomena, often considered as catastrophic, that occur in the tail of a distribution usually deviating from an expected, exponential decay. In optics, these events were first studied in the context of fibers, where they have been extensively analyzed, as optical rogue waves, in light of the well known analogy between optics and hydrodynamics, through the nonlinear Schroedinger equation. With the development and the broadening of the field, extreme events have been also studied in dissipative optical systems with or without spatial degrees of freedom.

In this Thesis we focused on the study of extreme events in three different active and dissipative optical systems, each presenting one or two spatial degrees of freedom, either in the transverse plane, perpendicular to the direction of propagation of light, or in the propagation direction. Localized structures of different nature represent an important possible solution in each one of the systems here studied, hence their interaction and the role played in the formation of extreme events have been also investigated into details. In the first system, a monolithic broad-area semiconductor laser (VCSEL) with an intracavity saturable absorber, we report on the occurrence of extreme events in the 2D transverse plane of the electric field intensity. In particular we highlight the connection between these objects and cavity solitons, both stationary and oscillatory, also present in the system. In the second system, a highly multimode laser with optical injection spatially extended along the propagation direction, we analyze the interaction and merging of phase solitons, localized structures propagating along the cavity carrying a 2π phase rotation. Extreme events have been investigated in two configurations: a first one where they emerge from the collision of phase solitons with other transient structures carrying a negative chiral charge, and a second one where high-peak events emerge from an unstable roll regime where phase solitons are not a stable solution. In both these systems we investigate the role of chirality in the extreme event formation. In the third system, a broad-area semiconductor laser (VCSEL) with optical injection, we study into details the interaction of cavity solitons in the transverse plane, described as two particles subjected to an interaction potential exponentially decreasing with the distance between the two objects: a possible analogy with hydrophobic materials is here suggested. Some preliminary results showing spatiotemporal extreme events in this system are also given.

KEYWORDS: dynamics of nonlinear optical systems -semiconductor lasers -extreme events -rogue waves -dissipative solitons.

Résumé

Les événements extrêmes sont des phénomènes, souvent considérés catastrophiques, qui se produisent dans la queue d'une distribution généralement en s'écartant d'une décroissance attendue exponentielle. En optique, ces événements ont été étudié dans le contexte des fibres, où ils ont été amplement analysés, comme des vagues scélérates, par analogie bien connue entre l'optique et l'hydrodynamique à travers l'équation de Schroedinger non linéaire. Avec le développement et l'élargissement du domaine, l'étude des événements extrêmes a été étendue à des systèmes dissipatifs avec ou sans degrés spatiaux de liberté. Dans cette thèse on se concentre sur l'étude des événements extrêmes dans trois différents types de systèmes optiques actifs et dissipatifs, présentant chacun un ou deux degrés spatiales de liberté, soit dans le plan transversal (perpendiculaire à la direction de propagation de la lumière) soit dans la direction de propagation. Des structures localisées de nature différente constituent une solution possible importante dans chacun des systèmes étudiés ; leurs interaction autant que leurs rôle dans la formation des événements extrêmes ont donc été analysés en détails. Dans le premier système, un laser à semiconducteur monolithique (VCSEL) à large surface avec un absorbant saturable, on présente la formation d'événements extrêmes dans le plan transversal à deux dimensions de l'intensité du champ électrique. En particulier, on met en évidence la liaison entre ces objets et les solitons de cavité, soit stationaires soit oscillatoires, aussi présents dans le système. Dans le deuxième système, un laser multimodal spatialement étendu dans la direction de propagation avec injection optique, on analyse l'interaction et la fusion des solitons de phase, des structures localisées qui se propagent dans la cavité en transportant une rotation de phase de 2π. Les événements extrêmes ont été étudié dans deux configurations : une première où ils émergent de la collision des solitons de phase avec des autres structures transitoires transportant une charge chirale négative, et une deuxième où des événements d'intensité élevée émergent d'un régime instable de motif en rouleau où les solitons de cavité ne sont pas des solutions stables. Dans les deux systèmes, on examine le rôle de la chiralité dans la formation des événements extrêmes. Dans le troisième système, un laser à semiconducteur avec injection optique, on étudie dans les détails l'interaction des solitons de cavité dans le plan transversal, décrits comme deux particules soumises à un potentiel d'interaction décroissant exponentiellement avec la distance entre les deux objets : une analogie possible avec les matériaux hydrophobes a été proposée. Des résultats préliminaires présentant des événements extrêmes spatiotemporels dans ce système sont aussi donnés.

MOTS-CLÉS : Dynamique des systèmes non linéaires -lasers à semiconducteur -événements extrêmes -vagues scélérates -solitons dissipatifs.

Introduction Preface

This manuscript is focused on the study of extreme and abnormal events in active, dissipative optical systems, such as lasers, with a focus on extended cavity lasers. Extreme events have been a topic of rising importance in the past decades and, even though initially their main field of study in optics regarded passive systems such as fibers, the playground has significantly broaden to include a multitude of different areas. In particular lasers are found to be an interesting workbench for the study of these phenomena, since they could allow us to shed some light on the deterministic nature of the extreme events observed and, in the case of extended cavities, to investigate the role of spatial effects. The aim of this introduction is to place the topic of this Thesis in the context of extreme events, a subject of study in many different fields (i.e. society, economy, geology, weather forecast and so forth). Then we will illustrate the main results that prompted the field in the optical context and focus on its broadening to nonlinear active dissipative systems. Finally we will discuss the role of localized structures and their interaction in this framework and illustrate the content of the present manuscript.

Extreme events in different fields

There are many phenomena that fall or could fall under the definition of extreme events, depending on the point of view of a person [Jentsch 2006]. For instance an event is usually considered extreme according to its rarity and catastrophic nature but these can be subjective concepts, for instance to an insurance company an event is to be considered rare if it presents a low probability of occurrence and catastrophic if it implicates big consequences. Hence, in this case, the extremeness of an event is considered in terms of both its impact and its frequency. To a scientist instead the main charm regarding extreme events consists in their uniqueness, and in the concept of whether or not it is actually possible to predict these objects. An event is considered extreme if it involves a big deviation from a series of measurement [Jentsch 2006]. In this view for example a magnetic storm or a particularly violent earthquake can be considered as extreme events even if they do not, respectively, affect any device or occur in a inhabited area.

Extreme events in nature and society can affect a single individual, for example the contraction of a rare disease, as well as multiple people or the en- vironment, for example natural disasters, such as floods or hurricanes, societal disasters, such as pandemics, or market crashes.

Given the variety of possible phenomena to be considered as extreme the most obvious approach for their characterization consists in the study of their statistics and dynamics, focusing in particular on the possible analogies between different extreme events, their generating mechanisms and thus their predictability.

In general extreme events occur in the tail of a probability distribution. In many, but not all, cases [Jentsch 2006] the tails of a distribution containing extreme events are heavy, or subexponential, with exponential tails characterizing any Gaussian distribution. Given the probability distribution of an event it is possible to evaluate their probability of occurrence or even the probability of having any event larger than a certain threshold. Some statistical theories have been developed on this topic [Gumbel 1958], but their working assumptions for the events to be independent and identically distributed are very strong and rarely in good agreement with reality.

The other main approach in the study of extreme events consists in looking at the dynamical properties of the phenomenon. All extreme events in nature are usually expression of the complex dynamics of a system [Jentsch 2006]: they do not occur randomly but due to a dynamical mechanism that introduces huge excursions in the system trajectory, bringing it far from its normal state. Some examples of possible extreme events scenarios are those leading to deterministic chaos or turbulence [Jentsch 2006].

Even if no universal theory can be developed for extreme events, we may expect some similarities between different phenomena. We also would like to point out that these similarities can be both in the cause, i.e. the generating mechanism or the initial conditions, or in the effects. Similarities between dif-ferent classes of extreme events in their effects are in fact much more easy to observe since they usually cause some sort of destruction (i.e. victims in the case of natural or financial disasters, system damage in case of power spikes and so forth). For these reasons prediction of extreme events is a rather important task. One of the main problems regarding predictability for these objects consists in the fact that many models are built on mean values [Jentsch 2006], hence when trying to describe extreme events such models either underestimate or fail in their prediction.

For a quantity to be predictable some amount of determinism is needed [Jentsch 2006]. Nevertheless we have to remind that also in deterministic systems chaos can emerge, causing the impossibility of long-term predictions [Erkintalo 2015] even if this does not imply, on a short scale, the total lack of possible precursors.

Fiber optics and the analogy with hydrodynamics

A particular kind of extreme events, often studied in optics are rogue waves (RWs) which are also studied in oceanography and can be compared with this field, due to the well-known analogy between optics and hydrodynamics. Rogue waves are isolated high-amplitude waves occurring more frequently than expected from Gaussian statistics [Kharif 2009, Onorato 2013].

The concept of rogue waves has been infamously well-known in oceanography since centuries, initially as anecdote or sailors' tale and then, after the first experimental observations, like the Drapner wave in 1995 (see for example [Haver 2000, Kharif 2009]), as an actual topic of study. In oceanography these objects are also called monster or freak waves and are known to destructively damage ships. There are different definitions for RWs, the main ones defining them as waves whose trough-to-crest height exceeds either two times the significant wave height (i.e. the mean of the highest third of the wave heights detected), or the mean plus eight times the standard deviation of the data. More details on the derivation of these definitions and rogue waves statistics are given in Chapter 1.

Experimental measurements of rogue waves in oceanography present some difficulties [Onorato 2013]. The most used instrumentation consist in a buoy, registering the surface elevation data in a single point of the sea. Some evolutions of this technique include the possibility of collecting also the direction information through the use of a directional buoy (see for example [De Pinho 2004]) and an alternative detection method consists in the use of a wave radar [START_REF] Wyatt | [END_REF]. In any case this kind of single point measurement presents some issues, since a simple temporal trace is not able to fully describe the spatiotemporal complexity of the sea surface. Since, to get a physical understanding of RW formation, some kind of space-time measurement is required [Onorato 2013], some alternative techniques have been proposed, among which the use of synthetic aperture radars (SARs) [Schulz-Stellenfleth 2004], but are still regarded with skepticism by the community [Janssen 2006]. More details on the experimental evidence of rogue waves in the ocean can be found in [Kharif 2009, Ono-Fig. 2 -Reprinted with permission from [Solli 2007]. Schematization of the experimental setup used by Solli et al. for the first reported experimental observation of rogue waves in fibers.

rato 2013].

To understand the analogy drawn between optics and hydrodynamics regarding rogue waves we have to consider the one dimensional nonlinear Schroedinger equation (NLSE)

i ∂A ∂z + c 1 ∂ 2 A ∂t 2 + c 2 |A| 2 A = 0 .
In particular in hydrodynamics [Osborne 2010], with

c 1 = - k 0 ω 2 0 , c 2 = -k 3 0 ,
where k 0 is the wave number and ω 0 is the carrier frequency, the NLSE describes the evolution of the surface wave group envelope in deep water. In fiber optics [Agrawal 2013], with

c 1 = |β 2 | 2 2 , c 2 = γ ,
where β 2 < 0 is the group velocity dispersion and γ the nonlinear coefficient, the NLSE describes the evolution of a light pulse envelope propagating inside a fiber. Note that usually the NLSE in hydrodynamics has time and space switched, but this poses no peculiar obstacle to the analogy, except the fact that the coefficients c 1 and c 2 need to be adapted [Osborne 2010, Dudley 2014].

From an historical point of view [Randoux 2016b] the discovery of multisolitons as solutions of the NLSE [START_REF] Zakharov | [END_REF]] prompted the optical field to the study of new solutions of the NLSE and eventually lead to the discovery of a new class of objects, describing the instability of a finite coherent background: these solutions are usually known as breathers, and include the Akhmediev breather [START_REF] Akhmediev | [END_REF]] and the Peregrine soliton [START_REF] Peregrine | Water waves, nonlinear Schrödinger equations and their solutions[END_REF]], that is the Akhmediev breather limiting case for a period going to infinity [Kibler 2010]. These two solutions, in particular, fully describe the modulational instability process and were then suggested as rogue wave prototypes in oceanography [Dysthe 1999].

Fig. 3 -Reprinted with permission from [Solli 2007]. Histogram of the pulse heights, displaying a large deviation from the distribution of most stochastic processes [Solli 2007].

With the seminal paper by Solli [Solli 2007] the concept of rogue waves started to be extended also to the optical context and the experimental observation of breather solutions in both water wave tanks [Chabchoub 2011] and fibers [Kibler 2010[START_REF] Frisquet | [END_REF]] brought forward the NLSE analogy in a rigorous way.

In the work by Solli, using the setup schematized in Fig. 2, the authors inject a mode-locked ytterbium-doped fiber laser into a nonlinear fiber generating supercontinuum: rogue waves where found in the heights of the electric field intensity at the output of the fiber, collected through the use of a photodetector, after a filtering and signal enlargement process. A statistical analysis of the pulse heights showed L-shaped histograms largely deviating from the distribution of most stochastic processes [Solli 2007] as illustrated in Fig. 3. This pioneering paper has the merit of introducing the concept of rogue waves to the optical community. Ever since the field has notably broadened and specialized. In particular, regarding rogue waves in fibers, detection without the use of signal filters, has become more precise, allowing for a single-shot, realtime measurements [Suret 2016[START_REF] Närhi | Real-time measurements of spontaneous breathers and rogue wave events in optical fibre modulation instability[END_REF]] in propagation.

Among other works in the fiber supercontinuum context we mention [Erkintalo 2009[START_REF] Kibler | [END_REF], Mussot 2009]. Rogue waves have also been studied in the case of optical fiber propagation with [START_REF] Taki | Third-order dispersion for generating optical rogue solitons[END_REF], Conforti 2015] and without [Walczak 2015] higher order dispersion, in the case of laser filamentation [START_REF] Kasparian | [END_REF], passive cavities [Conforti 2015] and photorefractive ferroelectrics [Pierangeli 2015]. In the context of fibers but in case of an active system we mention RW study in Raman fiber amplifiers [Hammani 2008] and mode-locked fiber lasers [Zaviyalov 2012]. For a full review on the study of rogue waves and extreme events in different optical contexts please refer to [Onorato 2013, Dudley 2014, Akhmediev 2016].

In the optical context in particular, Refs. [Akhmediev 2011a, Akhmediev 2011b] identify as an early warning for rogue wave occurrence, the emergence of a triangular shape (in logarithmic scale) in the spectrum related to the formation of the Peregrine soliton and other higher-order solutions of the NLSE. In Ref. [Birkholz 2015] the authors use a novel approach, applying nonlinear time series analysis to search for traces of predictability directly in the experimental data. Through the use of an adapted version of the method proposed by Grassberger and Procaccia [Grassberger 1983], the authors are able to distinguish stochastic processes from processes where determinism prevails, making possible a prediction. In particular the method is based on dividing the time series in a sequence of smaller segments and computing the distances between the points of different segments: an elevated number of small distances implies similarities between different segments and can be interpreted as a signature of determinism, while in the stochastic case most distances will show large differences between different segments. The authors search for similarities close to RW occurrence and identify multifilamentation as a mainly deterministic process where rogue waves appear with a warning. The application of the same method to the experimental data of the Draupner wave and to the power fluctuation data in the tail of the supercontinuum generation by Solli [Solli 2007] does not lead to the same conclusions but some issues are to be addressed [Erkintalo 2015].

In general we can observe that, within the scientific community two parallel ways of exploring the analogy between optics and hydrodynamics have been outlined.

The first method focuses on the shape of the possible rogue waves: this approach bases itself on the powerful possibility of finding the exact solutions of the NLSE and has fundamental historical reasons. As a general reminder we would like to point out that the NLSE in hydrodynamics, and also in optics, does not describe the shape of individual waves but their envelope, hence all the specific envelope solutions of the NLSE, among which breathers and solitons, are not, in general, to be considered as individual rogue waves [Dudley 2014]. In recent works it has been observed that the shape of rogue waves, even in system configurations well described by the NLSE, might differ significantly from the well-known breather and Peregrine soliton solutions [Randoux 2016a]. Furthermore these specific solutions are generated in optical systems under very specific initial conditions, which result unlikely to be observed in the ocean [Randoux 2016b].

The second approach consists instead in finding similarities between the two fields from a statistical point of view [Onorato 2004, Erkintalo 2010]. In fact as we will illustrate more into details in Chapter 1, under the assumption of a Gaussian distribution for the sea surface elevation it can be demonstrated that the amplitude of the water wave envelope (and with good approximation its height) would follow a Rayleigh distribution. In optics the main quantity under analysis is the intensity of the electric field and it can be demonstrated that, assuming a Gaussian distribution for the amplitude of the real and imaginary part of the electric field envelope, a negative exponential distribution for the intensity of the electric field envelope is to be expected. Hence if in oceanography a series of data acquisitions presents a Rayleigh distribution for the amplitude of the water wave envelope and in optics an experimental setup configuration displays a negative exponential distribution for the intensity of the electric field envelope, the two systems are to be consider analog. Furthermore a positive deviation from such distributions can be identified in both fields as a signature of rogue waves in the system, admitted it is large enough to pass the rogue waves thresholds.

The analogy between optics and hydrodynamics through the use of the nonlinear Schroedinger equation is certainly a powerful tool. In order to correctly exploit its potential it is nevertheless necessary to treat such analogy within its limit. Since the NLSE has been reported as not able to describe ocean waves on short scales [Dyachenko 2005] we might expect other models, other than the NLSE, or an analogy on a statistical level to give us some additional insight on this topic. Furthermore a statistical approach could also help establish a link with the wider theory of extreme events in physics [Dudley 2014].

Extreme events in active optical systems

In the previous Section we have illustrated the development of the extreme event field in the context of fibers and some other optical systems. Most of the work on fibers is developed in conservative (or weakly dissipative) systems with propagative geometry. Lasers on the other hand are active, dissipative systems with possibly spatial degrees of freedom. Extreme events in this context have been observed in the configurations of optically injected [Bonatto 2011], modelocked [Kovalsky 2011, Zaviyalov 2012, Lecaplain 2012], Raman [Randoux 2012] and pump-modulated [Pisarchik 2011] lasers.

The study of extreme events in (active) nonlinear optical cavities presents usually a different focus with respect to the fiber case. First of all, in the laser case, the main approach used in the investigation of extreme events is statistical, since none of the systems here presented are describable through the use of the one-dimensional nonlinear Schroedinger equation. Far from being a obstacle, this issue offers us the opportunity to reconnect to extreme event theory and even to the concept of rogue waves from a statistical perspective, exploring the analogy between the two field on a different level. Most of research interest in this field is gathered around the investigation of the general physical or dynamical mechanisms at the base of extreme event formation. In particular in the case of lasers these mechanisms include spatiotemporal chaos [Selmi 2016], vortex turbulence [Gibson 2016], external crisis [Bonatto 2011, Zamora-Munt 2013], dissipative solitons [Rimoldi 2017b] and soliton interaction [Walczak 2017]. The broadening of the extreme event field to the context of lasers allowed to study the deterministic nature of RWs and the consequences of the spatial effect [Bonazzola 2013], an issue very relevant also in light of the analogy with oceanography.

In this Section we propose a summary of the results obtained in the systems more relevant to this Thesis. In particular, regarding the subject of lasers with saturable absorber (LSA), we would like to mention the main results of [Bonazzola 2013, Selmi 2016, Coulibaly 2017].

In [Bonazzola 2013] the authors study the importance of spatial effects in the description of an all-solid-state LSA. In particular, using a rate-equation model with one equation for the electric field intensity, three equations for the populations in the laser levels and in the ground state and one equation for the difference of population in the saturable absorber, the authors analyze numerically different regimes, among which self-Q-switching, period doubling and deterministic chaos, well described by the model. Extreme events have been observed experimentally in some of the chaotic regimes for this type of laser [Kovalsky 2011, Hnilo 2011] but have no matching numerical description with the rate-equation model, leading to the conclusion that spatial effects, not included in the model, actually play an essential role in the formation of extreme events for this kind of laser. A few years later, in [Selmi 2016, Coulibaly 2017], the presence of extreme events is investigated experimentally and numerically in the transverse section of a semiconductor laser with an intracavity saturable absorber in a configuration of spatiotemporal chaos. The model there used to describe the laser [Bache 2005] includes spatial coupling through a diffractive term in the second derivative of the spatial coordinate in the equation for the electric field. Limiting their analysis to a single transverse spatial dimension, the authors observe the presence of extreme events in the experiment and in the simulations for certain choices of parameters. The origin of the observed extreme events is identified in spatiotemporal chaos, characterized through the Lyapunov coefficient analysis and the Kaplan-Yorke dimension: in particular [Selmi 2016] a growth in the percentage of the observed extreme events is associated to an increase of the Kaplan-Yorke dimension. As a general reminder we would like to stress that, while a certain degree of spatiotemporal complexity might be considered a requirement for the observation of extreme events in this kind of system, de- terministic chaos is in general, at least experimentally, not a sufficient condition [Kovalsky 2011]. More details on the results of these two articles in comparison with the results shown in this Thesis will be given in Chapter 1.

In the context of lasers with optical injection we are going to focus on the study of systems with a different number of spatial dimensions [Bonatto 2011, Zamora-Munt 2013, Gibson 2016]. This will help us understand better if indeed spatial dimensionality affects the generating mechanisms of extreme events. Furthermore works in purely temporal (or with a very small transverse section) systems confirm that space is actually not a requirement for the observation of extreme events, allowing for a simpler study of predictability and the role of noise.

The role of spatial dimensionality has been also investigated in other optical context for example in cavities with liquid-crystal light valve as gain media [START_REF] Montina | Nongaussian statistics and extreme waves in a nonlinear optical cavity[END_REF]] and in the passive case, in photorefractive media [Pierangeli 2015] and photonic crystals [Liu 2015]. In particular in [START_REF] Montina | Nongaussian statistics and extreme waves in a nonlinear optical cavity[END_REF]] the nonlocal coupling between different spatial regions of the electric field intensity gives rise to symmetry breaking in the system and the presence of extreme events and in [Pierangeli 2015] the combination of disorder and nonlinearity appears to lead to extreme events observation. In [Liu 2015] extreme events are triggered in a photonic crystal through the control of a spontaneous synchronization mechanism that enlarges the frequency bandwidth leading to high-peak pulses. A picture of rogue wave formation in the spatial distribution of the electromagnetic energy is illustrated in Fig. 4.

In [Bonatto 2011, Zamora-Munt 2013] the authors study the setup of a singletransverse mode optically injected semiconductor laser (in particular a VCSEL). Given the very small laser transverse section, it is possible to describe the setup as a zero-spatial dimension, purely temporal system. Extreme events are experimentally obtained when the slave laser emission frequency is unlocked from the frequency of the master laser and passes to a chaotic regime through period doubling. Note that extreme events are here observed at the beginning of the chaotic region but, since for more chaotic parameter choices the probability of large excursions in intensity increases, the extreme events threshold moves towards higher values of intensity displaying less extreme events in the system. The experimental results appear in good agreement with the numerical ones, achieved with a class-B rate-equation model for a semiconductor laser. Since extreme events do not appear every time the system is in a chaotic regime but in a very specific region of the phase space, the authors investigate the properties of the chaotic regions needed for the occurrence of extreme events. As illustrated in Fig. 5(a), the system presents three fixed points two of which (S2 and S3) are unstable foci and one (S1) is a saddle point. To S2 are associated a stable (blue dotted line) and an unstable manifold (green solid line). When an extreme event is about to occur, the trajectory of the system moves along the stable manifold of S2, as illustrated in Fig. 5(b) and in the zoom in (c). Due to the presence of the two-dimensional unstable manifold of S2 while approaching the blue dotted line the trajectory of the system spirals out. Finally, since while approaching S2, which is the low intensity frequency-locked solution for the injected laser, the system has accumulated carriers, a high-peak pulse is emitted, returning the system to low values of carrier density. Every time the system approaches the stable manifold of S2 an extreme event is likely to occur. It is important to note that [Zamora-Munt 2013], since the amplitude of the electric field intensity excursion depends on the system trajectory (varying) entry conditions along the stable manifold of S2, this process is not to be classified as excitable. Further numerical analysis draws the authors to suggest that the mechanism at the origin of extreme event formation in this system consists in the collision of the attractor associated to the fixed point S3 with the stable manifold of the saddle point S1. In particular, in the chaotic regions where there are no extreme events, the stable manifold of S1 assumes the role of a barrier, preventing the system from approaching S2. Instead in the chaotic regions with extreme events the collision of such manifold with the attractor provokes an expansion of the attractor itself, allowing for the system trajectory to approach S2 and trigger extreme events.

Finally in [Zamora-Munt 2013], the role of noise has been investigated, in order to characterize the deterministic nature of the extreme events observed. In particular the authors notice that in regions where extreme events occur without noise, and are then deterministic, noise can be detrimental because it diminishes the chances of the system to approach the stable manifold of the fixed point S2. Instead in cases where deterministic extreme events do not occur (i.e. there is no collision between the attractor of S3 and the stable manifold of S1), noise can actually trigger stochastic extreme events, helping the system to overcome the barrier set by the stable manifold of the fixed point S1.

Another important study in the context of lasers with injection is that of a singly resonant optical parametric oscillator with external seeding [START_REF] Oppo | Self-organization, pattern formation, cavity solitons, and rogue waves in singly resonant optical parametric oscillators[END_REF]] and a class-A laser with optical injection [Gibson 2016], both in the case of a two- In both cases it is possible to describe the system through a single (complex) equation for the electric field. In particular, in the case of the optical parametric oscillator, the system is described by a forced complex Swift-Hohenberg equation, while, in the case of the class-A laser, the model is a forced Ginzburg-Landau equation, both with external driving. The spatial dependence is introduced through a transverse Laplacian term, for the diffraction, and a Swift-Hohenberg correction in case of the parametric oscillator [Lega 1994, Oppo 2009]. In both these systems the authors find, numerically, the presence of extreme events in an irregular regime, emerging for low values of injection, when patterned solutions such as hexagons (see Fig. 6(a)) become unstable. This specific regime of optical turbulence results dominated by interacting vortices, as shown in Fig. 6(b-d). In particular, when the patterned stationary solution of the system becomes unstable the system exhibits at first phase instability, with the trajectory in the Argand plane spreading on a limit cycle. Amplitude instability then occurs, where the trajectory displays large fluctuations from the limit cycle. Finally, where the amplitude fluctuations reach zero intensity, pairs of vortices of opposite charge start to appear. Vortex nucleation and annihilation continues then to occur and, since the turbulent state displays globally almost constant intensity, when zero-intensity moving vortices occur, high-amplitude spikes appear at the same time. For large vortex density, collisions between different vortices can occur, giving rise short-lived spikes, among which extreme events can be found. The number and position of the interacting vortices defines the shape of the spikes. The regime here illustrated is associated by the authors to the defect-mediated turbulence studied in [Coullet 1989].

Predictability has been investigated through different techniques, from the simple comparison of the power fluctuation time traces in presence of an extreme event (see for example [Zamora-Munt 2013]), and to a method of symbolic ordinal time-series analysis [START_REF] Alvarez | [END_REF]].

In [Zamora-Munt 2013] the authors superimpose the time traces of different extreme events, centered on the extreme event maximum: when approaching the extreme event, the intensity temporal profile shows a behavior common to the different time traces, giving the ability to predict extreme event occurrence before it actually happens. The addition of noise does not destroy this common behavior but instead it just reduces the window of predictability. Another method [START_REF] Alvarez | [END_REF]] uses symbolic ordinal time series analysis to investigate predictability. In particular, considering the intensity peak heights for each peak above a certain threshold, the method studies the height of the previous three peaks. Comparing the height of the three different peaks, a symbol, taking into account just the relative temporal ordering of the values, is applied to each three-valued segment. All the pattern combination probabilities are then plotted as a function of the varying threshold: the pattern combination with the highest probability, also for high threshold values, is the more likely to precede an extreme event. One drawback of this method consists in the fact that it does not predict the height of the extreme event following a specific pattern combination. Furthermore we would like to stress that the most likely pattern combination identified through this method is to be considered just as a precursor: in fact nothing prevents such pattern to occur also in other places of the time trace that do not lead to the formation of an extreme event.

Another possible predictability technique is the one based on the Grassberger-Procaccia algorithm [Grassberger 1983], used in the context of fibers (and oceanography) [Birkholz 2015] and already illustrated in the previous Section.

Cavity solitons and their interactions

In the context of fibers we have pointed out that one of the possible approaches to study extreme events and, in particular, rogue waves, consists in focusing on the shape of the solutions of the NLSE. Some of these solutions (for example the Peregrine soliton and breathers such as the Akhmediev breather) have been suggested as prototypes of rogue waves. It is important to remark that the solutions of the NLSE regard a conservative system, hence, for instance, the Peregrine soliton is a conservative localized structure. Collisions between conservative solitons or breather solutions have been also studied as possible mechanisms for the formation of rogue waves [START_REF] Frisquet | [END_REF]].

Conservative solitons represent a continuous family where all amplitude sizes can be obtained, in the weakly dissipative case instead, due to broken scale invariance, these structures are discrete members of such continuous family, being stable only for specific amplitude sizes [START_REF] Fauve | Solitary waves generated by subcritical instabilities in dissipative systems[END_REF]].

In the case of dissipative systems, such as lasers, localized structures can also appear. Cavity solitons (CS) belong to this category: in particular localized structures present themselves in situations where a physical quantity self-organizes in one or more isolated regions, surrounded by a qualitative-different, homogeneous state. Localized structures have been widely studied theoretically and experimentally since the 1980s [Koga 1980, Laughlin 1983, Rosanov 1988, Lugiato 2003, Akhmediev 2005, Akhmediev 2008].

Cavity solitons appear as bright intensity peaks that coexist with a stable homogeneous background. The term "soliton" describes in general solitary waves propagating without any change in shape. Given that solitons appear also in conservative systems, with the name of "cavity" solitons we are able to address the dissipative character of the system under study. Physically speaking CSs emerge from the compensation, on one side, of the effects of nonlinearity and diffraction and, on the other side, of gain and losses, since we are dealing with a dissipative system. This double compensation makes it possible for CSs to maintain a specific shape in space and time, independently from the boundary condition. Spatial cavity solitons were theoretically predicted in the transverse plane of semiconductor microcavities in [Brambilla 1997[START_REF] Michaelis | [END_REF], Spinelli 1998, Spinelli 2001]. The main advantages of generating CSs in these devices instead of the previous experimental implementations in macroscopic cavities and systems with optical feedback (see for example [START_REF] Saffman | [END_REF], Taranenko 1997, Schreiber 1997]), consists in the fast response of the devices and their compactness. The plasticity of cavity solitons, that is the possibility of switching on and off these objects (due to bistability), manipulate them independently and set in controlled motion, through the application of an external gradient, together with the additional advantages of developing these structures on a semiconductor device, makes cavity solitons very good candidates for optical information processing [Firth 1996, Brambilla 1997]. In particular spatially modulating a parameter of the setup can trap CSs in specific locations of the transverse plane, basically providing an optical memory array. The first experimental evidence of cavity solitons in a vertical-cavity surfaceemitting laser (VCSEL) with optical injection is reported in [Barland 2002], where the authors independently manipulate CSs on a diameter of ≈200 µm, through an external optical perturbation. Furthermore in [START_REF] Pedaci | [END_REF]] an optical memory array is realized through the spatial modulation of the holding beam phase.

In absence of an optical injection, cavity solitons can be observed in lasers with saturable absorber. Self-sustained cavity soliton emission presents a great advantage in terms of compactness and stability of the setup. In the present case the laser threshold can be subcritical, allowing for the coexistence of two stationary solutions, one of which is strictly greater than zero. This configuration implies optimal conditions for the observation of cavity solitons, since the background intensity is simply given by spontaneous emission, setting the contrast between CSs and the surrounding homogeneous state at its maximum. In this layout the system acts as a Cavity Soliton Laser (CSL), that is, a laser, homogeneously pumped along its transverse section, that emits localized structures surrounded by spontaneous emission [Bache 2005, Prati 2007a, Aghdami 2008, Prati 2010a]. One of the main differences regarding CSs with and without optical injection consists in the fact that, in the forced system, phase symmetry is broken and the Introduction CSs are locked to the phase of the holding beam, which is not the case for the system without optical injection, where phase symmetry is preserved.

There are two different experimental realization of a CSL in semiconductor materials through the use of a saturable absorber. The first realization exploits an experimental setup with two VCSELs in a face-to-face configuration, one acting as amplifier and the other acting as absorber [Genevet 2008]. The second experimental realization consists in a monolithic VCSEL with intracavity saturable absorber [START_REF] Elsass | [END_REF]], which results very close to the theoretical model elaborated for this kind of system [Bache 2005]. Another type of CSL [Tanguy 2006, Tanguy 2008] exploits frequency-selective feedback [Giudici 1999] from a Bragg grating.

In comparison with the spatial localization here illustrated, temporal localization along the propagation direction of nonlinear cavities can lead to the possibility of temporal cavity solitons. These objects rise from the compensation of nonlinearity and chromatic dispersion [Grelu 2012] and can occur both in system with [Leo 2010] and without optical injection [Grelu 2012]. The possibility of generating solitons at arbitrary distances inside the cavity, presents once again an interesting development in the context of information storage, with temporal cavity solitons as optical bits living indefinitely inside the cavity. In particular the observation that temporal solitons in microresonators are described by frequency combs [Del'Haye 2007[START_REF] Herr | [END_REF], spanning large portions of an octave, in the frequency domain attracted ever since a lot of interest in the study of these objects.

Three-dimensional localized structures, also known as spatiotemporal solitons or light bullets, represent an interesting possible merging of the properties of the two kinds of localization here presented. Three-dimensional dissipative solitons have been predicted in prototypical nonlinear resonators [START_REF] Tlidi | Three-Dimensional Optical Crystals and Localized Structures in Cavity Second Harmonic Generation[END_REF], Brambilla 2004[START_REF] Jenkins | [END_REF]] and nonlinear mirrorless configurations [Vladimirov 1999]. Experimental observations of this phenomenon are, for now, still lacking.

Given the potential application of CS studies to the context of information encoding, the interaction between these structures is of fundamental interest: in order to be treated independently as bits of information, cavity solitons need to exhibit no interaction effect. Hence the presence of a critical distance below which CSs tend to merge defines the maximum information density of a potential optical memory array. Regarding temporal cavity solitons in a fiber laser with a passive resonator driven by an external field, the critical temporal distance at which two CSs start to interact was shown to be around 40 ps, implying an information storage density of 125 bit/m [Leo 2010]. In more general and earlier studies on spatial localized structures where the nonlinear medium was a bistable interferometer or a collection of two-level atoms [Rosanov 1990, Brambilla 1996] two critical distances, d 1 and d 2 , relevant for cavity soliton interaction, have been found. In particular for an initial distance d between the two objects smaller than d 1 the two CSs would merge or annihilate, for d larger than d 2 they would not interact and for d 1 < d < d 2 the Fig. 7 -Reprinted with permission from [Chaté 1999]. Simulation results of a forced complex Ginzburg-Landau equation corresponding to different forcing and detuning with as a common initial condition a phase jump of 2π.

two CSs would repel till they reach a distance d = d 2 and then stop interacting. These results have been later confirmed in [Tissoni 1999] for a bulk semiconductor model, except one or more (depending on the sign of the detuning) equilibrium distances d * in between d 1 and d 2 were found, making the interaction between the two CSs either attractive or repulsive until the condition d = d * was reached. Finally in [Vladimirov 2002] clusters of localized structures are studied in the transverse plane of driven optical cavities below transparency, through the calculation of an interaction potential. In particular it was shown that the interaction between two nearby CSs implies a non-Newtonian motion between the two structures, with a velocity proportional to the perturbation generated by one soliton onto the other.

A special kind of localized structures

A special kind of localized structures characterized by a phase jump, has been predicted in extended oscillators in presence of spatially periodic forcing [Coullet 1986]: the origin of these objects is, according to the author, the result of a transition between commensurate states, of uniform or modulated type, and incommensurate states, with soliton-like phase disturbances. The same kind of localized structure has been studied in the context of a forced complex Ginzburg-Landau model, in one-dimensional, symmetric space [Chaté 1999]. As an example some of the results in this paper are illustrated in Fig. 7. In particular Fig. 7 shows different simulation results obtained when setting as initial condition a phase jump of 2π. In the case on the left side the phase jump propagates in a stable behavior, in the case in the middle, for a larger forcing, the solution becomes unstable and, in the case to the right, for a larger detuning and forcing, spatiotemporal intermittency is observed. In [Coullet 1998] the authors describe this kind of structures as optical excitable waves, occurring in the 2D transverse plane of a laser with optical injection. A similar kind of excitable wave, along the propagation direction, has been predicted in [Longhi 1998], where the author shows that soliton solutions of the Sine-Gordon equation are responsible for the localized structures observed in forced extended oscillators.

Motivation and contents

The purpose of the work here presented consists in studying extreme and abnormal events as well as localized structures interaction in three different active optical systems. The first system considered, studied in Chapter 1, is a broad-area semiconductor laser with an intracavity saturable absorber, like the one used in experiments [START_REF] Elsass | [END_REF]. The previous and contemporary work by other groups on extreme events in this context, illustrated in [Selmi 2016, Coulibaly 2017], has been developed in one transverse spatial dimension. In our case we instead consider a system with a two-dimensional transverse plane where we will characterize all the possible solutions available, including cavity solitons. The study of extreme events in the chaotic regime has required the development, for the first time, of an algorithm for the detection of spatiotemporal maxima in three dimensions (2D+time). Through a detailed statistical analysis we will show the presence of extreme events in this system for specific parameter regions. The inclusion of 2D spatial effects in our simulations plays a crucial role in the formation of extreme events, as suggested in [Bonazzola 2013], and it is at the origin of the discrepancies with the one-dimensional studies. Finally the observed striking similar spatial and temporal shape of the extreme events and stationary and self-pulsing CS will imply the enhancement of extreme events when the system is close to the cavity soliton attractor. Future development of this work consist, on one side, in the investigation of predictability, through the application of techniques such as the Grassberger-Procaccia algorithm [Grassberger 1983, Birkholz 2015] and the method of symbolic ordinal time series analysis [START_REF] Alvarez | [END_REF]], and, on the other hand in the study of the conservative limit of this system and its possible implications in the observation of extreme events.

Chapter 2 is devoted to the study of a highly multimode semiconductor laser with optical injection, alternatively in a ring or Fabry-Perot configuration. This system, being 1D+time, spatially extended along the propagation direction, builds a bridge between the cases, illustrated in the present Introduction, of lasers with injection in zero-spatial dimension [Bonatto 2011, Zamora-Munt 2013] and with a two-dimensional transverse plane [Gibson 2016]. A full derivation of the model used to describe this kind of laser and its comparison with other possible system descriptions are fully analyzed in the text. Interestingly enough, of the two mechanisms for the formation of extreme or abnormal events observed in this laser, one presents striking dynamical similarities with the zero-dimensional case. Furthermore some parallels on the role of chirality in the formation of these high-peak events, with the emergence of extreme events in vortex turbulence of [Gibson 2016], can also be drawn. Interaction between chiral localized structures in the laser is also studied, leading to the merging of these objects in soliton complexes. Instead the collision of the same localized structures with some transitory state carrying an opposite chiral charge leads to the formation of another kind of extreme events, once again highlighting the role of chirality in the collision process [Gibson 2016]. Finally the last system considered, in Chapter 3, is a broad-area semiconductor laser with optical injection pumped above threshold. In this case the interaction of cavity solitons is extensively studied, opening the way to a possible analogy with hydrophobic materials. One main difference with respect to similar systems below transparency [Tissoni 1999] consists in the absence of a critical distance for which cavity solitons stop interacting. Interestingly a similar kind of long-range interaction has been observed in the case of the chiral localized structures illustrated in the previous Chapter. Extreme events in a class-A laser with optical injection have been already reported in [Gibson 2016]. In the present case we will consider instead a semiconductor class-B laser and search for extreme events in a chaotic regime for low values of optical injection. Some preliminary results of this analysis will be here shown but further investigation is still in order.

The study of extreme events in active optical systems represents an interesting challenge and opportunity to address the dynamical mechanisms behind the formation of these objects. This investigation will in time allow for a general deeper understanding of the systems considered and their dynamics. Finally the connection of this field to the wider topic of extreme events in physics presents an interesting development in terms of possible analogies, given the high capability to observe, control and study these events in an optical system.

Chapter 1

Broad-area semiconductor laser with saturable absorber Preface

This Chapter is devoted to the description of a semiconductor laser (specifically a VCSEL, Vertical Cavity Surface Emitting Laser) with an intracavity saturable absorber. In particular, after a brief description of the system, the model used to represent it and an analysis of the possible attractors for this kind of laser, we are going to concentrate our attention towards the analysis of the extreme events observed numerically in the transverse plane of laser, perpendicular to the direction of propagation, in a chaotic regime. In the final part of the Chapter we are going to draw a comparison of extreme events with the other states possible for this kind of system, namely stable, oscillatory and chaotic cavity solitons. The main results here illustrated have been published in [Rimoldi 2017b].

Contents

Introduction

A laser with a saturable absorber [Antoranz 1982] as the one here studied can be schematized as a laser whose cavity contains, together with the active medium, with population inversion, also an absorber, or passive medium. In this kind of laser there is no external optical injection as the coherent light is generated inside the same cavity. The presence of the absorber grants the possibility of obtaining bistability in the system between the trivial laser-off state, which results stable up to a threshold, and the nontrivial stationary solution. The single mode laser with saturable absorber allows, given certain choices of parameters, for passive Q-switching in the region where both stationary solutions result unstable. Semiconductor lasers [Agrawal 1986] are one of the classes of lasers most commonly used in applications. VCSELs (Vertical Cavity Surface Emitting Lasers) in particular have been of rising importance in the past forty years. These devices [Wilmsen 2001] are grown by placing layers of different materials on the substrate: the gain region lays between two mirrors (one at the top and one at the bottom), which are both distributed Bragg reflectors containing about twenty alternating λ/4 layers of semiconductor material to produce a reflectivity > 99%. In comparison with the other conventional semiconductor lasers, as edge-emitting lasers, VCSELs present the peculiarity of having a cavity axis perpendicular to the mirrors. It is for this reason that such a high reflectivity is required for the laser to work properly, as in a single passage the light gain results pretty small, since the active region, consisting of a few quantum wells, is crossed perpendicularly by light involving a length of just ≈ 10 nm. Typically the diameter of a VCSEL is around 10 µm, but large diameter VCSELs up to 200 µm have been produced and analyzed extensively [Lugiato 2015] for the study of cavity solitons (see [START_REF] Elsass | [END_REF]] and the European projects PIANOS and FunFACS 1 ). Comparatively the height of a VCSEL is around a few microns, the laser cavity being almost completely occupied by the two mirrors. There are multiple advantages in the use of VCSELs , namely their compactness, their low cost and low threshold. Furthermore it is very useful to observe 1 www.funfacs.org that the shortness of the cavity in these devices implies that no more than one longitudinal mode can oscillate, and, in particular for large-area devices, the transverse modes result similar to those of lasers with a planar cavity. Finally, regarding the way we model this kind of lasers, we can observe that the low transmission limit is here very well fulfilled, and, since the active region is much smaller than a wavelength, we can also neglect the grating of the active medium and reduce the equations to those of a ring cavity. The fact that large-area VCSELs can be well modeled in low-transmission limit, single longitudinal mode approximation and planar transverse mode approximation makes them the perfect candidates for the study of cavity solitons (CSs) in a broad-area cavity with planar mirrors [Lugiato 2015].

The model here analyzed has been widely used by our group for the study of cavity solitons [Bache 2005, Prati 2007a, Prati 2010a, Tissoni 2011, Vahed 2014, Turconi 2015]: in particular it has been casted as illustrated in the next Sections to well describe the experimental setup at the laboratory of C2N in Paris [START_REF] Elsass | [END_REF], Selmi 2016, Coulibaly 2017]. In this case, as depicted in Fig. 1.1 the VCSEL presents three quantum wells, two of which are pumped and act as gain and a third one with no pumping that acts as the saturable absorber. The system presents then an intracavity saturable absorber, in contrast with other experimental implementations where the saturable absorber element is inserted through the use of a semiconductor saturable absorber mirror (SESAM) outside the cavity [Keller 1996, Keller 2006], which implies for the introduction of a delay into the equations. Here on the contrary, no delay term has to be considered. The model we consider consists in a set of rate equations which is derived from the Maxwell-Bloch equations [Fedorov 2000] where we have adiabatically eliminated the polarization. The spatial dependence is introduced by the transverse Laplacian term that presents both a purely imaginary part, which accounts for the diffraction (in the paraxial and slow varying envelope approximation), and a real part, which is introduced phenomenologically to maintain a finite linewidth for the gain and plays the role of a diffusive term. Furthermore in the equations for the carrier population in the active and the passive medium a radiative recombination term has been introduced through the years [Prati 2007a] for a better description of the experimental setup, allowing to find CSs for a linewidth enhancement factor larger than zero in the absorber and in cases where the relaxation rates and the saturation intensities both in the carrier populations of the amplifier and the absorber present similar values. Here we are going to use the same model with a different purpose, which is analyzing and describing the spatiotemporal chaos observed both numerically and experimentally in this kind of lasers, in pursuit of characterizing the extreme events emerging from this regime.

Model

The model used to describe the system of a monolithic broad-area VCSEL with an intracavity saturable absorber, of which a schematic representation is plotted in Fig. 1.2, is given by the following set of rate equations [Fedorov 2000, Bache 2005, Prati 2007a], which is derived from the Maxwell-Bloch equations when the polarization is adiabatically eliminated:

Ḟ = [(1 -iα)D + (1 -iβ) D -1 + (d + i)∇ 2 ⊥ ]F , (1.1a) Ḋ = b[µ -D(1 + |F| 2 ) -BD 2 ] , (1.1b) Ḋ = rb[-γ -D(1 + s|F| 2 ) -B D2 ] . (1.1c)
Where F is the slowly varying envelope of the electric field, D and D are the population variables related respectively to the carrier density in the amplifier and in the absorber (β), µ (γ) is the pump (absorption) parameter and s is the saturation parameter. Time is scaled to the photon lifetime (≈ 10 ps), b is the ratio of the photon to carrier lifetime in the amplifier and r is the ratio of the carrier lifetimes in the amplifier and in the absorber. Space is instead scaled on the diffraction length (≈ 4 µm). d is a diffusion coefficient included into the equations to take into account the gain finite linewidth. Although this term has here been introduced phenomenologically, a detailed derivation can be found in [Fedorov 2000], and comes from a nonstandard adiabatic elimination of the polarization variable. Finally B is a coefficient of radiative recombination, assumed equal for the two carrier densities [Prati 2007a]. For more details on the definition of these parameters see Appendix A as well as [Vahed 2012].

The dynamical equations are integrated through a Fourier split-step method with periodic boundary conditions and the spatial grid considered measures 256 × 256 pixels. More details about the integration method are given in Appendix B.

2 Please note that here and throughout the Thesis we have adopted the convention for which dotted variables represent the time derivative of such variables, ∂ i A represent the partial derivative of a quantity A with respect to i and ∇ 2 ⊥ is the transverse Laplacian

∇ 2 ⊥ = ∂ 2 x + ∂ 2 y .

HSSs and linear stability analysis

The model in Eqs. (1.1) presents the following trivial homogeneous stationary solution (HSS)

F s,0 = 0 , D s,0 = 1 + 4Bµ -1 2B , Ds,0 = √ 1 -4Bγ -1 2B , (1.2)
and the nontrivial HSS

F s,nt = √ I s e -iω s t , (1.3a) D s,nt = (1 + I s ) 2 + 4Bµ -(1 + I s ) 2B , (1.3b) Ds,nt = (1 + sI s ) 2 -4Bγ -(1 + sI s ) 2B , (1.3c)
with I s such that

D s,nt + Ds,nt -1 = 0 , (1.4)
and laser frequency ω s as

ω s = αD s,nt + β Ds,nt . (1.5)
Deriving µ from Eq. (1.3b) and using Eq. (1.4), we can write the pump parameter µ as a function of the stationary homogeneous intensity I s as follows: Broad-area semiconductor laser with saturable absorber

For the linear stability analysis we are going to recast Eqs. (1.1) in terms of the variables ρ, φ such that F = ρe i(φ-ω s t) , obtaining

ρ = (D + D -1)ρ -2 ∂ x ρ∂ x φ + ∂ y ρ∂ y φ -ρ∇ 2 ⊥ φ +d -ρ (∂ x φ) 2 + ∂ y φ 2 + ∇ 2 ⊥ ρ , (1.7a) ρ φ = -ω s ρ -(αD + β D) ρ -ρ (∂ x φ) 2 + ∂ y φ 2 + ∇ 2 ⊥ ρ +d 2 ∂ x ρ∂ x φ + ∂ y ρ∂ y φ + ρ∇ 2 ⊥ φ , (1.7b) Ḋ = b[µ -D(1 + ρ 2 ) -BD 2 ] , (1.7c) Ḋ = rb[-γ -D(1 + sρ 2 ) -B D2 ] .
(1.7d)

We now study the stability of the stationary homogeneous solution considering small spatiotemporal perturbations from the stationary values of F s , D s and Ds , modulated in the transverse plane by a wavevector k = (k x , k y ) and exponentially growing (decaying) in time, so that:

    ρ φ D D    =     ρ s φ s D s Ds     + e λt+i(k x x+k y y)     δρ δφ δD δ D    , (1.8)
with ρ s = √ I s , φ s , D s and Ds HSSs of the model in Eq. (1.1), satisfying Eqs. (1.2) for the trivial HSS and Eqs. (1.3), (1.4) and (1.5) for the nontrivial HSS with φ s = -ω s t. Inserting (1.8) into Eqs. (1.7) we obtain the following linearized equations for the perturbations:

λδρ = (δ D + δD) ρ s + (D s + Ds -1) δρ + k 2 x + k 2 y (ρ s δφ -dδρ) , ρ s λδφ = (αδD + βδ D) ρ s -k 2 x + k 2 y (dρ s δφ + δρ) , λδD = -b 2D s ρ s δρ + 1 + 2BD s + ρ 2 s δD , λδ D = -rb 2 Ds sρ s δρ + 1 + 2B D + sρ 2 s δ D , in matrix form, λ     δρ δφ δD δ D    = J LSA     δρ δφ δD δ D    (1.9)
where J LSA is the Jacobian matrix

     D s + Ds -1 -dk 2 ρ s k 2 ρ s ρ s -k 2 ρ s -dk 2 -α -β -2bD s ρ s 0 -b 1 + 2BD s + ρ 2 s 0 -2rb Ds sρ s 0 0 -rb 1 + 2B Ds + sρ 2 s      , with k 2 = |k| 2 = k 2 x + k 2 y .
Eq. (1.9) admits non-zero eigenvectors if and only if det(J LSA -λI ) = 0 , (1.10) with I is the identity matrix. Eq. (1.10) gives a fourth degree polynomial in λ (third degree in case of the trivial solution) (1.11) where the coefficients c i , with i = 0, 1, 2, 3, depend on the squared modulus of the wave vector k 2 , the system parameters, the electric field intensity.

λ 4 + c 3 λ 3 + c 2 λ 2 + c 1 λ + c 0 = 0 ,
A system results stable if and only if the real parts of all the eigenvalues λ of the characteristic equation (1.10) are negative. The Routh-Hurwitz stability criterion [Routh 1877, Hurwitz 1895], illustrated more in details in Appendix C, gives a necessary and sufficient condition for the real part of the eigenvalues to be negative, exploiting the coefficients c i . In particular it states that the necessary and sufficient condition for the real eigenvalues to be negative is that the coefficient c 0 satisfies c 0 > 0 , (1.12) instead, regarding the complex eigenvalues, the necessary and sufficient condition for their real part to be negative is that

c 2 1 + c 0 c 2 3 -c 1 c 2 c 3 > 0 , (1.13)
and for all the coefficients to be c i > 0. If for certain choices of parameters and/or for certain values of k 2 and ρ 2 s the condition (1.12) is not satisfied then we have a stationary instability, where the system will evolve toward a new stationary state. On the other hand if the condition (1.13) is not satisfied we talk about Hopf or dynamical instability, because the final state will not be stationary but instead oscillatory or chaotic.

Bistability

To study the bistability of the system we are going to consider the laser very close to the threshold µ thr in Eq. (1.14). The Taylor series expansion of Eq. (1.4) at the first order in I s is the following

D s + Ds -1 = 0 D s,0 + Ds,0 -1 + ρ 2 s 2B -(1 + s) + 1 1 + 4Bµ + s √ 1 -4Bγ = 0 ,
where we remind that if we set D s,0 + d s,0 -1 = 0 alone then we obtain µ = µ thr with µ thr as in Eq. (1.14). The coefficient of the laser intensity is generally small: for this reason inside the square parenthesis we can approximate µ with its threshold value µ thr which is equivalent to write

1 + 4Bµ ≈ 2(1 + B) -1 -4Bγ .
After some algebra we can then observe that the value of the pump parameter close to the threshold can be approximates as follows µ = µ thr -aI s

Broad-area semiconductor laser with saturable absorber with µ thr as in Eq. (1.14) and a as

a = 1 -Ds,0 + 1 2B -(1 + s) + 1 2(1 + B) - √ 1 -4Bγ + s √ 1 -4Bγ = (1 + s) 1 + B - √ 1 -4Bγ 2 + (1 + B) [2s -(1 + s)(1 + B)] 2B √ 1 -4Bγ + (1 -s) 2B .
We can notice [Prati 2007a] that for s = 1 the expression for a simplifies as

a = 1 + B √ 1 -4Bγ 2 -B(1 + B) B √ 1 -4Bγ .
If a > 0 the laser bifurcation is subcritical, otherwise if a < 0 it is supercritical.

In particular we find that the bifurcation is subcritical for γ > γ min with γ min given by

γ min = [2s + (1 + 2B) (1 + s)] (1 + 2B) 2 (1 + s) 2 -4s 8B(1 + s) 2 + (1 + 3s)(1 -s) 8B(1 + s) 2 - B 2 -1 + 1 2(1 + s)
, which, for s = 1 gets reduced to 1 2

(1 + B) 3 B - B 2 - 3 4 .
Furthermore the presence of the square root in the expression for a introduces an upper limit for bistability, given by

γ max = 1 4B .
We observe then that the quadratic recombination term introduces the possibility of having bistability even for s = 1, in contrast with the condition s > 1 + 1/γ that applies when such term is neglected [Bache 2005].

Plane-wave instability

Let us now consider first the case of plane wave instabilities: to do so we neglect any spatial dependence (in the form of k 2 ) in the perturbations. We observe that in this case, the second column of the Jacobian matrix J LSA contains only zero elements: this means that, even if we have introduced a perturbation on the phase δφ, such perturbation does not affect the evolution in time of the other variables ρ, D, D, furthermore the phase evolution in time depends only on the perturbations δρ, δD and δ D. Hence we can completely neglect the second row in the Jacobian matrix in the stability analysis of both the trivial and the nontrivial solutions, and J LSA gets reduced to a 3×3 matrix.

For the trivial solution, substituting Eqs. (1.2), with ρ s = 0 in J LSA , we obtain

J LSA,0 =   D s,0 + Ds,0 -1 0 0 0 -b (1 + 2BD s,0 ) 0 0 0 -rb (1 + 2B Ds,0 )   ,
with D s,0 and Ds,0 as in Eqs. (1.2). The eigenvalues of J LSA,0 are immediately found as

λ 1 = D s,0 + Ds,0 -1 , λ 2 = -b (1 + 2BD s,0 ) , λ 3 = -rb (1 + 2B Ds,0 ) .
The eigenvalues λ 2 and λ 3 are always negative, instead for λ 1 > 0 we obtain the laser threshold [Prati 2007a]

µ thr = (1 + B) 1 + 1 - √ 1 -4Bγ B -γ , (1.14)
which gets reduced to µ thr = 1 + γ in the case B = 0 [Bache 2005], hence the trivial solution (1.2) is stable for µ < µ thr .

For the nontrivial solution, substituting Eqs. (1.3) in J LSA , we obtain the following

J LSA,nt =   0 ρ s ρ s -2bD s ρ s -b 1 + 2BD s + ρ 2 s 0 -2rb Ds sρ s 0 -rb 1 + 2B Ds + sρ 2 s   ,
where we have rewritten for simplicity D s,nt = D s and Ds,nt = Ds , defined as in Eqs.

(1.3) with ρ s = √ I s such that Eq. (1.4) results valid. The characteristic equation (1.10) gives a third degree polynomial

λ 3 + c 2 λ 2 + c 1 λ + c 0 = 0 ,
with the following coefficients

c 0 = 2rb 2 D s + Ds s + 2BD s Ds (1 + s) + sρ 2 s ρ 2 s , c 1 = 2b(D s + Ds rs)ρ 2 s + rb 2 1 + 2BD s + ρ 2 s 1 + 2B Ds + sρ 2 s , c 2 = b 1 + r + (1 + rs)ρ 2 s + 2B(D s + Ds r) .
Exploiting the Routh-Hurwitz criterion (see Appendix C), we search for the plane-wave stationary instability boundary setting c 0 = 0. We observe that c 0 presents the same sign as the derivative of Eq. (1.6) with respect to the intensity

I s = ρ 2 s : in fact c 0 dµ (ρ 2 s ) /dρ 2 s = 2rb 2 ρ 2 s (1 + sρ 2 s ) 2 -4Bγ ,
which is always positive. In particular c 0 presents then a negative value as long as µ ρ 2 s has a negative slope, the stationary instability boundary lies always on the left turning point of the stationary homogeneous curve in Fig. 1.3(a) and the negative slope branch of the nontrivial solutions results always unstable. Regarding the plane-wave Hopf instability of the nontrivial solution we observe that for B = 0 the stationary solutions for the carrier densities become

D s,B=0 = µ 1 + ρ 2 s , Ds,B=0 = - γ 1 + sρ 2 s .
Following the Routh-Hurwitz criterion, we observe the the plane-wave Hopf instability boundary is given by the condition c 1 c 2c 0 = 0. Neglecting the terms ≈ rb 2 [Powell 1971] we obtain:

c 1 c 2 -c 0 = 0 2b 2 µ -r 2 γs ρ 2 s = 0 µ = r 2 γs ,
In particular for µ r 2 γs = µ h the system results Hopf unstable. In the case of B = 0 the expression for µ h is not as simple and an approximated expression valid for B γ is [Prati 2007a]

r c ≈ 1 + 3 2 B γ 1/3 + 1 8 B γ 2/3 + 53 96 B γ .
We observe then that, since B and γ are always positive, differently from the model without radiative recombination [Bache 2005], for r = 1 and s = 1 the system can be bistable without presenting any plane wave Hopf instability, thus allowing cavity solitons.

Pattern-forming instabilities

We are now going to analyze the stability of the nontrivial HSS in Eq. ( 1.3) with respect to spatially modulated perturbations. From the characteristic equation (1.10) we obtain a forth degree polynomial in λ as in Eq. (1.11), whose coefficients are:

c 0 = rb 2 1 + d 2 F (ρ s ) k 4 + -2sρ 4 s (αD s + β Ds ) -2 [αD s + β Ds s + 2B(α + sβ)D s Ds ] ρ 2 s + 2dG (ρ s ) k 2 , c 1 = b 1 + d 2 H(ρ s )k 4 + -2 (αD s + β Ds rs) ρ 2 s +d 2(D s + Ds rs)ρ 2 s + 2brF(ρ s ) k 2 + 2brG(ρ s ) , c 2 = 1 + d 2 k 4 + 2bdH(ρ s )k 2 + 2b(D s + Ds rs)ρ 2 s + rb 2 F(ρ s ) , c 3 = bH(ρ s ) + 2dk 2 with F(ρ s ) = 1 + 2BD s + ρ 2 s 1 + 2B Ds + sρ 2 s , G(ρ s ) = D s -Ds + Ds (1 + 2BD s )(1 + s) + sρ 2 s ρ 2 s , H(ρ s ) = 1 + r + (1 + rs)ρ 2 s + 2B(D s + Ds r) .
Exploiting the Routh-Hurwitz criterion we investigate first the stationary instability associated to the real eigenvalues setting c 0 = 0. We obtain the following condition on the square modulus of the transverse wave vector

k 2 = 2ρ 2 s (αD s + β Ds )(1 + sρ 2 s ) + 2BD s Ds (α + sβ) + β Ds (s -1) (1 + d 2 ) (1 + 2BD s + ρ 2 s ) (1 + 2B Ds + sρ 2 s ) -2dρ 2 s D s + Ds s + 2BD s Ds (1 + s) + sρ 2 s (1 + d 2 ) (1 + 2BD s + ρ 2 s ) (1 + 2B Ds + sρ 2 s )
which for s = 1 gets reduced to

k 2 = 2ρ 2 s (αD s + β Ds )(1 + ρ 2 s ) + 2BD s Ds (α + β) (1 + d 2 ) (1 + 2BD s + ρ 2 s ) (1 + 2B Ds + sρ 2 s ) -2dρ 2 s 1 + 4BD s Ds + ρ 2 s (1 + d 2 ) (1 + 2BD s + ρ 2 s ) (1 + 2B Ds + sρ 2 s )
, that, for d = 0, corresponds exactly to the expression found in [Prati 2007a]. This boundary condition is plotted in Fig. 1.4(a) for a specific choice of parameters, which will be the similar for most of the simulations here considered: as we can observe the upper branch of the stationary homogeneous solution is stationary unstable. The (Hopf) instability associated with the complex eigenvalues for r > r c is the only one present in the plane wave limit, leading to Q-switching. More in general for k 2 = 0 this instability is related to the transition from stationary cavity solitons to oscillatory cavity solitons, both of which we are going to study more into details in the next Sections. In Fig. 1.4(b) we show the Hopf instability domain obtained with the Routh-Hurwitz criterion c 2 1 + c 0 c 2 3c 1 c 2 c 3 < 0, for r = 1.75, d = 0.01 (black solid line) and the same parameters as in 1.4(a): we can notice that the grey shaded region is already stationary unstable. Finally in Fig. 1.4 the pattern instability boundaries for d = 0 are also plotted, respectively in red dashed and solid line, together with the boundaries at d = 0.01: as the reader can observe, at least for this set of parameters there is just a small difference between the two boundaries, showing that the introduction of the diffusion term d does not change the stability scenario.

Zoology of the stationary and nonstationary solutions

Let us now explore, numerically, what is the typical scenario regarding the stationary and nonstationary solutions of the model in Eqs. (1.1).

It is well known that, due to the combined effect of diffraction, nonlinearity, cavity feedback and saturable absorption this kind of laser system can develop self-pulsations and modulational instability. Between the possible stationary and nonstationary solutions that can occur in the transverse plane of our laser we can list: cavity solitons, both stationary [Bache 2005] and oscillatory, chaotic solitons [Vahed 2014], and spatiotemporal chaos [Selmi 2016, Coulibaly 2017].

All of these solutions are illustrated in Fig. 1.5 in function of the pump parameter µ.

Cavity solitons

Cavity solitons (CSs), of which an example is illustrated in Fig. 1.6, are bright intensity peaks that occur in configurations where a patterned solution, due to pattern-forming instability, coexists with a stable homogeneous background. Lasers with saturable absorber, when presenting a subcritical bifurcation, represent a privileged framework for the study of CSs because it is here possible to obtain CSs against the homogeneous background given by the trivial I 0 laser-off solution as in Eq. (1.2) and illustrated in Fig. 1.3(a): therefore the contrast between the patterned solution and the homogeneous background results at its maximum, and, also experimentally, CSs emerge from a background of pure noise [Prati 2007a]. Furthermore the absence of an injected external field, presents its own experimental advantages in terms of compactness of the setup.

Physically speaking these objects arise from the double compensation of, on one side, the effects of nonlinearity and diffraction, and, on the other side, the gain and the losses: this compensation makes it possible for CSs to maintain a spe- cific shape in time and space. Due to the stability condition for I 0 , CSs are stable only below the laser threshold in this kind of system: in particular, for the same parameters as in Fig. 1.5(a), the cavity soliton branch results stable for 4.572 < µ < 5.16. In the simulations these objects are obtained starting from the laser-off solution and injecting a Gaussian beam of low intensity (of amplitude ≈ 3 and width ≈ 8 µm) for 3 ns.

Extended spatiotemporal chaos

Extended spatiotemporal chaos is a stable regime that appears both above and below the laser threshold: in particular, for the same parameters (r = 1) as in Fig. 1.5(a), its branch is stable down to µ = 4.5 below the laser threshold, and, above the laser threshold there seems to be no limit to its stability. This regime has been the subject of extensive study in [Selmi 2016, Coulibaly 2017] and it has there been characterized in the case of one spatial dimension plus time as spatiotemporal chaos through a Lyapunov coefficients analysis and a study on the Kaplan-Yorke dimension. In Fig. 1.7 we present the intensity (a) and the phase (b) of the electric field in the transverse plane for this solution. In the simulations this regime is obtained starting from values of µ above the laser threshold and giving some random initial condition to the system: given our choice of parameters, the system ends up on the only stable solution, which is the extended chaotic regime. Then to obtain the same regime also below the laser threshold we set µ < µ thr and we use as initial condition the final state of the preceding simulation, displaying spatiotemporal chaos: the system stays on that attractor allowing us to follow and plot the whole branch as in Fig. 1.5(a). The intensity of the electric field in the transverse plane continues to evolve and no periodicity has been observed either in space or time. Regarding the phase we can qualitatively observe that many full 2π phase rotations take place at any given time and, as we will see The black dashed line represent a possible spectrum profile with a power law decay that best fits the spectrum profile of the extended chaotic regime. We highlight that all these spectra have been acquired for the same value of pumping current µ = 4.8.

trate in the following Section. First of all we can observe from Fig. 1.9(a) that the amplitude of each mode interested by the extended chaotic dynamics is considerably higher with respect to the cavity soliton cases. We can also observe from Fig. 1.9(b) and (c), where all the spectrum profiles have been normalized to their maximum value, a clear enlargement of the spectrum profile in the extended chaotic case, showing a partition of energy on many different spatial scales. The black dashed line depicts a possible spectrum profile with a power law (ak -n ) decay that best fits the spectrum profile of the extended chaotic regime. It is well known in turbulence studies [Frisch 1995] that one of the requirements for a regime to be considered turbulent, in the sense of Kolmogorov, is for its spectrum to follow a power law decay. In our case we can observe that the enlargement of the spectrum in the extended spatiotemporal chaotic regime makes the system more "turbulent" than, for example, in the chaotic soliton case, nevertheless, given the comparison with the power law decaying spectrum profile in Fig. 1.9(b), we do not have any definitive arguments to claim the chaotic regime to be turbulent, especially since in our case the spectrum clearly decays faster than a power law for high modes. In any case the broadening of the spectrum confirms the chaotic nature of the studied regime [Eckmann 1985, Manneville 1995].

Chaotic cavity solitons

Chaotic cavity solitons (CCS) appear in a region of coexistence between the extended spatiotemporal chaotic regime and the stable laser-off solution (hence below the laser threshold): these are localized chaotic solutions, similar to the chaoticons described in [Verschueren 2013]. At difference with chaoticons, here the spatial size is fixed and rather small (corresponding to two or three times the diameter of a CS), whereas in [Verschueren 2013] chaoticons can exist with any size. In the simulation these objects are obtained starting from the laser-off solution and injecting a Gaussian beam of high intensity (of amplitude ≈ 5 and width ≈ 8 µm) for 3 ns: the amplitude of the injection is critical for the observation of CCS, as for lower amplitudes the system relaxes on the stable CS solution.

In the case here presented the chaotic cavity soliton branch results stable in an inset of the stable cavity soliton branch: in particular for the same parameters as in Fig. 1.5(a), the stability region amounts to 4.67 < µ < 4.84. Below the left boundary of this domain the system falls on the cavity soliton branch and above the right boundary the chaotic regime extends to the whole transverse plane. In The averaged spectrum profile depicted in Fig. 1.9 results in the end quite symmetrical, and, even if the spectrum appears larger than the one of CSs it is presents a smaller width with respect to the extended chaotic regime. In particular, in the first case depicted in Fig. 1.10(c), the radius of the ring corresponds to the maximum for modes other than zero in Fig. 1.9(c). In Fig. 1.11 we depicted the phase of the electric field at the same times as in Fig. 1.10. It is interesting to observe that in this regime optical vortices, highlighted by the white circles in Figs. 1.10 and 1.11, are clearly observed in the transverse plane: in these points the intensity of the electric field is zero or almost zero (consistently with some small numerical noise) and the phase develops a full 2π rotation around them. Since all the dynamics in the transverse plane is, in this regime, limited to a very small area, no more than two or three vortices are observed at the same time, giving us the possibility to observe with more clarity the evolution of each vortex. This observation allows us to consider the chaotic cavity soliton regime as an unstable vortex regime somehow similar to the vortex turbulence observed, in class A lasers for very low values of coherent injection [Gibson 2016]. Furthermore this leads to a possible similar interpretation of the extended spatiotemporal chaotic regime, where many optical vortices are also present, so much that it is difficult to understand their actual effect on the dynamics. Further investigations are to be addressed in some future work.

Oscillatory cavity solitons

Oscillatory cavity solitons (OCSs), like the one depicted in Fig. 1.12, are a last possible scenario. In this case the cavity soliton solution is unstable and the system develops a localized state whose peak intensity oscillates in time, as depicted in Fig. 1.13(b). This kind of behavior presents features similar to those of passive Q-switching [Dubbeldam 1999], but while passive Q-switching occurs there in a laser with a saturable absorber in a situation where no stable homogeneous solution is found, here instead this peculiar form of localized passive Q-switching occurs in a case where there is bistability with the stable background. In particular in Fig. 1.5(b), for r = 1.75, the OCS branch results stable for 4.91 < µ < 5.17. In spectrum profiles taken when the oscillation of the soliton in time reaches its minimum (blue profile) and its maximum (red profile). The Fourier spectrum at maximum oscillation closely follows the same shape of the stationary CS spectrum depicted in Fig. 1.9, except for a decrease in the amplitude of the central mode, a common thread also to the spectrum profile at minimum oscillation. These objects can be obtained starting from a stable cavity soliton solution in the middle of the stable branch, possibly not too close to the boundary of it, since the stability region for CS and OCS might differ slightly, and increasing the parameter r. For instance the simulation whose results have been depicted in Figs. 1.12 and 1.13(a-b) has been obtained starting from a stable cavity soliton solution at µ = 5 and r = 1 and increasing r up to r = 1.75. The OCS branch results also stable for a certain range of r [Vahed 2014]. Finally in Fig. 1.13(c) we illustrated the behavior of the pulsation period of the OCSs as a function of the pump parameter µ, for the same branch depicted in Fig. 1.5(b). In good agreement with the results reported in [Vahed 2014], the pulsation period of these objects appears to diverge logarithmically for low values of µ as was also observed in [Turconi 2015]. From the simulations, we observed that the diffusive term d results irrelevant when we are dealing with localized structures such as CS, CCS and OCS but it needs to be introduced in presence of STC, to avoid the formation of filaments due to high spatial frequencies excitation. The basic effect of d consists in setting a filter to the number of spatial frequencies involved in the dynamics. This phenomenon is well illustrated in Fig. 1.14, where we plotted the Fourier spectrum profile of the chaotic regime for different values of the diffusion coefficient d. Without d we observe that the Fourier spectrum flattens and the spatial structures in the intensity narrow to single independent points of the grid, due to the transportation of energy from the low unstable wave vectors to the higher ones. All these phenomena make the simulations unreliable without d when studying STC. The necessity of this term has been proven theoretically in Ref. [Fedorov 2000]. The specific value of d has been chosen phenomenologically as the smallest possible value to avoid self-collapsing without evidently disrupting or changing the dynamics of the model and the stability of its solutions.

Role of the diffusion coefficient d

Extreme event detection

The main focus of this study regards the investigation of the extended chaotic regime to search for extreme events and figure out how the occurrence of these objects might be affected by the parameters. To this aim, first of all, we have to decide what quantity gives the best measure of the amount of extreme events present in the system, and ultimately is suitable to a statistical analysis. Many different studies in optics (see for instance [Walczak 2015] and almost all the works developed in optical fibers) consider the entire matrix of intensity, i.e. all the values explored by the intensity of the electric field during the acquisition time or the simulation, as a good candidate for this kind of analysis. The main issue to be addressed when dealing with this variable is that any kind of statistics developed on it does not allow us to identify directly an extreme event. Possible deductions from a statistical analysis on the total intensity are, whether or not there might be extreme events in the system and from what value of intensity the system is most likely to exhibit one. This analysis does not discriminate between points relative to the same or to different peaks in the intensity. For this reason, if in a system we have extreme events of different durations, the statistics on the total intensity would be largely biased towards the extreme events of larger duration instead then towards the likely more disruptive short-lived extreme events. An obvious advantage of this kind of statistics is instead its simplicity.

Since in our numerical simulations it is possible to sample the transverse plane with enough accuracy to be able to follow the evolution of each individual peak we decided for a different kind of analysis where the quantity of interest is each spatiotemporal structure observed in the intensity of the electric field: we then identify as events in our analysis the local spatiotemporal maxima of the transverse field intensity. This allows us to discriminate between different peaks during the analysis and, in particular, it assigns one single value (the peak value) to each spatiotemporal structure occurring during the simulation, differently from what would occur in the total intensity analysis. Extreme events are, in this way directly identifiable by the statistics as all those spatiotemporal structures with a peak value above a certain threshold.

Maxima individuation method

The method of individuation of the spatiotemporal maxima consists in a nontrivial extension of the method introduced in Ref. [Coillet 2014] for a system of 1D+1 dimensions. If we refer to the zooms of the transverse plane intensity in Fig. 1.15, the spatial maximum highlighted by the red arrow in (a), (c) will not be selected as long as it is only a maximum in space, it will be selected just when it reaches also its maximum in time, at the point highlighted by the green arrow in (b) and only at that point it will be registered. This procedure is applied to all the local maxima in the spatial transverse plane throughout the entire simulation. The pictures in Fig. 1.15 are separated 9 ps from each other just to give a good idea of how the selection method works. Nevertheless we want to remark that the sampling of the transverse plane occurs every 1 ps (and the integration time unit is 100 fs). The points highlighted by the red circles in Fig. 1.15 are other spatiotemporal maxima occurring at the time of the snapshot.

From a more numerical point of view this procedure is accomplished as follows: each point of the 3D spatiotemporal grid is compared with its 26 neighbors, including those on the diagonal directions, as depicted in Fig. 1.16, somehow resembling the concept of looking at the spatial and temporal derivatives of the intensity. The basic requirement for each point of the grid to be considered a local maximum if its value of intensity is the highest among all its neighbors, underestimates the number of local maxima due to the numerical discretization, since there are many cases in which a peak presents two points with the same intensity value. The requirement we make is then for a point to be a possible local maximum if its intensity value is higher or equal to the values of intensity of all its neighbors and strictly higher than that of a neighbor in at least one of the directions depicted in Fig. 1.16. With this definition we are able to avoid the selection of points with values of intensity that lie on a plateau but we still have to re-move the points at the border of intensity plateaux where the intensity starts to decrease in one of the directions. We highlight that there is no real possible distinction between these border plateau points and the case of two or three points on a peak event having the same value, so the strategy consists in removing all the border plateau points except one, so that the statistics is not biased by the possibly large amount of points lying at the edge of same-valued regions in intensity and, at the same time we do not loose the peak points that represent the events we are looking for. To accomplish this goal for each candidate local maximum we select all the other local maxima, if any, that are closer than five units (note that some points will be obviously selected more than one time but this will not matter). Then, considering just the points with the problem of having other possible local maxima in their proximity, we divide them in groups of all the points that are closer than thirty units in any direction of the grid. This is done so that even if the border plateau points cover a big area, all points have a good chance to be identified at once. The distance of thirty units has also to be chosen as not too large because it would risk considering other areas that are not related to the same plateau. Finally for each group we select all the points that have the same intensity values and we erase them all except one.

The remaining points are the spatiotemporal maxima on which the statistical studies have been developed. A threshold of detection at I = 0.5 has also been introduced to avoid numerical noise to bias the statistics.

Finally since for each simulation we are usually dealing with a big amount of data (matrices of ≈ 256 × 256 × 25000 units 3 ) we preferred to divide the data and detect the local maxima in subset grids of ≈ 256 × 256 × 200 units 3 . Doing so, to avoid missing any local maximum on the border of each grid, every time we select a new grid we shift back the time coordinate of 10 units. Some local maxima might be selected two times but after the full procedure only one of the two points with the same coordinates will be kept into account.

Statistical analysis

The statistical analysis has been developed on all the spatiotemporal maxima observed during 25 ns long simulation, with an acquisition rate of one image of the transverse plane every 1 ps. The integration time step, as previously explained, is 100 fs. The chosen spatial size is 256 × 256 pixels which roughly corresponds to 256 × 256 µm 2 , being the spatial step 0.25 units and space rescaled on the diffraction length (which is about 4 µm).

PDF of the spatiotemporal maxima

Three typical probability density functions (PDFs) of the spatiotemporal maxima are illustrated in Fig. 1.17. To study the presence of extreme events we introduced three threshold definitions. The first two threshold definitions are the most commonly used in rogue wave (RW) analysis. Fig. 1.17 -PDFs of the spatiotemporal maxima detected numerically during the simulations, as described in the previous Section. (a) corresponds to a simulation for µ = 7 and r = 2.5, (b) corresponds to a simulation for µ = 5 and r = 2.4 and (c) corresponds to a simulation for µ = 5 and r = 1. All simulations lasted 25 ns. The three vertical dashed lines correspond to the three different thresholds described in the text. The green and magenta lines are, respectively a Gumbel and a Weibull PDF computed through the average and standard deviation of the data. [START_REF] Oppo | Self-organization, pattern formation, cavity solitons, and rogue waves in singly resonant optical parametric oscillators[END_REF], Liu 2015, Gibson 2016], corresponding to the vertical black dashed line in Fig. 1.17, consists in the mean intensity, averaged on all the points of the transverse plane during the whole acquisition time, plus eight times the standard deviation. The second threshold definition, corresponding to the vertical red dashed line in Fig. 1.17, traditionally used in oceanography, consists in two times the significant wave height H s , which is the mean of the highest third of the spatiotemporal maxima intensity values. We point out that through this threshold definition it is possible to avoid spatially large events to count more in the statistics because their width covers more grid points. The third threshold definition, corresponding to the vertical blue dashed line in Fig. 1.17, has been introduced specifically in this analysis for self-consistency with the data and it somehow resembles the first definition for its form: it consists in the average of the intensity values of just the spatiotemporal maxima plus eight times their standard deviation. As noticeable in Fig. 1.17 this third definition is by far the most restrictive, which remains true for all the set of parameters considered in our analysis. For this reason, in this Chapter, we will consider the third threshold as the main indicator of the presence of extreme events in the data. We would like also to highlight that all the RW thresholds here introduced have been computed on the data above I = 0.5, to avoid any numerical noise to bias the statistics, as it was done for the maxima detection method: this approach makes the thresholds even more restrictive and is commonly used also in the experiments [Selmi 2016]. In Fig. 1.17 we also introduced two different probability density functions, the first one being the Gumbel PDF, plotted in green

Gum PDF (z) = 1 β exp[-z -exp(-z)] with z = I max -I max β + γ ,
with γ Euler's constant, and the second one being the Weibull PDF, plotted in magenta

Wei PDF (I max ) = k λ I max λ k-1 exp - I max λ k . (1.15)
Both of these probability density functions and the reason for their introduction will be described more into details in the next Subsection. For now we just want to remark that these two PDFs are not a fit of the our data: the parameters β, k and λ are instead computed from the data average, median and standard deviation values and they represent what behavior the data would follow if the statistics predicted by these distributions was the one to which the system obeyed. We also observe that, for this reason, the probability density functions here shown have no degree of freedom, hence the fact that the data well follows one or another behavior is a stronger statement than what would have been in case of a simple fit. In Fig. 1.18 we also plotted, for comparison, the PDF of the total intensity, that is, all the values in intensity of all the points in the spatial grid at all the recorded times, which is actually the standard quantity studied by other groups. The blue solid line corresponds to a negative exponential PDF computed as exp(-I/ I )/ I , once again this PDF is not a fit, it is instead computed on the mean intensity value of the data. The three plots in Fig. 1.18 have been shown for comparison with other works on extreme events but we would like to point out once again the few issues that this representation of the data can lead us into: looking at the PDF of the total intensity we can get an idea of whether the data under analysis is likely to exhibit extreme events given a positive deviation from the negative exponential decay rate. Nevertheless this data analysis treats equally all the intensity values assumed by all the points of the grid, regardless of the presence of spatiotemporal structures within the grid: this means that any spatiotemporal structure that is somehow larger in space or longer in time is forced to weight more, no matter its intensity, in the data analysis than smaller, and potentially more destructive, structures. One final observation about the two different kinds of PDFs here shown is the following. For the PDF of the spatiotemporal maxima we are looking at the distribution of events that could be extreme for certain values of intensity: any given point of the plot above the threshold directly corresponds to a certain known number of extreme events; in the case of the PDF of the total intensity we are instead not able to associate directly the points that deviate from the negative exponential to any event.

From this data analysis we can conclude that for certain parameter values Figs. (see 1.17(b)-1.18(b)) the system clearly exhibits extreme events, with a large number of events above the third threshold in Fig. 1.17(b) and a clear deviation from the negative exponential PDF In Fig. 1.18(b). In Fig. 1.17-1.18 we used specific probability density functions for comparison with the data. In the following we are going to explain the reason for this choice.

Rayleigh and negative exponential PDF

The use of a negative exponential PDF in comparison with the data of total intensity is a known procedure in extreme event analysis, especially in fiber optics [Suret 2016, Randoux 2016b] and comes directly from the analogy with oceanography. In this Subsection we will summarize the reason for its application. Let us consider the sea surface elevation variable η(t) [Onorato 2013, El Koussaifi 2017] and assume that such variable follows a Gaussian distribution such that

P g (η(t)) = 1 σ √ 2π e -η 2 (t) 2σ 2 .
Considering the complex signal z(t) = η(t) + i η(t), with both η(t) and η(t) following a distribution given by P g (η(t)), the envelope of the wave is, in good approximation, given by the amplitude of z(t) as

ρ(t) = |z(t)| = η 2 (t) + η2 (t) .
The probability distribution of the wave envelope is then computed as follows: let us first consider the joint probability of η(t) and η(t) (uncorrelated by construction),

P(η(t), η(t)) = P g (η(t))P g ( η(t)) = 1 2πσ 2 e -η 2 (t)+ η2 (t) 2σ 2 = 1 2πσ 2 e -ρ 2 (t) 2σ 2 .
Changing variables from (η(t), η(t)) to amplitude and phase (ρ(t),θ(t)), the transformation is given by the Jacobian

J = cos θ(t) -ρ(t) sin θ(t) sin θ(t) ρ(t) cos θ(t)
and its determinant is det (J) = ρ(t). Hence

P(η, η)dηd η = ρ(t) 2πσ 2 e -ρ 2 (t) 2σ 2 dρdθ = ρ(t) σ 2 e -ρ 2 (t) 2σ 2 dρ (1.16) = P r (ρ(t))dρ
with P r (ρ(t)), the distribution followed by the amplitude ρ(t) and the wave envelope, being a Rayleigh distribution:

P r (ρ(t)) = ρ(t) σ 2 e -ρ 2 (t) 2σ 2 .
Finally, in oceanography, the distribution of heights of the wave envelope (from minimum to maximum) results more relevant for statistical and practical purposes than the amplitude distribution. In good approximation, for a narrow band process [Onorato 2013], we can assume here that any excursion from trough to crest in the surface elevation is close to ρ(t) = H/2 where H is the height. Hence the height distribution follows a Rayleigh distribution as well, unless a factor 1/2

P h (H) = H 2σ 2 e -H 2 8σ 2 .
(1.17) Now, in optics, the easier variable to record experimentally is usually the intensity of the electric field: the distribution followed by the intensity I(t) = ρ(t) 2 is then given by changing the integration variable in Eq. (1.16)

ρ(t) σ 2 e -ρ 2 (t) 2σ 2 dρ = 1 2σ 2 e -I(t)
2σ 2 dI = P e (I(t))dI

with P e (I(t)), the distribution followed by I(t), being a negative exponential distribution

P e (I(t)) = 1 2σ 2 e -I(t) 2σ 2 .
(1.18)

Finally, for data treatment, we can observe that the mean value of the distribution P e (I(t)) is 2σ 2 and then:

P e (I(t)) = 1 I(t) e -I(t) I(t) . (1.19)
As a side note, from Eq. (1.17), we can observe also the origin of the (first) extreme event threshold definition containing 8σ, used earlier in the text: the term in the exponent for P h (H) can be considered as two times the squared ratio between the wave height and a certain height threshold

- H 2 8σ 2 = -2 H 2 16σ 2 = -2 H 2 H 2
thr with H thr = 4σ. Similarly to the case of the significant wave height H s , oceanographers consider a wave as rogue if its height is larger than 2H thr = 8σ, which is the same threshold definition we represented with the dashed vertical black line in Fig. 1.17 (here for the intensity instead than for the amplitude), when considering that the average of the sea surface elevation is zero, given the assumption of Gaussianity.

Weibull PDF

The Weibull PDF is often used in oceanography [Forristall 2000, Bitner -Gregersen 2012, Toffoli 2016] and in optics [Dudley 2008, Onorato 2013, Gibson 2016].

The reason for such an extensive use lies in the versatility of such density function

Wei PDF (x) = k λ x k-1 e -x k ,
where, with respect to Eq. (1.15) we just set x = I max /λ, in fact, we can observe that for

k = 1 Wei PDF (x) = 1 λ e -x is a negative exponential distribution k = 2 Wei PDF (x) = 2 λ xe -x 2 is a Rayleigh distribution .
As explained in the previous Section, in oceanography, under the assumption of a Gaussian distribution for the sea elevation, the distribution of the water wave envelope heights should in good approximation follow a Rayleigh distribution: the Weibull distribution for k > 0 gives then an idea of whether and how much the data deviates from a Gaussian-related behavior, leading to possible rogue waves in the system.

In optics instead we can observe that any Weibull distribution, for k > 1, well describing the data, indicates a deviation from the negative exponential statistics with the parameter k quantifying the magnitude of such deviation.

In our specific case we chose to use the Weibull PDF as a comparison for the PDF of the spatiotemporal maxima and not in the PDF of the total intensity, because in the second case it did not seem to well describe the data for any parameter choice we made, thus not giving us any insight of how much the data deviated from the negative exponential. Instead the application to the PDF of the spatiotemporal maxima seemed more interesting because the Weibull PDF was able to well describe the data in the case of few extreme events in the system (see for example Fig. 1.17(a) and (c) where, respectively k ≈ 1.52, 1.54). In the case of many rogue waves, depicted in Fig. 1.17(b), the deviation from the negative exponential tail seems too big even for the Weibull distribution to adapt.

Gumbel PDF

The use of the Gumbel PDF as a comparison for the spatiotemporal maxima data requires a larger explanation, especially for its connection with classical extreme value theory [Gumbel 1958[START_REF] Coles | An Introduction to Statistical Modeling of Extreme Values[END_REF], Albeverio 2006].

The main goal of classical extreme value theory, given a set of random, independent and identically distributed (i.i.d.) variables, is the study of the extremes (both maxima and minima) of the time series provided by such variables. One of the main results of this theory consists in the following theorem [Albeverio 2006]:

Theorem (Fisher-Tippett-Gnedenko). Let {X i } i∈N be a sequence of i.i.d. real random variables, distributed according to some cumulative distribution F, and consider the partial maxima M n := max{X i : i = 1, ..., n} with n ∈ N.

If a sequence of pairs of real numbers {a n , b n } n∈N exists such that a n > 0 and

M n -b n a n d -→ H for n→ ∞, (1.20)
where d -→ indicates convergence in distribution, then H belongs to one of three nondegenerate limit distribution, namely Fréchet, Weibull and Gumbel distributions.

Less rigorously what the theorem states [Albeverio 2006] is that if the maximum M n has a limit distribution then such distribution belongs to one of three standard types.

Considering now [Albeverio 2006] the probability P that M n -b n a n will assume a value less or equal to x, called also cumulative probability, we can observe that

P M n -b n a n ≤ x = P{X i ≤ a n x + b n , ∀i = 1, ..., n} = F(a n x + b n ) n
where in the first equality we have observed that asking that for the n maxima of the sequence {X i } i∈N to be ≤ x is equivalent to ask for all the terms of the sequence to be ≤ x. In the second equality we observed that, since the variables of the sequence X i are F-distributed, then the probability that one variable of the sequence X i for a certain i is less or equal to a n x + b n is F(a n x + b n ). Furthermore the probability that every variable of the sequence is less than a n x + b n is F(a n x + b n ) n because the variables are assumed independent. At this point Eq. (1.20) can be expressed as follows:

lim n→∞ F(a n x + b n ) n = H(x) x ∈ R
with H(x) the cumulative distribution functions (CDF) of one of these three standard extreme value distributions

Fréchet : Φ α (x) = 0 x ≤ 0 e -x -α x > 0 Weibull : Ψ α (x) = e -(-x) α x ≤ 0 1 x > 0 Gumbel : Λ(x) = e -e -x x ∈ R with α > 0.
It is important to notice that the Weibull extreme value distribution here casted is more properly called reversed Weibull distribution, as in fact the CDF here represented corresponds to the complementary cumulative distribution function.

As it turns out in fact while the Gumbel and the Fréchet CDF relate to maxima, the ordinary Weibull CDF

Wei CDF = 1 -Ψ α (x) = 1 -e -(-x) α x ≥ 0 0 x < 0
would relate to minima, as the smallest extreme value instead of the greatest.

In order then to make a comparison with the other two standard extreme value distribution it is needed to reverse Wei CDF , so that also this distribution would relate to maxima.

Up to a linear transformation the Fréchet, Weibull and Gumbel distributions can be rearranged in the following common form, also known as generalized extreme value (GEV) distribution [Albeverio 2006]

P γ (x) = exp -1 + γ x -m σ x -1 γ + for x ∈ R
where ( * ) + = max(0, * ) and γ is the extremum value index. For γ > 0 the GEV distribution corresponds to the Fréchet distribution Φ α (x), for γ < 0 it corresponds to the reversed Weibull distribution Ψ α (x) and for γ → 0 it corresponds to the Gumbel distribution Λ(x).

We would like to remark [Albeverio 2006] that there are many distributions that do not belong to the three mentioned domains of attraction: for instance the maximum distribution for geometric and Poisson distributions cannot be well approximated by any of the standard extreme value distributions, instead Pareto-like distributions usually fall in one of the three given categories.

Another remark on the extreme value theorem [Albeverio 2006] is that its validity is for statistically independent data, but in reality the data under analysis is often statistically dependent. Even if this is beyond the purpose of this digression and this Thesis, in this case we need to wonder to what extent the approximations of the GEV distribution are still valid for dependent data. The problem is usually solved introducing appropriate mixing conditions, for instance, given a stationary time series X t with t ∈ N, with a marginal distribution function F, we can still state that

P(X t ≤ a n x + b n , ∀i = 1, ...n) = F θn (a n x + b n )
where θ is called the extremal index. If θ ≈ 1 then it means that the data considered is almost independent. Some examples for the use of the Gumbel distribution in oceanography, as directly related to classical extreme value theory are [Bitner -Gregersen 2012],

where the assumption of wave crests as independent events is made, and [Fedele 2012a, Fedele 2012b]. In general we can observe that for a sufficient number of data the wave crest distributions in these papers seem indeed to follow a Gumbel CDF.

In our present case, we apply, for comparison, the Gumbel PDF as an extreme value PDF that can be considered to analyze data with large deviations from their average. The other two distributions indicated in the extreme value theorem do not seem to well describe the data, for any choice of parameter made. To give an idea about a direct application of the theorem, we can assume that for each spatiotemporal maximum it is possible to find a time series or a set of grid points of which such spatiotemporal maximum presents the largest value of intensity. These sets of grid points are assumed to be independent from each other and together they could play the role of the sequence of i.i.d. real random variables (each variable able to assume all the values in its set of grid points). Then if the maxima of each variable, that are our spatiotemporal maxima, have a limit distribution, that has to belong to one of the three standard extreme value distributions. As we will see further into details in the following Section, for the most extreme choices of parameters (see for example Figs. 1.17(b) and 1.18(b)), the Gumbel distribution does not well describe the data, which implies that the theorem does not hold in this case. This occurs because either there is no limit distribution or, most probably, the sets of grid points is not to be approximated anymore as independent: in fact the data shows the most extreme behavior when, globally, a self-pulsating, Q-switching-like phenomenon takes place. As observed in the introduction of this Thesis, in the general theory of extremes the data of extreme phenomena is often found to be dependent.

Dependence on the laser parameters

One of the aspects of this data analysis we are particularly interested in is identifying what parameter region is more likely to exhibit extreme events, also from an experimental point of view. The main parameters that can affect the dynamics of our model in Eqs. (1.1) are µ, the pump parameter and r, the ratio of the carrier lifetimes. To compare the extreme nature of the intensity values in each simulation we considered two indicators. The first indicator is the kurtosis of the PDFs, both of the spatiotemporal maxima and of the total intensity. The kurtosis is the ratio of the forth moment about the mean to the square of the variance of the data. The second indicator, applied mainly to the PDFs of the spatiotemporal maxima, is the rogue wave ratio, which is the ratio of extreme events, according to one of the three rogue wave threshold definitions, to the number of all the spatiotemporal maxima. In the case of the PDF of the total intensity this indicator simply consists in the ratio of the number of grid points with an intensity value above the first threshold to the total number of grid points on which the simulation has been computed. In Fig. 1.19 we depicted four density plots relative to the PDF of the spatiotem- The first two RW thresholds, which do not introduce a particularly restrictive condition on the data, seem to show a similar behavior in the density plots (a) and (b), identifying, as the region with the higher ratio of extreme events, the one located to the higher left corner of the map. We would like to stress out that the grey region at the top left corner is not to be considered part of the analysis, since extended spatiotemporal chaos is not stable in this area: there is only the laser-off solution. The third RW threshold, more restrictive for our data, and the kurtosis density plots For the sake of comparison, in Fig. 1.20 we depicted the density plots of the total intensity PDFs for different values of µ and r, relative to the (first threshold) RW ratio (a) and the kurtosis (b), here illustrated after subtracting the value of kurtosis for the negative exponential PDF K exp = 9. A comparison with the results depicted in Fig. 1.19 shows that all the data plots seem visually similar, with a typical structure at the left boundary of the extended spatiotemporal chaos regime. Some small differences are in any case specific to the different quantity under study: in particular We believe that this small difference between the density plots of the total intensity and the spatiotemporal maxima is to be related to the global behavior of the system, in particular the time trace of the spatially averaged intensity: in the peak parameter region highlighted in Fig. 1.19(c) and (d) such time trace presents a pulsating behavior, as depicted in Fig. 1.21, somehow reminiscent of the passive Q-switching obtainable in a purely temporal LSA system. This pulsating behavior implies that the transverse plane presents, most of the time, low values of intensity, with some, but not many, spatiotemporal maxima, and, every once in a while, it bursts everywhere to high values of intensity with a big number of spatiotemporal maxima. Hence the long time intervals where low values of intensity are displayed, will count more in the statistics of total intensity than in the statistics of the spatiotemporal maxima. For this reason, according to the statistics of the total intensity, and to the threshold computed on such data, the system seems to present a lower number of extreme events than it actually contains. This Q-switching-like behavior of the average intensity then stops for higher values of r and the bias on the total intensity extreme event estimation method stops as well. To show better the kurtosis and RW ratio behavior of Fig. 1.19 and their agreement, we also plotted in Fig. 1.22 the values of the two indicators for fixed r = 2.2 (a) and µ = 5 (b). In particular in Fig. 1.22(a) we preferred to display a plot for a lower value of r, instead than for the most favorable situation of or in the bistable region, it will be more likely to observe extreme events also experimentally.

Comparison with other works in the field

We would like now to compare our results with those of Refs. [Selmi 2016, Coulibaly 2017], concerning the same system. In Ref. [Selmi 2016] the authors observe extreme events experimentally in the setup of a VCSEL with an intracavity saturable absorber. The data experimentally collected consists in the intensity time trace, averaged on a small area of diameter ≈ 5 µm [Coulibaly 2017], due to the size of a pinhole set before acquisition and the finite size of the detector, a fast avalanche photodiode. The model is numerically studied [Selmi 2016, Coulibaly 2017] in 1D+time, without radiative recombination (B = 0) and with a linewidth enhancement factor for the absorber equal to zero (β = 0) as in Ref. [Bache 2005]. In order to take into account the finite linewidth of the gain, the authors prefer to consider a tophat finite pump profile, as the one in Ref. [Prati 2010a]: a spatial confinement in the transverse plane introduces in fact a limitation in the number of transverse modes involved in the dynamics.

The authors in [Selmi 2016, Coulibaly 2017] characterize the extended spatiotemporal chaos in 1D+time through Lyapunov coefficients analysis and the Kaplan-Yorke dimension, and relate the presence of extreme events to an increase of the chaoticity of the system [Selmi 2016] and to quasiperiodic extended spatiotemporal intermittency [Coulibaly 2017]. We would like to highlight that one of the main differences between the numerical results in Refs. [Selmi 2016, Coulibaly 2017] and the present study consists in the fact the authors have limited themselves to 1D spatial dimension instead of considering the full 2D transverse plane: this is justified in order to describe the experimental results where a gold mask, set on the VCSEL surface, assigns to the pump a rectangular geometry. Furthermore all the statistical analysis in [Selmi 2016] is developed on the (spatially averaged) total intensity time trace. In the subsequent paper [Coulibaly 2017] the theoretical statistical analysis has been developed on the spatiotemporal heights in the 1D+time system.

We computed a set of 25 ns-long simulations with the same parameters as in [Bache 2005], referenced in [Selmi 2016, Coulibaly 2017]: that is b = 0.005, r = 1, s = 10, B = 0, α = 2, β = 0, γ = 0.5 and variable µ. In Fig. 1.23 we depicted the homogeneous stationary solution (black line) together with the timeaveraged maximum intensity of the extended spatiotemporal chaotic branch, plotted in orange. It is noticeable that for this set of parameters the chaotic branch is not stable below the laser threshold µ thr = 1.5 (its left boundary lies at µ ≈ 1.49). This implies that there is almost no multistability between the chaotic regime and the laser-off state.

Applying our method of analysis on the spatiotemporal maxima to this set of simulations, we observed a clear increase of both the indicators of kurtosis and (third threshold) RW ratio of the spatiotemporal maxima PDFs as depicted in Fig. 1.24. This result is in agreement with what we have already observed in the previous Section and it confirms that the system is more likely to exhibit extreme events for low values of µ, possibly below the laser threshold, or, if that is not allowed by the stability of the chaotic branch, as close as possible to µ thr . This result somehow disagrees with what has been observed in [Selmi 2016], where the authors report an increase of the extreme events of the total intensity time trace for higher values of µ and related this phenomenon to the similar increase of the Kaplan-Yorke dimension (see for instance Fig. 1.25). We have already noticed that, in the regions where our RW indicators are the largest, the spatially averaged intensity might actually well behave in time, even showing a Q-switching-like behavior, as in Fig. 1.21 where no extreme events could be found: the discrepancy between the results in Ref. [Selmi 2016] and this Thesis is probably due to the quantity under study, somehow reductive of the nature of the system, and also to 1D geometry. In fact in [Coulibaly 2017], an almost opposite behavior for the kurtosis and RW ratio, in better agreement with our results, is observed.

Finally we compared into details our usual method of analysis with the method used in the experiments, which consists in studying the intensity time trace, averaged in a small area of the transverse plane. A typical time trace of this data is depicted in Fig. 1.26 for µ = 1.49. In Fig. 1.27 we depicted the PDF of the spatiotemporal maxima obtained from the analysis of the data for (a) µ = 1.49 (b) µ = 3, which are respectively the left and right boundaries of the chaotic branch in Fig. 1.23. Furthermore in Fig. 1.27(c) we plotted the PDF of the temporal maxima observed in the time trace (spatially averaged on a small area) as depicted in Fig. 1.26 for µ = 1.49, and the same PDF for µ = 3. The RW thresholds have been computed for (c) and (d) on the respective intensity averaged time traces. A small difference with the method used in [Selmi 2016] consists in the fact that we are evaluating the statistics of the maxima instead of that of the heights, but given that the baseline in Fig. 1.26 is set to zero, the discrepancy between the data is assumed to be small. It is clear from Fig. 1.27(a) and (c) that the experimental method under- estimates the presence of extreme events in the system in their most favorable regime, presenting (c) just one point above the third RW threshold in comparison with the large amount of events with high values of intensity observed in the spatiotemporal maxima PDF (a).

Waiting time statistics

The statistics of the waiting times between consecutive extreme events has been reported to follow a log-Poissonian statistics in [Arecchi 2011, Dal Bosco 2013].

To verify this observation in our system we considered a 250 ns-long simulation for µ = 6 and r = 1: this is not one of the most favorable regimes for the occurrence of extreme events but, even if the number of events above the third threshold is small, given the longer simulation time, we can any way obtain a good amount of extreme events for our study. For each extreme event registered at a time t w in the position (x w , y w ) of the transverse plane we verified if any other extreme events occurs at a further time Fig. 1.28 -PDF of the logarithmic waiting times between two consecutive extreme events in the same 5 pixel area for a simulation with µ = 6 and r = 1. In green we presented a fit of the data with a negative exponential.

t w + t in an area of 5 pixels centered on the same position (x w , y w ). We would like to highlight that, due to this procedure, a longer simulation time is also needed in order to have enough cases where an extreme event occurs at least a second time in the same spatial location. In Fig. 1.28 the resulting PDF of log[t w + t/t w ] seems to be in good agreement with a negative exponential fit, plotted in green, leading us to conclude that the extreme events in our system are, in good approximation, independent variables as in Ref. [Arecchi 2011, Dal Bosco 2013]. This result is in contrast with what has been observed in [Coulibaly 2017], where the waiting time distribution seems to deviate from a Poissonian statistics, following a stretched negative exponential. Once again the discrepancy between the two studies seems to lie in the different quantity under analysis. In Appendix D we summarize the derivation of the log-Poissonian distribution performed in [Sibani 1993, Sibani 2003].

Profiles and cavity soliton comparison

We are now going to analyze the spatial and temporal profiles of the extreme events detected in a favorable parameter region.

In Fig. 1.29(a) we depicted an example of the kind of extreme events occurring in the intensity in the transverse plane, here relative to a simulation for µ = 4.8 and r = 2.2. Furthermore in Fig. 1.29 we reported the spatiotemporal maxima PDF (b) and the total intensity PDF (c) for this specific simulation.

In Fig. 1.30 we illustrated the corresponding profiles in time (a) and space (b) for the intensity I (red solid line) and the carrier densities in the active D (green dotted line) and passive D (yellow dash-dotted line) media. For the sake of simplicity we limited our data analysis to the x and y axes. We can observe that the occurrence of a high peak in the intensity implies, as expected, the presence of a dip in the carrier density of the active medium and a peak in the carrier density (c) Total intensity PDF. of the passive medium. In particular in the temporal profile of Fig. 1.30(a), acquired registering the intensity values throughout the simulation in the spatial point where the extreme event occurs, we can also appreciate the slower dynamics of the carrier density in the recovery after the extreme event occurrence.

We would like to remark that, due to the movement of the spatiotemporal struc-tures along the transverse plane in time, in some cases it is difficult to acquire an extreme event profile where only one structure is involved: what often occurs is that the profile along one of the axis contains the merging between two or more structures thus overestimating the size of a single extreme event. For this reason we limited our analysis to the extreme event profiles whose closest minima (to the left and to the right of the peak) have intensity values below I = 1. The typical FWHM in space for the extreme events observed is 6 µm and the typical FHWM in time is 16 ps. These values remain almost constant across all the different parameter choices made throughout this study, suggesting the presence of a typical spatial and temporal size for this kind of phenomena. 1.5(a) and here studied, has been simulated for r = 1 and does not coexist with OCSs [Vahed 2014]. We would like to point out that the parameter choice for the comparison with stationary CSs is not particularly favorable to extreme events (see Fig. 1.17(c) and 1.18(c)), nevertheless some extreme events still occur in the system allowing for this shape comparison to take place. Extreme events and cavity solitons, both stationary and self-pulsing, show a very similar spatial and temporal size, seemingly implying the same generating mechanisms to lie at the formation of both phenomena. In particular the generating mechanism of these structures in space is related to the modulational instability of the homogeneous solution and their scale is determined by the diffraction length (which depends on the cavity length and the wavelength of light) [Spinelli 1998]. Instead the formation of these structures in time is due to the effect of the Hopf instability on the stationary solutions (both the homogeneous and the CS ones), leading to Q-switching [Dubbeldam 1999].

Finally we observe that, for values of r where CSs exist, extreme events are more frequent when the system presents multistability with soliton solutions. Extreme event occurrence may therefore be enhanced by the presence of cavity solitons and related to the existence of the dissipative soliton attractor.

Future perspectives

This work presents two main future perspectives: the first one consists in the investigation of predictability both in space and time, for instance comparing the shape of the field intensity when approaching an extreme event to check for possible regularities (see for example [Zamora-Munt 2013, Birkholz 2015[START_REF] Alvarez | [END_REF]). This study requires a specific nontrivial technique since the structures observed in the chaotic regime move in time along the transverse plane and to develop a rigorous study it is necessary to be able to identify and "follow" a structure also when it does not present any peak. This development would be in any case very interesting, also in light of the hydrodynamical analogy and, more in general, of the possible analogies with other extreme event fields, since the presence of a typical shape either in space or time could be identified as a precursor of the extreme event itself, thus enabling us to reduce its possible damages.

Another possible perspective related to the study of this system is the investigation of its conservative limit, which we are going to explain in more detailed in the following Subsection and in Appendix E.

Conservative limit of a LSA

It is possible to obtain a conservative limit for the model described in Eqs. (1.1). In particular this system can be described as a particle of unitary mass subject to a Toda potential. This alternative description, in the purely temporal case, has already been proven for a Class B [Arecchi 1984] laser in [START_REF] Oppo | Toda potential in laser equations[END_REF]] and verified experimentally in [Cialdi 2013] but it has been extended only recently to the case of a passive medium inside the laser cavity [Rimoldi 2014].

Let us consider first Eqs. (1.1) without the radiative recombination terms B = 0 and neglecting the transverse Laplacian ∇ 2 ⊥ : the equations that describe a purely temporal semiconductor laser with an intracavity saturable absorber are the following

Ḟ = [(1 -iα)D + (1 -iβ) D -1] F , (1.21a) Ḋ = b µ -D 1 + |F| 2 , (1.21b) Ḋ = rb -γ -D 1 + s|F| 2 .
(1.21c)

We can recast Eq. (1.21a) for the electric field F (and the corresponding equation for its complex conjugate) as two equations, one for intensity I and the second for the phase φ such that F = √ Ie iφ . We then obtain

İ = 2I(D + D -1) , φ = -(αD + β D) , Ḋ = b [µ -D (1 + I)] , Ḋ = rb [-γ -D (1 + sI)] .
The trivial stationary solution is given by

I s,0 = 0 , D s,0 = µ , Ds,0 = -γ .
The nontrivial stationary solution is

φ s = -ω s t , D s = µ 1 + I s , Ds = - γ 1 + sI s ,
with I s such that D s + Ds -1 = 0. We also define ω s the reference frequency of the emitted field

ω s = αD s + β Ds = α + γ α -β 1 + sI s ,
which implies [Bache 2005] that the laser frequency changes with the field intensity, provided that α = β. In any case the stability of the nontrivial solution is not influenced by the value of ω s .

Let us now consider just the equations for I, D and D

İ = 2I(D + D -1) , Ḋ = b [µ -D (1 + I)] , Ḋ = rb [-γ -D (1 + sI)] .
Introducing the following variables

x = log I I s , D = D s (1 + n) , D = Ds (1 + n) , we obtain ẋ = 2 (D s n + Ds n) , (1.24a) ṅ = -b [I s (e x -1) + n(1 + I s e x )] , (1.24b) ṅ = -rb [sI s (e x -1) + n(1 + sI s e x )] .
(1.24c)

If we now compute ẍ, we obtain the following expression

ẍ + 2b (D s + s Ds r) I s (e x -1) (1.25) = -2b (D s n + Ds r n) -2b (D s n + Ds rs n) I s e x , which for r = 1 gets reduced to ẍ + 2b(µ -γ -1)I s (e x -1) = -b ẋ -2b(D s n + Ds s n)I s e x ,
where we have observed that

D s -Ds -1 = 0 µ 1 + I s -I s 1 + I s -γ 1 + sI s -sI s 1 + sI s -1 = 0 µ 1 - I s 1 + I s -γ 1 - sI s 1 + sI s -1 = 0 µ -γ -1 = I s µ 1 + I s - γs 1 + sI s µ -γ -1 = I s (D s + s Ds ) .
Hence we can deduce that for similar values of the relaxation times for D and D the first term on the right hand side in Eq. (1.25) becomes purely dissipative above the laser threshold and the second one acts as a forcing or damping term, depending on the values of the parameters involved. The conservative limit is then given by ẍ + 2b (D s + s Ds r) I s (e x -1) = 0 (1.26)

ẍ + dV LSA (x) dx = 0
where we introduced the equation that describes the motion of a unitary mass oscillator.

If we identify ω 2 LSA = 2b (D s + s Ds r) I s as the frequency of the relaxation oscillations on the stable branch to the right of the Hopf bifurcation, we can obtain an expression for the potential

V LSA V LSA (x) = ω 2 LSA V(x) with V(x) = e x -x where V(x) = e x -
x is actually the Toda potential as studied in [START_REF] Oppo | Toda potential in laser equations[END_REF]] and illustrated in Fig. 1.32. Eqs. (1.24) in the conservative limit can be recast as

ẋ = 2 (D s n + Ds n) , (1.27a) ṅ = -b [I s (e x -1)] , (1.27b) ṅ = -rb [sI s (e x -1)] ,
(1.27c) which in the original variables corresponds to It is possible to obtain an analytic approximation for x, n and n the variation of energy along a period in the limit of large energies. In this limit it is also possible to give an analytically expression for the period and width of the pulses, as defined in [Cialdi 2013], which results in very good agreement with the results of numerical simulations both in the Q-switching and the damped oscillation regimes. All these results have been illustrated into details in Appendix E.

İ = 2 [(D -D s ) -( D -Ds )] I , Ḋ = -bD s (I -I s ) , Ḋ = -rb Ds s (I -I s ) .
The reintroduction of the radiative recombination term B = 0 does not affect the conservative limit since B appears only in terms ∝ n, n, n 2 , n2 . In fact the complete Eqs. where the term ∝ n 2 , n2 are neglectable. Computing ẍ we obtain

ẍ + 2b (D s + s Ds r) I s (e x -1) = -2b (D s n + Ds r n) + 2B D 2 s n + D2 s r n -2b (D s n + Ds rs n) I s e x
and the conservative limit is still the same as in Eq. (1.26).

Reintroducing the equation for the phase, which does not affect the conservative nature of the system, and recasting the equations for F, D and D we obtain

Ḟ = F [D(1 -iα) -D s + D(1 -iβ) -Ds ] Ḋ = -bD s (I -I s ) Ḋ = -rb Ds s (I -I s ) ,
and, observing that D s + Ds -1 = 0, this gets reduced to

Ḟ = F [D(1 -iα) + D(1 -iβ) -1] , (1.29a) Ḋ = -bD s (I -I s ) , (1.29b) Ḋ = -rb Ds s (I -I s ) , (1.29c)
where we can notice that the conservative limit affects only the equations for the carrier densities.

Since the introduction of spatial coupling through the transverse Laplacian ∇ 2 ⊥ term is also not supposed to affect the conservative nature of the system we can now rewrite Eqs. (1.29) as the the conservative limit of a broad-area semiconductor laser with saturable absorber

Ḟ = F D(1 -iα) + D(1 -iβ) -1 + (δ + i)∇ 2 ⊥ , (1.30a) Ḋ = -bD s (I -I s ) , (1.30b) Ḋ = -rb Ds s (I -I s ) .
(1.30c)

These results

have not yet been tested numerically or experimentally but they represent an interesting possible development of the model in Eqs. (1.1), both to observe how the stability of the system and its possible solutions would be modified in the conservative limit and to verify if a Toda potential description can actually well approximate the evolution of the 2D+time system as the one in Eqs. (1.1).

Conclusions

In this Chapter we have studied the model for a monolithic broad-area VCSEL with an intracavity saturable absorber, analyzing into details all the possible solutions available for this kind of system. We then focused on the solution of extended 2D spatiotemporal chaos and investigated the presence of extreme events in this regime. This study required the development, for the first time, of an algorithm of detection for the spatiotemporal maxima in three dimensions (2D+time) in the transverse profile of the field intensity. Through statistical analysis we have shown the numerical existence of extreme events in the system, according to different possible definitions, illustrated the most suitable parameter region to observe them, also experimentally, and characterized their occurrence. A study on the spatial and temporal profile allowed us to assign to these events a typical size. The comparison of our work with others in the community [Coulibaly 2017] showed some differences which are fully justifiable when considering the numerical and experimental analysis have been, in those cases, developed in 1D + time systems. We believe, as suggested in Ref. [Bonazzola 2013], that 2D spatial effects play in fact a crucial role in the formation of extreme events. Finally we focused on a comparison in spatial and temporal size of these objects with the other solutions present in the system, namely cavity solitons, both stationary and self-pulsing. The striking similar shape between all these structures may imply an enhancement of extreme events when close to the dissipative soliton attractor. Future perspectives of this work consist in the investigation of predictability and in the study of the conservative limit for this model as well as its consequences on the observed extreme events.

Chapter 2

Semiconductor ring laser with injection Preface

This Chapter is dedicated to the study of semiconductor lasers with optical injection and one spatial dimension along the propagation, as in the case of optical fibers. The introductory part of the Chapter is devoted to the derivation of the model and its characterization. The first main section of the Chapter consists in the study of the results achieved in semiconductor lasers with injection in two different cavity configurations, where phase solitons have been experimentally observed [START_REF][END_REF], Gustave 2016b, Walczak 2017]. In particular we will focus on the formation of phase soliton complexes and on the possible collisions between phase solitons, that can lead to the development of extreme events. These results have been published in [Gustave 2017, Walczak 2017]. The second main section of the Chapter is dedicated to the study of abnormally high events in an unstable roll regime of the system. This regime has been extensively studied experimentally by F. Gustave [Gustave 2016a]. The results obtained numerically in this study and their comparison with the experiment have been published in [Rimoldi 2017a]. The main focus of the numerical and theoretical work carried out in previous papers regards the confirmation of the experimental observation of phase solitons in this system [START_REF][END_REF], Gustave 2016b]. The complete model is based on a set of effective Maxwell-Bloch equations to take into account the effective macroscopic polarization for a semiconductor material. This model considers the free propagation in the empty cavity space, and consists of three differential equations, for the electric field, the carrier density and the polarization of the semiconductor material. In particular in the complete model, instead of introducing any spectral filtering section to describe the material gain dispersion, a phenomenological formula for the gain and refractive index dependence on the frequency and carrier density has been introduced to fit the microscopic susceptibility: this leads to an evolution equation for the macroscopic polarization coupled with the equations for the electric field and the carrier density dynamics [Prati 2007b]. As we will see in the following the complete model is then reduced through the adiabatic elimination of the polarization. This reduction corresponds to assume a flat gain curve and even if it produces a model able to well describe phase solitons in the system, when the dynamics observed becomes more complex, the same process of "pixelation", observed in Chapter 1, takes place and a corrective diffusion term is to be phenomenologically inserted in the equation for the electric field to avoid self-collapsing.

More recently our group has started to focus its attention on different aspects of the interaction dynamics between phase solitons, including merging in phase soliton complexes [Gustave 2017] and collisions [Walczak 2017]. It is in this last regime that we found experimental and numerical evidence of extreme events rising from collisions. Furthermore we have been able, through the use of the same model, to observe experimentally and numerically abnormally high events in an unstable roll regime, where phase solitons are not a stable solution. Both kinds of highpeak events observed in this system seem to be strictly related to the chirality of the objects involved, as we will see in the following.

Model

Let us consider a ring resonator as the one schematized in Fig. 2.2, with a semiconductor active medium of length l = 1 mm placed in a unidirectional ring cavity of total length L = 1 m, driven by coherent optical injection. This schematization is slightly different from that illustrated as an example in Fig. 2.1 in order to grant an easier model derivation. We suppose that the two beam splitters (BS) have a transmittivity T = 0 and that the two mirrors (M) are perfect reflectors with T = 0. In the hypothesis of linear polarization, orthogonal to the direction of propagation, the electric field injected into the cavity at z = 0, Ẽ(0, t), and the electric field propagating in the active medium, Ẽ(z, t), can be written as follows

ẼI (0, t) = E I 2 exp (-iω 0 t) + c.c. Ẽ(z, t) = E (z, t) 2 exp [i (k 0 z -ω 0 t)] + c.c.
with E I ∈ R, k 0 is the wave vector such that k 0 = ω 0 n/c, where n = √ b is the background refractive index, ω 0 is the injected frequency and b is the dielectric constant. Similarly, for the macroscopic material polarization we can write

P(z, t) = P (z, t) 2 exp [i (k 0 z -ω 0 t)] + c.c. .
In the paraxial, slowly varying envelope and rotating wave approximations, the radiation-matter interaction can be described through the following set of partial differential equations [Prati 2007b]

∂E ∂z + 1 v ∂E ∂t = gP , (2.1a
)

∂P ∂t = 0 b A(N)E -B(N)P , (2.1b) ∂N ∂t = I eV - N τ e - i 4h (E * P -E P * ) , (2.1c)
where N is the carrier density in the semiconductor and τ e is the carrier density nonradiative decay time. I, V and e are respectively the pump current, the sample volume and the electric charge. The electric field propagation velocity in the medium is given by v = c/n b and g = iω 0 Γ c / 2 0 nn g c , with n g the group index and Γ c the field confinement factor. We suppose here that the Fourier components of the macroscopic polarization P (ω) and the electric field Ê (ω) are linked through the relation

P (ω) = 0 b χ 1 (ω) Ê (ω) , with χ 1 (ω) = A(N) B(N) -iω
where 0 is the vacuum permittivity (and b is the relative permittivity in the background medium). The functions A(N) and B(N) are complex and can be expressed in terms of the linewidth enhancement factor α, the transparency carrier density N 0 , the gain width Γ(N) and the detuning δ(N) between the gain peak and the reference frequency (for |δ(N)| Γ(N)). In particular we have [Prati 2007b]

A(N) B(N) = -f 0 (α + i) (N/N 0 -1) , with f 0 ∈ R B(N) = [Γ(N)(1 -iα) + 2iδ(N)] /τ d
where τ d is the dipole dephasing time. Hence Eqs. (2.1) can be recast as

∂E ∂z + 1 v ∂E ∂t = gP , τ d ∂P ∂t = [Γ(N)(1 -iα) + 2iδ(N)] [-iζ(1 -iα) (N/N 0 -1) E -P ] , τ e ∂N ∂t = Iτ e eV -N - iτ e 4h (E * P -E P * ) ,
where we highlight that ζ = f 0 0 b is the differential gain, with f 0 measuring the maximum gain [Prati 2007b].

Introducing the new variables

E = ηE , P = i ηP , D = ζ(N/N 0 -1)
where η2 = ζτ e /(2hN 0 ), we obtain the following set of equations

∂E ∂z + 1 v ∂E ∂t = -i gP , τ d ∂P ∂t = [Γ(D)(1 -iα) + 2iδ(D)] [(1 -iα)ED -P] , τ e ∂D ∂t = ζ I I 0 -1 -D - 1 2 (E * P + EP * ) ,
with I 0 = eV N 0 /τ e being the transparency current. The functions Γ(D) and δ(D) can be phenomenologically derived by a linear fit of the gain curve calculated with a microscopic model for different values of D and have the following expressions, as reported in [Prati 2007b],

Γ(D) = 0.276 + 1.016D δ(D) = -0.169 + 0.216D .
Setting then g = -i g, so that g ∈ R, and µ = ζ(I/I 0 -1) we obtain

∂E ∂z + 1 v ∂E ∂t = gP , (2.2a
)

τ d ∂P ∂t = [Γ(D)(1 -iα) + 2iδ(D)] [(1 -iα)ED -P] , (2.2b
)

τ e ∂D ∂t = µ -D - 1 2 (E * P + EP * ) . (2.2c)
Our objective is now to make the boundary conditions periodic and isochronous.

To do so [Lugiato 2015], referring to Fig. 2.2 we consider the field envelope at z = 0 = L E(0, t) = √ TE I + RE(l, t -∆t)e -iδ 0 , with R = 1 -T (in the hypothesis of no absorption by the beam splitters) and E I = ηE . Furthermore ∆t = (Ll)/c is the time required to the electric field to travel along the cavity outside the active medium and δ 0 = (ω cω 0 )Λ/c is the cavity frequency closest to ω 0 where Λ = Ll + nl is the effective cavity length.

If we introduce the transformation [Gustave 2016b]

η = z l , t = t + z l ∆t
the boundary condition assumes the isochronous form

E(0, t ) = √ TE I + RE(1, t )e -iδ 0 .
Furthermore since

∂ ∂z = 1 l ∂ ∂η + ∆t l ∂ ∂t ∂ ∂t = ∂ ∂t , Eq. (2.2a) becomes 1 l ∂E ∂η + ∆t l + n c ∂E ∂t = gP 1 l ∂E ∂η + Λ cl ∂E ∂t = gP ∂E ∂η + Λ c ∂E ∂t = glP , hence Eqs. (2.2) becomes ∂E ∂η + Λ c ∂E ∂t = glP , ∂P ∂t = 1 τ d [Γ(D)(1 -iα) + 2iδ(D)] [(1 -iα)ED -P] , ∂D ∂t = 1 τ e µ -D - 1 2 (E * P + EP * ) .
Introducing now the new field envelopes

E (η, t ) = E(η, t )e (ln R-iδ 0 )η + √ TE I η (2.3) P (η, t ) = P(η, t )e (ln R-iδ 0 )η (2.4) we obtain ∂E ∂η + Λ c ∂E ∂t = (ln R -iδ 0 ) E - √ TE I η + √ TE I + glP , ∂P ∂t = 1 τ d [Γ(D)(1 -iα) + 2iδ(D)] (1 -iα)D E - √ TE I η -P , ∂D ∂t = 1 τ e µ -D - 1 2 e -2 ln Rη E * P + E P * - √ TE I η P + P * ,
where we have used e (ln R-iδ 0 )η = (E -√ TE I η)/E, and the periodic isochronous boundary condition becomes

E (0, t ) = √ TE I + RE(1, t )e -iδ 0 = E (1, t ) .

Low transmission limit

We then apply the low transmission limit, defined as

T 1 , gl 1 , |δ 0 | 1
with the pump parameter A and the cavity detuning θ such that

A = gl T = O(1) , θ = δ 0 T = O(1) .
In this limit E and P as defined in Eqs. (2.3) and (2.4) can be approximated by E and P, respectively. Hence the set of equations gets reduced to c Λ ∂E ∂η

+ ∂E ∂t = 1 τ p -(1 + iθ) E + E I √ T + AP , ∂P ∂t = 1 τ d [Γ(D)(1 -iα) + 2iδ(D)] [(1 -iα)DE -P] , ∂D ∂t = 1 τ e µ -D - 1 2 (E * P + EP * ) ,
with the boundary condition E(0, t ) = E(1, t ), where we have introduced the photon decay time τ p = Λ/(cT) and used the fact that ln R/T = ln(1 -T)/T ≈ -1. Finally, setting the dimensionless time τ = t /τ d and the following substitutions

σ = τ d τ p , b = τ d τ e , y = E I √ T , we obtain cτ d Λ ∂E ∂η + ∂E ∂τ = σ [-(1 + iθ) E + y + AP] , ∂P ∂τ = [Γ(D)(1 -iα) + 2iδ(D)] [(1 -iα)DE -P] , ∂D ∂τ = b µ -D - 1 2 (E * P + EP * ) ,
and, renaming AP, AD, Aµ as P, D and µ respectively, the model becomes

cτ d Λ ∂E ∂η + ∂E ∂τ = σ [y -(1 + iθ) E + P] , (2.5a
)

∂P ∂τ = [Γ(D)(1 -iα) + 2iδ(D)] [(1 -iα)DE -P] , (2.5b) ∂D ∂τ = b µ -D - 1 2 (E * P + EP * ) , (2.5c)
with the boundary condition

E(0, τ) = E(1, τ) .
In particular, in summary, the variables E, P are defined respectively, unless renormalizations, as the slowly varying envelopes of the electric field and of the effective macroscopic polarization, and D is the excess of carrier density with respect to transparency. y is the amplitude of the injected field, θ is the frequency mismatch between the injected field and the closest empty cavity resonance and µ is the pump parameter.

As a final remark we also want to highlight that Maxwell-Bloch equations for multilongitudinal mode emission, similar to the ones here illustrated, were introduced for the first time for a two-level unidirectional ring laser with injected field in [Narducci 1985].

Rate-equation model

The model in Eqs. (2.5) can be reduced to a rate-equation model, if we assume a flat gain (that is a gain of infinite linewidth). This is performed by eliminating adiabatically of the polarization P, setting ∂P ∂τ = 0, so that the evolution of P results completely enslaved by the evolution of D and E P = (1iα)DE and reducing Eqs.

+ ∂E ∂τ = σ [y -(1 + iθ) E + (1 -iα)DE] , (2.6a
)

∂D ∂τ = b µ -D 1 + |E| 2 , (2.6b)
which can be written as

∂E ∂η + Λ cτ d ∂E ∂τ = T [y -(1 + iθ) E + (1 -iα)DE] , ∂D ∂τ = b µ -D 1 + |E| 2 .
Rescaling then time to the cavity roundtrip time Λ/c instead of the dephasing time τ d , we obtain ∂E ∂η

+ ∂E ∂τ = T [y -(1 + iθ) E + (1 -iα)DE] , (2.7a
)

∂D ∂τ = bT σ µ -D 1 + |E| 2 .
(2.7b)

This reduced model has been tested and compared with the complete one and has been proven capable to well describe the dynamics of the system in the regime of phase solitons and homogeneous solutions as we will also observe further into details in the next Section.

From a more numerical point of view, Eqs. (2.7) are less stiff than Eqs. (2.5) and their numerical integration is much faster.

Rate-equation model with a term of diffusion

A slightly different model that we will use for the larger part of this Chapter is the following

∂E ∂η + ∂E ∂τ -d ∂ 2 E ∂η 2 = T [y -(1 + iθ) E + (1 -iα)DE] , (2.8a
)

∂D ∂τ = bT σ µ -D 1 + |E| 2 . (2.8b)
where a diffusion term in the electric field has been introduced phenomenologically to take into account the finite gain linewidth. We have in fact observed that, simulating Eqs. (2.7) in more complex regimes than those of phase solitons or homogeneous states, the system tends to self-collapse, similarly to what we have already observed for the two-dimensional model of a laser with saturable absorber, studied in the previous Chapter. Even though the value of the coefficient d is arbitrarily assigned as the smallest value capable of avoiding self-collapsing, we will show in the next Subsection that it scales as the squared ratio of the cavity free spectral range to the gain linewidth, as proved in [Fedorov 2000] for a two-level medium.

Modified forced Ginzburg-Landau model

In order to determine the order of magnitude of the diffusion term in Eq. ( 2.8a), we will show here how a similar (complex) term can be derived in a rigorous way from the complete model with a nonstandard adiabatic elimination of the polarization valid when the laser is close to threshold.

We start from the complete model in Eqs. (2.5), reported here for convenience

cτ d Λ ∂E ∂η + ∂E ∂τ = σ [y -(1 + iθ) E + P] , ∂P ∂τ = ξ(D) [(1 -iα)DE -P] , ∂D ∂τ = b µ -D - 1 2 (E * P + EP * ) ,
where we have defined

ξ(D) = [Γ(D)(1 -iα) + 2iδ(D)]
. Dividing all the equations by

√ b, we obtain c √ τ d τ e Λ ∂E ∂η + 1 √ b ∂E ∂τ = σ √ b [y -(1 + iθ) E + P] , (2.9a) 1 √ b ∂P ∂τ = ξ(D) √ b [(1 -iα)DE -P] , (2.9b) 1 √ b ∂D ∂τ = √ b µ -D - 1 2 (E * P + EP * ) . (2.

9c)

Considering now = µ -1 ≈ 10 -2 as the smallness parameter, being µ thr = 1 the pump value at threshold, we can observe that in the weak injection limit also y = O( ). The physical value of polarization dephasing time is τ d ≈ 10 -13 s, and the photon and carrier lifetime result as τ p ≈ 10 -7 s and τ e ≈ 10 -9 s. Hence σ = τ d /τ p ≈ 10 -6 and b = τ d /τ e ≈ 10 -4 . We can then set

√ b = b 0 , σ √ b = σ 0 2 with b 0 = O(1) , σ 0 = O(1) .
Regarding the complex function ξ(D) we can observe that

ξ(D) √ b = [Γ(D)(1 -iα) + 2iδ(D)] √ b = ξ 0 (D) with ξ 0 (D) = Γ 0 (D)(1 -iα) + i ∆(D)
where

Γ 0 (D) = Γ(D) b 0 , ∆(D) = 2δ(D) √ b = O(1) .
Thus the model in Eqs. (2.9) can be approximated with

η 0 ∂E ∂η + 1 b 0 ∂E ∂τ = σ 0 2 [y -(1 + iθ) E + P] , 1 b 0 ∂P ∂τ = ξ 0 (D) [(1 -iα)DE -P] , 1 b 0 ∂D ∂τ = b 0 µ -D - 1 2 (E * P + EP * ) .
where we have defined η 0 = cτ d /(Λb 0 ). Rescaling the time variable as τ = √ bτ = b 0 τ we have

η 0 ∂E ∂η + ∂E ∂τ = σ 0 2 [y -(1 + iθ) E + P] , (2.10a) ∂P ∂τ = ξ 0 (D) [(1 -iα)DE -P] , (2.10b
)

∂D ∂τ = b 0 µ -D - 1 2 (E * P + EP * ) . (2.10c)
We assume that

∂E ∂τ = O(1) , ∂E ∂ (η/η 0 ) = η 0 ∂E ∂η = O(1)
which is reasonable considering that, for the time scales introduced and a cavity length Λ = 3 m, we have η 0 = 10 -3 . Slightly modifying the equation for P, we have

η 0 ∂E ∂η + ∂E ∂τ = σ 0 2 [y -(1 + iθ) E + P] , Γ 0 (D)(1 -iα) + i ∆(D) + ∂ ∂τ P = ξ 0 (D)(1 -iα)DE , ∂D ∂τ = b 0 µ -D - 1 2 (E * P + EP * ) .
Given the assumptions already stated, the equation for the electric field results We can then observe that at order 0 the effects of an injected field of order are negligible: D = µ/(1 + |E| 2 ) = 1 as in the free running laser, where |E| 2 = µ -1 = O( ) and the polarization follows adiabatically the dynamics of E. Then, in presence of injection such that y of order , we assume |E| 2 = µ -1 + O( 2) and it follows that D = 1 + O( 2), hence we can write D = µ -|E| 2 . Furthermore we can also set

∂E

ξ 0 (D) = ξ 0 (1) + O( 3 ) = Γ 0 (1)(1 -iα) + i ∆(1) + O( 3 ) ,
where we assumed that ξ 0 (D) varies slowly with D (with dξ 0 (D)/dD ≈ O( 2)).

Renaming for convenience ξ 0 (1) = ξ 1 , Γ 0 (1) = Γ 1 and ∆(1) = ∆ 1 and neglecting the terms of order 3 and higher, the equation for P becomes

Γ 1 (1 -iα) + i ∆ 1 + ∂ ∂τ P = ξ 1 (1 -iα)DE 1 + ∂/∂τ + i∆ 1 Γ 1 (1 -iα) P = 1 -iα + i ∆ 1 Γ 1 DE .
Then the solution for P truncated at order 2 is

P = 1 -iα + i ∆ 1 Γ 1 1 - ∂/∂τ + i∆ 1 Γ 1 (1 -iα) + 2 (∂/∂τ + i∆ 1 ) 2 Γ 2 1 (1 -iα) 2 DE = 1 -iα - ∂/∂τ Γ 1 + 2 ∂ 2 /∂τ 2 + i∆ 1 ∂/∂τ Γ 2 1 (1 -iα) DE = (1 -iα)DE -Γ 1 ∂E ∂τ + 2 (1 + iα) Γ 2 1 (1 + α) 2 ∂ 2 ∂τ 2 + i∆ 1 ∂ ∂τ E = (1 -iα) µ -|E| 2 E -Γ 1 ∂E ∂τ + 2 (1 + iα) Γ 2 1 (1 + α) 2 ∂ 2 ∂τ 2 + i∆ 1 ∂ ∂τ E ,
where we have observe that D is constant and D = 1 ≈ µ -|E| 2 . Finally, noticing that, from the equation for the electric field, at order

∂E ∂τ = -η 0 ∂E ∂η ,
we can also write

P = (1 -iα) µ -|E| 2 E -Γ 1 ∂E ∂τ + 2 (1 + iα) Γ 2 1 (1 + α) 2 η 2 0 ∂ 2 ∂η 2 -iη 0 ∆ 1 ∂ ∂η E .
This is the term given by the nonstandard adiabatic elimination of P and it is to be inserted in Eq. (2.10a). The remaining first temporal derivative is negligible since its contribution, once inserted in Eq. (2.10a), would be of order 3 and a temporal derivative is already present in the equation. We then obtain

η 0 ∂E ∂η + ∂E ∂τ = σ 0 2 y -(1 + iθ)E + (1 -iα) µ -|E| 2 E + 2 (1 + iα) Γ 2 1 (1 + α) 2 η 2 0 ∂ 2 ∂η 2 -iη 0 ∆ 1 ∂ ∂η E = σ 0 2 y -(1 + iθ)E + (1 -iα) µ -|E| 2 E + b(1 + iα) Γ(1) 2 (1 + α) 2 η 2 0 ∂ 2 ∂η 2 -iη 0 ∆ 1 ∂ ∂η E ,
where we can observe that the terms on the first line are given by a standard adiabatic elimination of P and D with the approximation that D ≈ 1/(1 + |E| 2 ) ≈ 1 -|E| 2 , while the term on the second line comes from the nonadiabatic elimination of P. In the original variables of space η and dimensionless time τ we have then

cτ d Λ ∂E ∂η + ∂E ∂τ = σ y -(1 + iθ)E + (1 -iα) µ -|E| 2 E + d c 2 τ 2 d Λ 2 ∂ 2 ∂η 2 -2i cτ d δ(1) Λ ∂ ∂η E ∂E ∂η + Λ cτ d ∂E ∂τ = T y -(1 + iθ)E + (1 -iα) µ -|E| 2 E (2.11) + d σ 2 T ∂ 2 ∂η 2 -2iσδ(1) ∂ ∂η E
where we have substituted η 0 with its expression in terms of the different timescales

η 0 = cτ d /(Λ √ b) and ∆ 1 = 2δ(1)/ √ b, set d = 1 + iα Γ(1) 2 (1 + α 2 )
and remembered that σ = τ d /τ p and τ p = Λ/(cT). Eq. ( 2.11) is a modified forced Ginzburg-Landau equation, that has been introduced in [START_REF][END_REF]] to describe the dynamics of a phase soliton in this kind of system. Unfortunately, the description of this system through just one equation for the electric field, as in Eq. ( 2.11), has been shown to lack the system asymmetry typically observed experimentally, where only structures carrying positive chiral charges result stable: in fact, even if Eq. ( 2.11) lacks the parity symmetry η → -η of the standard Ginzburg-Landau, due to the term in δ(1), unless this parameter assumes an unrealistically large value, the equation still supports both chiral charges [Gustave 2015].

Ordering the equation so that all the derivatives are on the left hand side we have

1 + 2i dσδ(1) ∂E ∂η + Λ cτ d ∂E ∂τ - dσ 2 T ∂ 2 ∂η 2 E = T y -(1 + iθ)E + (1 -iα) µ -|E| 2 E ,
and, rescaling the time on the dimensionless ratio of cavity roundtrip time to the polarization dephasing time τ d , we obtain the equation in the same variables of Eqs. (2.8)

1 + 2i dσδ(1) ∂E ∂η + ∂E ∂τ - dσ 2 T ∂ 2 ∂η 2 E (2.12) = T y -[1 -µ + i (µα + θ)] E -(1 -iα)|E| 2 E ,
Here we can observe that the diffusion coefficient is complex but its order of magnitude corresponds more or less to the square of the ratio of dephasing time to roundtrip time, unless an extra factor of order 10 -1 given by T/ √ 1 + α 2 ≈ T/α.

Rate-equation model with a Ginzburg-Landau term

To preserve the asymmetry of the experimental system we observed that it is necessary to keep the equation for the dynamics of the carrier density D, as the main mechanism for the chiral selection and the strong η → -η symmetry breaking lies in the non-instantaneous medium dynamics. To do so we develop the same kind of approximation as in the previous Subsection but we avoid to adiabatically eliminate D, thus obtaining 1 + 2i dσδ 0 (1) ∂E ∂η

+ ∂E ∂τ - dσ 2 T ∂ 2 ∂η 2 E = T [y -(1 + iθ)E + (1 -iα)DE] ∂D ∂τ = Tb σ µ -D(1 + |E| 2 ) . (2.13)
The model in Eq. ( 2.13) is equivalent to the complete model of Eqs. (2.5) since it maintains the chirality asymmetry and sets a filter on the spatial modes for 2.17) for the same parameters. more complex dynamics. The only problem consists in the fact that such filter, being that of the complete model, is not very efficient and would force us to run longer simulations on larger grids (or consider larger values of σ) in order to well describe the system under study, especially in complex regimes. It is for this reason that we preferred to develop the studies in the next Sections through the model in Eqs. (2.8), which provides both asymmetry and a more efficient filter.

HSS and linear stability analysis

The complete model in Eqs. (2.5) admits the following homogeneous stationary solution (HSS)

E s = ρ s e iφ s , P s = (1 -iα)DE , D s = µ 1 + ρ 2 s (2.14) with y 2 = ρ 2 s (1 -D s ) 2 + (θ + αD s ) 2 (2.15) φ s = arctan θ + αD s D s -1 (2.16)
where µ = 1 corresponds to the threshold value for the pump parameter in the free running laser. In Fig. 2.3 we plotted the stationary intensity I = ρ 2 s of the electric field as a function of the injection amplitude y as indicated in Eq. ( 2.15) for values of θ and µ similar to those studied in the rest of this Chapter. The value of linewidth enhancement factor α = 3 has been kept constant in all the results here shown. The equation for the bistability boundaries is [Gustave 2016b]

θ ± = - µα ± µ 2 ρ 4 s (1 + α 2 ) -(1 + ρ 2 s ) 2 -µ 2 (1 + ρ 2 s ) 2 , ( 2
.17) from which we can derive that the rightmost point of the bistability domain is given by θ = -µα/(1 + ρ 2 s ) 2 . Hence in order to have bistability the detuning θ needs to be negative. In Fig. 2.3 we have plotted Eq. ( 2.17) for the same parameters of the stationary intensity in (a). As in the previous Chapter, we now study the stability of the stationary homogeneous solution considering small spatiotemporal perturbations from the stationary values of E s , P s and D s , modulated by a longitudinal wavevector k z and exponentially growing (decaying) in time, so that:

      E E * P P * D       =       E s E * s P s P * s D s       + e λt+ik z z       δE δE * δP δP * D       , ( 2
      δE δE * δP δP * δD       = J (2.5)       δE δE * δP δP * δD       (2.19)
where J (2.5) is the Jacobian matrix

      -i cτ d Λ k z -σ(1 + iθ) 0 σ 0 0 0 -i cτ d Λ k z -σ(1 -iθ) 0 σ 0 ξ(D s )(1 -iα)D s 0 -ξ(D s ) 0 ξ(D s )(1 -iα)E s 0 ξ * (D s )(1 + iα)D s 0 -ξ * (D s ) ξ * (D s )(1 + iα)E * s -b 2 P * s -b 2 P s -b 2 E * s -b 2 E s -b       , with ξ(D s ) = Γ(D s )(1 -iα) + 2iδ(D s ).
The rate-equation model in Eq. ( 2.7) and the rate-equation model with diffusion in Eq. (2.8) present the same homogeneous stationary solutions for the electric field and carrier density of the complete model in Eqs. (2.14), (2.15) and (2.16). Considering small spatiotemporal perturbations from the stationary values of E s , E * s and D s , modulated by the longitudinal wavevector k z , we have

  E E * D   =   E s E * s D s   + e λt+ik z z   δE δE * D   , ( 2
-ik z + T [D s -1 -i (αD s + θ)] 0 T(1 -iα)E s 0 -ik z + T [D s -1 + i (αD s + θ)] T(1 + iα)E * s -Tb σ D s E * s -Tb σ D s E s -Tb σ (1 + |E s | 2 )   ,
and the Jacobian matrix J (2.8) for the rate equation model with diffusion defined as follows

  -ik z -dk 2 z + T [D s -1 -i (αD s + θ)] 0 T(1 -iα)E s 0 -ik z -dk 2 z + T [D s -1 + i (αD s + θ)] T(1 + iα)E * s -Tb σ D s E * s -Tb σ D s E s -Tb σ (1 + |E s | 2 )   ,
An analytic stability analysis on the complete model, using mode expansion, is reported in [Gustave 2016b] where, under the reasonable assumptions of b = O( ), σ = O( 2) and the modes α n = 2πncτ d /Λ = O( ), the system is solved perturbatively in . This assumption implies that σ b and, together with the fact that µ thr = 1, brings, in the single mode limit, the model close to its class-A limit. Noticing that the eigenvalues λ can be written as λ = λ 1 + λ 2 2 with λ 1 imaginary, the instability condition results simplified to Re(λ 2 ) > 0, which gives a biquadratic equation in α n for the boundary of the stability domain

c 4 b -4 α 4 n + c 2 b -2 α 2 n + c 0 = 0
with coefficients in a convenient and short form as reported in [Gustave 2016b].

Setting c 0 = 0, we observe that in the single mode limit n = 0

c 0 = ρ 2 s + 1 2 -µ 2 dy 2 d (ρ 2 s )
and, for a typical S-shape stationary curve as that in Fig. 2.15, c 0 = 0 where the derivative is zero, that is in the turning points which means, that the negative slope part of the stationary curve is, as usual, unstable. Still in the single mode limit, the other solution ρ 2 s = √ µ -1, which exists even when the stationary curve is single valued, corresponds to the injection locking threshold in the class-A laser limit: in fact the stationary curve is unstable when the injected amplitude is so small that ρ 2 s < √ µ -1.

Comparison of the different models

In the present Subsection we will focus on a numerical comparison between the three different models given by Eqs. (2.5), (2.7) and (2.8). We will not consider the model of rate equations with a Ginzburg-Landau term as in Eq. ( 2.13) since the instability domain regions do not basically differ from the complete model, at least for the range of modes we are interested in.

The numerical method we use to compare the different models consists in computing, for each mode k z,n , the eigenvalues λ of the Jacobian matrix and select the largest real part between these eigenvalues max(Re(λ)). For each n-th mode, the value of the stationary intensity for which max(Re(λ)) changes its sign, gives the instability boundary. In all the results here shown we have kept constant the parameters α = 3, σ = 3 × 10 -6 , and T = 0. which we have observed to be the minimal value of d necessary to avoid selfcollapsing. As we can observe the model of pure rate equations does not introduce a limit to the number of unstable sidemodes, such limit is instead present in both the complete model and the model with diffusion, even if with different ranges.

The stationary curve relative to this particular choice of parameters is depicted in Fig. 2.5. For a chosen value of injection y in the triple-valued region of the stationary curve, as highlighted in Fig. 2.5, there are three fixed points A, B and C: in particular A is a node, B is a saddle and C is a focus. From Eqs. (2.14),(2.15) and ( 2.16) it is clear that for the free running laser (y = 0) the stationary equation is reduced to D s = 1 (which implies |E s | 2 = µ -1) and θ = -α, which defines the frequency of the free running laser. For definition

θ = ω c -ω inj k , -α = ω c -ω L k ,
with ω c cavity frequency, ω inj injection frequency and ω L , laser frequency. It then results that

ω inj -ω L = -k(θ + α) .
Hence θ + α represents the frequency mismatch between the driving field and the free running laser [Spinelli 1998]. As we will see further into detail in the next Section, the necessary condition for the existence of a phase soliton, consists in the stability of the whole upper branch: in particular for θ + α large and positive the upper branch results stable, instead for θ + α small or even negative, a part of the upper branch results unstable against a band of sidemodes as occurs in Fig. 2.4 [Gustave 2016b]. Furthermore, when θ + α is negative, the instability of the upper branch starts from the turning point, instead, for θ + α positive (but small), it starts from a different point in the upper branch. In Fig. 2.6 we depicted the stationary curves (a) and (b) and the instability regions (c) and (d) of the rate-equation model with (in red) and without diffusion (in blue) for different values of parameters, that will be also considered later in the Chapter. In particular, Figs. 2.6(a) and (c) were depicted for µ = 1.1, θ = -2.7 and b = 0.0005, and Figs. 2.6(b) and (d) were depicted for µ = 2, θ = -2.7 and b = 0.0001. In both these cases, the necessary condition for the existence of phase solitons in the system is satisfied, since the upper branch results stable. In Fig. 2.6(c) and (d) we did not report also the instability domain of the complete model since, for the range of modes here depicted it overlaps with the rateequation instability domain. In Fig 2 .7 we depicted a comparison between all three different models for a smaller value of the diffusion coefficient d = 10 -12 . The other parameters are the same as in Fig. 2.4. We can clearly observe that for lower values of d the model of Eq. (2.8) approximates better the instability domain of the complete model: this is because, similarly to what observed in Subsection 2.2.5, the order of magnitude for the diffusion coefficient in the complete model corresponds to

T √ 1 + α 2 τ d Λ/c 2 ≈ 1 10 10 -13 10 -8 2 = 10 -11
In fact, more rigorously, for a cavity length of 1 m, when assuming a gain linewidth of 10 THz, the number of modes under the gain line should be much larger than what depicted in Fig. 2.4(c), leading to a value of d smaller by some order of magnitude. In this work we chose to keep d = 10 -6 in order to filter the higher order modes that are not relevant for the dynamics: for this reason we could then limit the number of grid points in the numerical simulations to a reasonable value (1024 points) and at the same time avoid the self-collapsing effect that we would observe for lower choices of d. Finally we would like to stress out that our choice for d is also justified when considering that all the simulations here presented have been run for values of injection close to the turning point in the stationary curve: as we can observe from the instability plots in Fig. 2.4 and 2.6, in the proximity of the turning point, a larger value of d (for which the instability domain would appear more similar to that of the rate equation model) would affect only slightly the stability of the upper branch.

Finally in Fig. 2.8 and 2.9 we reported, respectively the stationary curves and the instability domains of Eqs. (2.7) and (2.8) for four different values of µ and θ. On the two rows θ is decreasing from left to right (i.e. moving farther from the resonance condition θ = -3) while on the two columns µ is increasing from top to bottom. We can observe that the farther we move from the threshold value µ thr = 1 the more modes are falling under the gain line for the model with diffusion.

Experimental setup

The experimental setup is a macroscopic ring laser of length Λ ≈ 1 m with a 4-mm-long, 980-nm semiconductor optical amplifier (SOA), acting as gain medium and set inside a free space optical cavity [START_REF][END_REF]]. This setup has been studied into details by F. Gustave in his Ph.D. thesis [Gustave 2016a]. A schematization of this experimental setup is depicted in Fig. 2.10. The ring cavity is eight-shaped to allow a lower incidence angle of the beam on the cavity mirrors. An anti-reflection coating on the semiconductor optical amplifier pre- vents from the self-lasing of the semiconductor element and an optical isolator (OI), placed to the opposite side of cavity with respect to the gain medium ensures the unidirectionality of the beam. The ring cavity itself consists of three high reflectivity (R > 99%) mirrors M and one beam splitter (BS) of reflectivity R ≈ 90%, that allows the injection of the external coherent field. The beam splitter used to provide the ring laser output has instead a reflectivity R ≈ 10%. The coherent field comes from a tunable single-mode laser where it is possible to broadly control the injection frequency and thus the detuning between the injection and the ring laser. Furthermore the presence of a piezoelectric actuator (piezo) set on one of the mirror mounts inside the cavity allows to adjust the cavity length, hence giving the possibility of a fine control of the cavity length and thus of the detuning [Gustave 2016b]. The laser output is measured by means of fast photodiodes connected to a 12.5 GHz bandwidth real-time oscilloscope with 100 GS/s sampling rate. A heterodyne setup based on a 3×3 fiber coupler allows to access both the power and the relative phase between the forcing and the ring laser [START_REF][END_REF], Gustave 2016b]. The second optical setup, studied, alternatively to the main one, by the experimental team at INPHYNI is depicted in Fig. 2.11. In this configuration the slave laser is (a semiconductor laser) in a Fabry-Perot configuration. The active medium is antireflection-coated only on one side, with the other side being highly reflective. Once again the introduction of an optical isolator (OI) prevents any reflection from the slave to the master laser (which provides the input beam). M are the mirrors that enclose the cavity and present a reflectivity R > 99% and BS are two beam splitters with reflectivity R ≈ 90%, inserted to provide an input for the forcing field and an output for the detection of the emitted field. The output can be observed through the same technique of the previous setup, allowing access to the information on the phase of the field or, when this is not necessary, through a high sensitive photodiode connected to a fast oscilloscope, to obtain information just on the laser power [Walczak 2017]. We will compare our numerical results with the experimental ones obtained in these two configurations. Even if these two setups differ by their geometry in a relevant way, they have been found to present the same qualitative features [Gustave 2017]. In particular both of them are very strongly multimode semiconductor lasers with an external forcing. Furthermore even if in the Fabry-Perot configuration the presence of forward and backward fields is evident, no qualitative differences have been found with the case of ring geometry, at least in the regimes studied.

Phase solitons and complexes

All the simulations here shown have been run with a split-step numerical method, similar to the case of the previous Chapter. Further details on the numerical method can be found in Appendix F. As anticipated in the linear stability analysis Section, a necessary condition for the existence of phase solitons (PSs) is the stability of the upper branch of the stationary curve, as the one depicted in Fig. 2.5: this requirement is satisfied if θ + α is large and positive, provided that the amplitude of the injected field is such that the fixed points A and B are sufficiently close (that is, the system is close enough to the left turning point of the stationary curve) [Gustave 2016b]. Numerically, a phase soliton can be obtained setting as initial condition the locked state corresponding, for fixed injection, to the point A on the upper branch, and superimposing a positive phase kink of 2π along z for the field E: in particular the form we used for the phase kink is

Φ + (η) = 4 tan -1 [exp (-βη)],
with β large to have a steep kink [Gustave 2016b]. When the phase soliton is stable, the final state corresponds to a traveling pulse with a 2π chiral charge.

Here and in most of the work presented in this Chapter we depicted our results in a spatiotemporal representation. In Fig. 2.12 we plotted the spatiotemporal diagram of 2000 cavity roundtrips in a generic unstable roll pattern in two different reference frames. The first reference frame, in Fig. 2.12(a), coincides with the experimental one, where the data is acquired in a specific point of the cavity (here at η = 1, that is z = L) as it would be when placing a detector right at the output of the cavity. The second reference frame, in Fig. 2.12(b) is obtainable only numerically and corresponds to acquiring, at each roundtrip τ = nτ c with n ∈ N, the entire spatial intensity profile along the cavity, and then stacking all the recorded profiles (much like snapshots taken every roundtrip time τ c ) onto each other, so to obtain the spatiotemporal diagram illustrated. We can clearly notice that the two spatiotemporal diagrams are mirror images one of the other, as it is expected in the case of any propagating structure. Even if the true spatiotemporal diagram is given in Fig. 2.12(b), given this similarity, from now on we will just represent the numerical results in the reference frame (a), to facilitate the comparison with the experimental data.

In [START_REF][END_REF], Gustave 2016b] the numerical analysis on phase solitons was made under the assumption that the laser was just a few percent above threshold, with a pump parameter µ ranging from 1.01 to 1.03. Here we studied phase solitons in a configuration more similar to the experimental case, where the laser was about 10% above threshold [Gustave 2017]. Furthermore farther from the laser threshold the dynamics tends to be faster and eventual merging or collision of structures occur on a shorter temporal scale, which made it easier to observe these phenomena. In this Section we will use the pure rate-equation model as in Eq. (2.7), since it was shown [Gustave 2016b] that phase solitons do not require a refined description of the active medium, which would be provided with the complete model, and they are correctly described when the macroscopic polarization is adiabatically eliminated.

Parameter choice

In order to search for the correct parameters to observe phase solitons in this new configuration we are going to make some theoretical considerations [Gustave 2017]. Starting from Eqs. (2.6) we assume that the excess pump with respect to the threshold µ -1 = remains small. We then introduce the field amplitude F and the injection amplitude y 0 defined so that

E = (µ -1) 1/2 F = 1/2 F , y = (µ -1) 3/2 y 0 = 3/2 y 0 .
Then Eq. (2.6a) for F becomes

cτ d Λ ∂F ∂η + ∂F ∂τ = σ [(µ -1)y 0 -(1 + iθ) F + (1 -iα)DF] cτ d Λ ∂F ∂η + ∂F ∂τ = σ [(µ -1)y 0 -(1 + iθ) F + (1 -iα)(D -µ + µ)F] cτ d Λ ∂F ∂η + ∂F ∂τ = σ {(µ -1)y 0 + [µ -1 -i(θ + µα)] F + (1 -iα)(D -µ)F} cτ d Λ ∂F ∂η + ∂F ∂τ = σ y 0 + (1 -i∆) F + (1 -iα) D -µ F , (2.21) with ∆ = θ + µα µ -1 = θ + µα .
It is reasonable to assume that, in the limit of a laser very close to the threshold, D deviates little from the unsaturated pump value µ by some correction of order :

D = µ + D 1 with D 1 = O(1) .
We can then rearrange Eq. ( 2.21) as where we can observe that all the terms in the square brackets are of order unity. From Eq. ( 2.6b) the equation for D 1 becomes

cτ d Λ ∂F ∂η + ∂F ∂τ = σ [y 0 + (1 -i∆) F + (1 -iα)D 1 F] , ( 2
∂D 1 ∂τ = b µ -(µ + D 1 ) 1 + |F| 2 (2.23)
and, at the leading order, under the assumptions made, the original set of rate equations in Eqs. (2.6) can be written as

cτ d Λ ∂F ∂η + ∂F ∂τ = σ [y 0 + (1 -i∆) F + (1 -iα)D 1 F] ∂D 1 ∂τ = -b D 1 + |F| 2 ,
where we have used that µ ≈ 1 in the second equation. We can here notice that, if the effective detuning ∆ and the injected field y 0 are kept constant then the only role of the smallness parameter is to determine the decay rate of the electric field: this observation gives us a hint on how to change the parameters to find PSs for any value of µ, provided we know their stability range for a given pump value. For instance in [START_REF][END_REF]] it was reported that a stable PS exists for the parameters µ = 1.01, α = 3, θ = -2.97 and for a small range of injection values around y = 1.4 × 10 -3 , which corresponds to a value of ∆ = 6 and y 0 = 1.4. For a different value of µ we expect to find stable PSs for θ = -α + (6α) in an interval around y = 1.4 3/2 . Hence, for instance, for µ = 1.1 and α = 3 we expect stable PSs for θ = -2.7 and a range of injection values around y = 4.43 × 10 -2 [Gustave 2017].

All the simulations in this Section have been run setting µ = 1.1, θ = -2.7, b = 0.0005 and, as usual, α = 3, σ = 3 × 10 -6 , T = 0.3.

Single-charge phase soliton

In Fig. 2.13 we plotted the homogeneous stationary solution for the laser intensity as in Eq. ( 2.15) and superimposed the phase soliton stability domain, in particular indicating the maxima (square markers) and minima (triangle markers) of each PS profile for different values of injection. The instability domain for this specific set of parameters has been reported in Fig. 2.6(c), where we can observe that, as required, the whole upper branch of the stationary curve results stable. From Fig. 2.13 we observe that the interval of stability for phase solitons is 4.45 • 10 -2 ≤ y ≤ 4.69 • 10 -2 : the small difference with the theoretical prediction may be associated to the fact that, first of all, the prediction is accurate up to a order , which implies some small correction for = 0.1 as in our case. Furthermore we have observed that is involved in the definition of the decay rate of the electric field: since the stability of PS depends also on the ratio of the decay rates of the electric field and the gain, an increase in , even if small, can any way partially affect the stability of the phase solitons. It is interesting to observe, in the zoom depicted in Fig. 2.13(b), that, at the boundaries of the stability domain, the minimum intensity of the phase soliton results very close to that of the lower homogeneous state. In Fig. 2.14 we depicted the spatiotemporal diagrams (a,b), the trajectories in the complex plane (Re(E), Im(E)) (c,d) and the intensity temporal profiles (e,f) at the boundaries of the PS stability domain in y: in particular the left column of Fig. 2.14 corresponds to y = 0.0445 and the right column corresponds to y = 0.0469. In Fig. 2.14(c,d) we also plotted the three fixed points, defined for a given value of y in the triple-valued region of the stationary curve, as already illustrated in Fig. 2.5: in particular for this choice of parameters A is a stable node, B is a saddle and C is an unstable focus. We highlight that, due to the conditions for the existence of PSs, these three fixed points are always defined in the presence of stable phase solitons: during a single roundtrip, the trajectory of the system, in the Argand plane (Re(E), Im(E)), starts from the stable homogeneous stationary solution, given by the point A, then gets repelled by the point B and circles around, or passes near, C (in this projection).

Due to the different pump value with respect to previous studies [Gustave 2016b], we were here able to observe that the trajectory described by the system in the complex plane does not necessarily need to include the origin since the point that has an actual role in the dynamics is C and not (Re(E), Im(E)) = (0, 0). For a lower value of pump this difference was less evident since the value of intensity for C was even closer to zero. Furthermore we can observe that, if the dynamics was restricted to the Argand plane, the trajectory would always circle around C, nevertheless from Fig. 2.14(d) we can notice that there are exceptions to this rule, because the actual phase space includes a third dimension represented by the variable D as illustrated in Fig. 2.15. It is in fact for this reason that, for the right extremum of the stability domain, the minimum intensity of the PS results very close to that of C on the lower homogeneous branch. From Fig. 2.15 we can also observe that the main role of the point C is of repulsive nature, since the trajectory of the PS remains actually bounded to lower values of carrier density. The phase portraits in Fig. 2.14(c,d) seem in good agreement with the results by Chaté et al. in [Chaté 1999] for phase-charged solitons observed with the forced complex Ginzburg-Landau equation in the condition of small forcing. The results here shown for single-charge phase solitons are in very good agreement with the experiment [START_REF][END_REF], Gustave 2016b], where PSs carrying a topological charge with these specific characteristics have been observed, as depicted in Fig. 2.16, where, in particular, on the right panel the (logarithmic) color scale indicates the frequency of the observations. In the spatiotemporal diagrams depicted in Fig. 2.14(a,b) we can observe that a structure propagating on a vertical line would correspond to a pulse that moves at velocity c, that is, a structure that covers the cavity length Λ in a roundtrip time τ c = Λ/c. Hence the more a line describing the propagation of a structure is inclined towards the right side, the slower the object is propagating. From Fig. 2.14(a,b) and (e,f) we can then notice that for higher values of injection y the PS increases its velocity and becomes shorter in time. In particular, a delay of ∆τ accumulated over N roundtrips implies, for the PS, a velocity

v = NΛ Nτ c + ∆τ = Λ/τ c 1 + ∆τ/(Nτ c ) = c 1 + ∆τ/(Nτ c ) , so that v c ≈ 1 - 1 N ∆τ τ c .
We can define the width of a PS as the temporal interval between points 1 and 2 in Fig. 2.14(e,f). Points 1 and 2 are defined, respectively, as the points on the PS temporal profile with intensity I min + I bg /2 and I max + I bg /2, with I min , I max and I bg the minimum, maximum and background intensity values. In Fig. 2.17(a) we can observe the behavior of the PS velocity and width as a function of the injected amplitude y: as we can observe the two quantities present an almost linear dependence on the injection. As depicted in Fig. 2.17(b), this implies also a linear decrease of the velocity with the width of the PSs. Walczak in our team, using the Fabry-Perot experimental setup illustrated in Fig. 2.11.

Attractive interaction

As depicted in Fig. 2.18, in the experiment, there are configurations where chaotic regions coexist with multiple single-charge PSs. In this case, especially when two PSs are separated by just a few times their width, a weak attractive interaction may occur, causing the PSs to get closer to each other, experiencing a stronger interaction the more they approach, and leading to their merging into a double-charge PS. One particular feature illustrated in Fig. 2.18 consists in the fact that after the merging the double-charge PS propagates at a smaller velocity than the previous single-charge phase solitons, as highlighted by the black dashed line: a fact that will be better investigated in the next Subsection.

The numerical simulations exhibit the same kind of merging, thus confirming the experimental description of the phenomenon. In Fig. 2.19 we reported nine cases of consecutive merging of single-charge phase solitons with other multiple-charge phase solitons for a fixed value of injection y = 0.0458 in the middle of the PS stability domain: in all cases a single-charge soliton is gener- ated at a similar distance (corresponding approximately to 20% of the roundtrip time) from a soliton of charge n -1. Due to the different speed of the PSs, the two objects merge, forming a PS of charge n. The objects resulting from the interaction appear stable after the merging and propagate in the cavity without any change over a long time. We can also observe, as in the experimental case, that, of the two structures interacting, it is usually the one that is generated at a following time that undergoes the most attraction to the previous PS, thus bending towards it. The times indicated on the spatiotemporal diagram of Fig. 2.19 do not give any indication regarding the interaction time, since, depending on the figures, the initial distance between the two objects may not be necessarily exactly the same. In Fig. 2.20 we show two special cases of merging, for instance in Fig. 2.20(a) we show that, for this specific choice of parameters, the merging of a PS of charge 10 with a single-charge PS does not lead to a PS of charge 11 but to a PS of charge 17, as six additional charges are generated as a result of the interaction. Furthermore (b) trying to merge the PS of charge 17 with another single-charge PS, the system finally looses its stability and a more complex dynamics arises from the collision, as the system falls locally on a turbulent state, which corresponds to the (unstable) lower branch of the stationary curve, coexisting for this parameter choice with the stable upper branch. Eventually this regime spreads to the all cavity. The resulting dynamics is then not anymore well described by the pure rate-equation model.

(a) PS 1 + PS 1 → PS 2 (b) PS 2 + PS 1 → PS 3 (c) PS 3 + PS 1 → PS 4 (d) PS 4 + PS 1 → PS 5 (e) PS 5 + PS 1 → PS 6 (f) PS 6 + PS 1 → PS 7 (g) PS 7 + PS 1 → PS 8 (h) PS 8 + PS 1 → PS 9 (i) PS 9 + PS 1 → PS 10
In Fig. 2.21 we depicted three cases of merging (for y = 0.0458) between two single-charge PSs starting at an initial distance (dist ini ) of, respectively, 5% (a), 7% (b) and 10% (c) of a roundtrip time τ c from each other: as we can observe the interaction between the two objects occurs at longer times for larger distances. Furthermore, as in Fig. 2.19, we can distinguish a long initial stage where the interaction results very weak, the distance between the two PS is barely affected and the change in velocity of any of the two objects is almost not detectable: it is only when the two PSs are sufficiently close that the velocity of the second pulse is abruptly affected giving rise to the merging.

In Fig. 2.22, we reported the dependence of the merging times on the initial distance at which the two phase solitons are generated, varying in particular the initial distance from 5% to 10% of the roundtrip time τ c . The merging time seems to increase exponentially, as shown in semi-logarithmic scale in Fig. 2.22(b).

In this study no critical distance has been found that causes the interaction between the two PSs to stop, but the numerical validation of this statement would require extremely long simulations. For instance, from Fig. 2.22(a) we can notice that more than 10 6 roundtrips would be necessary to observe the interaction between two PSs set at an initial distance of 11%. Weak and strong coupling and their role in the formation of 2D spatial soliton complexes have been analyzed in [Rosanov 2005]: the extension of this concept to our 1D longitudinal soliton complexes could be an interesting development of the work here presented. However there would for sure be a few conceptual differences: for instance in [Rosanov 2005] the system was phase-invariant, allowing the phase to freely vary during the evolution of the system. In our case instead the presence of a coherent injection fixes the reference phase, allowing only multiples of 2π phase jumps in the system. Furthermore the non-instantaneous active medium response permits just one sign for the chiral charges of phase solitons, thus limiting the degrees of freedom of our system (see also the end of Section 2.7.4).

Multiple-charge phase solitons

As already observed in Fig. 2.19, multiple-charge phase solitons are stable in this system. In Fig. 2.23 we illustrated some examples of the possible multiplecharge PSs that can be observed numerically: in particular (a) is a PS of charge 6, (b) is a PS of charge 10 and (c) is a PS of charge 17. In particular PSs stable up to charge 17 have been obtained and n-charge PSs with n = 2, ..., 10 can be obtained from collisions of (n -1)-charge PSs with a single-charge PS. These observations are in very good agreement with the experimental results illustrated in Fig. 2.24 where PSs of charge 6 (a) and 11 (b) are depicted. We would like to point out that the two experimental measurement were not performed for exactly the same parameters, as it can be also noticed by the different background intensity values. This can be due to the possible small alterations of the cavity length, thus the different slopes of the experimental PSs are not to be considered a good indicator for the velocity for these objects. The fact that in the numerical simulations it was not possible to obtain a PS of charge 11 is due to the possibly different value of injection used, which may affect the stability of multiple-charge PSs more than that of a single-charge phase soliton. From Fig. 2.23 we can observe that there is a clear dependence of the PS velocity on the number of charges carried. This dependence is illustrated in Fig. 2.25, where we can observe that PSs with larger charge lag behind in the copropagating frame, reaching a plateau speed at approximately charge 4: which is intuitive, in terms of radiation-carrier interaction, since we are to expect larger structures to be dragged more by the active medium inertia.

Given the uncertainty on the exact setup cavity length it results difficult to compare the velocity of the experimental PSs with different charges when they are not measured simultaneously: for this reason it was preferred to observe the evolution of the same multiple-charge PS during the course of time, when affected by thermal fluctuations of the active medium (see [Gustave 2017]). These fluctuations provoke an enlargement or a compression of the multiple-charge phase soliton size, which, depending on its width will move at a different veloc-ity inside the cavity, allowing the experimentalists to find an excellent correlation between structure size and velocity [Gustave 2017].

In Fig. 2.26 we depicted the phase portrait (b) and intensity temporal profile (c) of a double-charge PS observed experimentally, which is stable during more that 2×10 4 roundtrips. In (b,c) the color scale illustrates the number of observations, hence the yellow points in (b) are the most visited ones by the trajectory of the system in the phase space. In the case of PSs these points identify the approximate location of the fixed point A. The different position of this fixed point with respect to the previous numerical (in Fig. 2.15) and experimental (in Fig. 2.16) results is simply due to the fact that the experimental phase is actually the phase difference between the injection and the laser relative to an arbitrary reference phase, hence the trajectory in the phase portrait could be arbitrary rotated around the origin without any change in shape and still be consistent with the measurement. One interesting feature of the PS trajectory, that we will address later in this Subsection, is that during a single roundtrip the field does not reapproach the stationary solution in A until the end of the second rotation. The intensity profile shown in Fig. 2.26(c) confirms the shape preservation property of the PS, with basically no fluctuations appearing during the propagation over 2×10 4 roundtrips. In Fig. 2.27 we can find a confirmation of the numerical results on the PS stability domain illustrated in Fig. 2.14: in particular we can observe that PSs of double (a,c) and triple (b,d) charge do not seem to necessarily include the origin in all their rotations as it was already observed numerically for a single-charge PS in Fig. 2.14(d). To better illustrate this point, in Fig. 2.28 we plotted the trajectory of the system in the full 3D phase space (a,b) and its projection on the (Re(E), D) plane during a single roundtrip in the case of a phase soliton of charge 2 (a-c) and 3 (d-f). As we already observed in the case of a single-charge PS in Fig. 2.15, the 2π phase rotations occur at a certain distance from the point C, whose main role in the dynamics is of repulsive nature. Furthermore we can notice that n -1 of the phase rotations performed by the trajectory lie almost on the same plane of the phase space, set for high values of carrier density: the only rotation that does not lie in such area is the one that connects the bounded trajectory to the point A and gets repelled by the point B. Given the inclination of the plane where the n -1 phase rotations tend to lie, it is clear that there might be some values of injection y for which the point C would not appear inside the rotations circle in the projection of the trajectory on the Argand plane: this does not imply that the point C terminates its important role in the dynamics, instead it simply entails that the projection onto the complex plane does not perfectly capture the whole dynamics of the system. Finally from Fig. 2.28(a,b) we can also observe that the reason why, experimentally (in Figs. 2.26 and 2.27) and numerically the trajectory of the system during a single roundtrip does not reapproach the stationary solution in A until the end of the roundtrip, is due to the fact that, when performing the n -1 phase rotations the system is around another value of carrier density, and actually not very close to the point A. (f) Same as (a), in the (Re(E), Im(E)) plane. Fig. 2.28 -Trajectory of the system in the phase space (a,b) and its projection onto the (Re(E), D) (c,d) and (Re(E), D) plane (e,f) for a phase soliton of charge 2 (a-c) and of charge 3 (d-f). A, B and C indicate the three fixed points given a value of injection.

Extreme events from phase soliton collisions

As illustrated in Fig. 2.20, when more complex dynamics is involved, the pure rate-equation model of Eqs. (2.7) is not sufficient anymore to describe the system. In this case the rate-equation model with diffusion of Eqs. (2.8) results more appropriate while still preserving the main characteristics of the model in Eqs. (2.7). The inclusion of the diffusion term can, in general alter a little the instability boundaries, as illustrated in Section 2.3. Nevertheless, as illustrated in Fig. 2.6, when close to the left turning point of the stationary curve, the stability boundary is just slightly affected allowing us to observe most of the features of the complete model in Eqs. (2.5) while preventing the self-collapsing, given by the infinite gain linewidth in the rate-equation model, and avoiding the huge computational time necessary for the simulations based on the complete model. It is through this model that we are able to observe, in the PS stability region, for high pump values, a regime where collisions give rise to extreme events.

Experimental results

In Figs. 2.29 and 2.30 we illustrated the main experimental results for this kind of regime, achieved by P. Walczak in the experimental team with the Fabry-Perot setup illustrated in Fig. 2.11. In particular in Fig. 2.29(a) we depicted the spatiotemporal diagram relative to this specific regime. Initially (for the first ≈ 7.7 × 10 3 roundtrips) the system results in the homogeneous stationary state, then, a different regime occurs, where double charge phase solitons are spontaneously nucleating from the homogeneous solution and colliding with other structures propagating in the opposite direction in this reference frame (which does not imply an opposite propagation direction in the laboratory frame). Around roundtrip 10 4 this regime ceases to exist, probably due to some change in the experimental parameters (most probably temperature fluctuations leading to changes in the detuning parameters), and the only object remaining is a stable double-charge PS which propagates in the cavity without any collision. In Fig. 2.29(b) we reported a zoom of Fig. 2.29(a) on the highest-peak event observed (highlighted by the black dashed lines in (a)), where we can observe more clearly the collision between the dissipative PS and a different transient structure. After the collision, the PS complex seems to loose one charge, which is then recovered around ten roundtrips later. The same PS incurs then in some further collisions with other counterpropagating structures, and, after moving past the collisional regime, propagates in a stable manner inside the cavity. The collisional regime observed induces extreme events in the system as can be observed in the experimental PDF of the intensity illustrated in Fig. 2.30. Since during this regime we can have, at the same time, coherently locked and dynamically complex states, the best approach to understand the behavior of the PDF consists in isolating, where possible, these different states, and identifying which part of the PDF function they are building: in this way we can also use an alternative definition for the extreme events present in the system as those events escaping the statistics of a more simple bounded dynamics. In particular in Fig. 2.30 the black solid line represents the total intensity PDF on the whole spatiotemporal diagram represented in Fig. 2.29. This PDF results asymmetric (with a kurtosis of 19.95) and shows the presence of a heavy tail. The dashed red line identifies the stationary homogeneous state, where the PDF presents a mostly Gaussian shape (with a kurtosis of 2.95): the peak is at the locked intensity value and the data are spread around this value due to physical and detection noise. The blue line with square markers isolates the PDF of the stable PS regime, from around roundtrip 10 4 (with a kurtosis of 5.53), while the green line with circle markers is relative to the PDF of the collisional regime and presents a kurtosis of 11.1. For the definition of kurtosis please refer to the previous Chapter. We can notice that the emergence of extreme events in the collisional regime strongly modifies the PDF, deviating from the Gaussian statistics with an heavy tail. Furthermore we can observe that the first bump of the PDF, at low values of intensity is clearly related to the PS shape: in particular, given the PS shape, with a depression on the trailing edge (see for example Fig. 2.26), the probability of low intensity values is increased. The black dashed line and the yellow dash-dotted line are two rogue wave thresholds, computed as the average intensity plus eight times the standard deviation. The black dashed line in particular is computed on the whole intensity time trace while the yellow dash-dotted line is limited to the data in the collisional regime. These two thresholds confirm the collisional events in the tail of the PDF to be extreme.

Numerical simulations

In the simulations of the model in Eq. (2.8) we kept almost the same parameters as in the previous Section, fixing α = 3, T = 0.3 and θ = -2.7. We instead introduced a diffusion coefficient d = 10 -6 , and set the pump current well above threshold µ = 2 (µ th = 1) and b = 0.0001. The detuning θ = -2.7 and the injection amplitude y = 0.11 were chosen for the homogeneous stationary solution to be triple-valued, where the upper state is stable, thus satisfying the conditions for the presence of PSs in the system. Furthermore the parameter choices for µ, θ and y result very close to the experimental values. The stationary curve and full instability domains for this choice of parameters are illustrated in Section 2.3 respectively in Fig. 2.6(b) and (d).

In Fig. 2.31 we show an example of the kind of collisions observed numerically in a zoom of the spatiotemporal diagram. We can notice that the dynamics here presented shows some clear similarities with the experimental data. In particular we can observe the presence of phase solitons and the emergence of some peculiar transient structures propagating at a different velocity that eventually collide with the PSs, giving rise to an event of high intensity. The initial condition of this simulation was given by the homogeneous stationary solution with a superimposed phase kink of 4π in the electric field (but the same kind of regime can develop also starting from a random initial condition). Lots of collisions can occur in this regime and we can notice that the kind of dynamics generated by a collision can be rather complex, giving rise to three or four peaks as a result of the interaction. In the three collision events represented here the number of charges carried by the PS complex are conserved also after the collision but we would like to point out that this is not always the case. In this work we will not compute the PDF of the total intensity since, even if the regime is rather similar to the experimental one, the larger number of PS complexes inside the cavity and the more complex behavior observed in the simulations, modify the statistics to a point that any comparison with the experimental data would be pointless. Furthermore in the experiment it was possible, in the same acquisition to isolate the different regimes responsible for the whole shape of the PDF: such technique is not available in the numerical case since the collisional regime is not spontaneously alternating with other regimes, as in the experiment, and becomes even more tangled for larger times.

The main focus of this study consists in assessing the nature of the counterpropagating structure that gives rise to collisions in the cavity. To this purpose in Fig. 2.32(a) we show a zoom of the spatiotemporal diagram during the same simulation of Fig. 2.31. In Fig. 2.32(b) we plotted the trajectory of the system in the Argand plane, relative to the horizontal cuts highlighted in Fig. 2.32(a), with time growing from blue to yellow. The points A, B and C represent the usual three fixed points on the homogeneous stationary solution for a given value on injection amplitude. In the first frame, relative to roundtrip section 1, we can observe a phase soliton complex of charge 3 propagating inside the cavity: the trajectory of the system consists in three counterclockwise phase rotations in the Argand plane, passing only once close to point A, as already illustrated for PS complexes in Fig. 2.28. On roundtrip section 2 a new object emerges on the right side of the PS, together with a clockwise phase rotation (of charge -1). The interaction between these two structures gives rise to a high intensity event on roundtrip section 3, where the clockwise rotation has already been lost to a positive charge rotation. We can notice, right after roundtrip section 3 that the interaction has also the effect of altering the phase soliton complex velocity inside the cavity. On roundtrip section 4 the PS complex has regained its shape, with one additional charge, coming from the interaction.

From this phase description we can argue that the appearance of these transient pulses with clockwise phase rotations, that collide with stable phase solitons, is the basic physical mechanism responsible for the extreme events observed in this system.

High-peak events in unstable roll regime

This second main section of the results focuses on the study of abnormally high events in an unstable roll regime. This kind of events has been observed experimentally by F. Gustave in the experimental team in [Gustave 2016a], with the experimental setup illustrated in Fig. 2.10, and is considered particularly interesting due to its connection to the phase dynamics. The model we used to describe this regime is, as in the previous Section, the rate-equation model with diffusion illustrated in Eqs. (2.8).

Roll patterns

The parametric region here considered is chosen so that the stationary homogeneous solution results unstable and a nonstationary irregular roll pattern develops instead. This regime differs from the one reported in the previous Section and in [Walczak 2017], since, there, the upper branch of the stationary curve was stable and PSs were stable as well: the extreme events developed in such regime were because of collisions between localized structures of different nature and velocity while here high-peak events emerge due to the interplay of opposite chiral charges with no clear collisions.

In the numerical results here illustrated we set as in the previous Section the parameters α = 3, T = 0.3 and b = 0.0001. We also fix the pump value to µ = 1.1 and vary the detuning θ and the injection amplitude Y. Setting θ = -3.04 we can observe from Fig. 2.4(c) that the upper branch results stable from values of stationary homogeneous intensity I 0.25. To obtain a roll pattern we set as initial condition the stable homogeneous solution at Y = 0.3 and then decrease the injection amplitude to lower values. Note that to obtain a roll regime some random initial condition of small amplitude has to be introduced to jump on the rolls branch, such condition has no particular importance and is not necessary in the further simulations without affecting the system dynamics. As an example at Y = 0.1 we obtain the spatiotemporal diagram of the electric field intensity depicted in Fig. 2.33(a). Rolls are generated as a result of a multimode instability, in particular when the system is not locked anymore and there is beating between the fundamental mode and other high-order side- modes. The trajectory in the phase space for such regime describes a closed loop around the fixed point A in a bounded region, which does not include the origin, as illustrated in Figs. 2.33(b) and (c), respectively in the phase space and in the Argand plane.

In Fig. 2.34 we depicted the rolls stability region (green dashed line with upward and downward triangle markers) together with the stationary homogeneous solution for the choice of parameters made. In particular the orange and green markers identify the roundtrip minima and maxima of the intensity time trace over 1000 τ c . If all the maxima (and all the minima) converge in one point of the plot for a given value of injection, then the roll pattern is stable: in particular the rolls stability region corresponds here to 0.006 ≤ Y ≤ 0.2. If instead the maxima and the minima spread along the vertical axis it means that the pattern is unstable: this is what occurs for Y = 0.005 and lower values of injection in Fig. 2.34(b), where we illustrated the unstable roll regime with orange markers. This region of more complex dynamics is where we are going to search for abnormally high events. We can observe that rolls start to emerge near the tip of the instability region in Fig. 2.4(b).

In Fig. 2.35 we illustrated more into details the roll pattern dynamics in the phase space and its eventual breaking for low enough values of injection. On the first column we represented the spatiotemporal diagrams relative to simulations for different values of Y. Each row corresponds to a new simulation, in particular we depicted here the final 5000 roundtrips of 10000 τ c -long simulations to be sure of removing any transient effect due to the initial change in the injection value. In the second and third columns we depicted the system trajectory during the last roundtrip of the simulation, respectively in the Argand plane and in the (D, I) plane. The points A, B and C are the usual fixed points on the stationary homogeneous solution for the given value of injection. The right turning point in the stationary curve, illustrated in Fig. 2.5, is, for these parameters, at Y = 0.0307 while the left turning point is at Y = 0.0037. Figs. showing every once in a while some large excursion (yellow line in (h)) around the origin in the Argand plane.

Finally below the (left) turning point the system describes some sort of bounded trajectory around the origin. We would like to stress that, for low values of stationary homogeneous intensity, the model in Eqs. (2.8) becomes less accurate since it starts to present bigger differences with the complete model in Eqs. (2.5), especially in terms of the number of sidemodes involved in the dynamics (see Fig. 2.4), hence any consideration on this final regime is to be considered in its limits.

From Fig. 2.35 it results clear that the most promising regime to observe high-peak events is represented by the simulation depicted in (g-i), when the system results very close to the left turning point, the three fixed points A, B and C are still present, affecting the dynamics, and the roll regime becomes unstable. This condition results very similar to the experimental case, where abnormally high events on a roll background have been observed, as illustrated in Fig. 2.36. Additionally to Fig. 2.5, in Fig. 2.37 we depicted the stationary homogeneous solution of the system as in Eq. ( 2.15) for the intensity as a function of the detuning θ, for the choice of parameters made: the vertical red dash-dotted line identifies the value of θ = -3.04. This plot gives us an idea of the values of detuning that still grant the presence of the three fixed points on the homogeneous stationary solution and it will be useful in the statistical analysis.

Spectral analysis

In order to characterize better the unstable roll regime observed in Fig. 2.35(g), in Fig. 2.38 we make a comparison between the Fourier spectrum for a simulation at θ = -3.03, where the roll pattern is stable, and a simulation at θ = -3.04 (same as in Fig. 2.35). In Fig. 2.38(a) and (b), on a zoom of the frequency-time diagram, we can notice that, while in the case (a) of the stable roll pattern the Fourier spectrum remains the same for different roundtrips, in the case (b) of unstable rolls it shows a more complex behavior. In Fig. 2.38(c) we show both spectra averaged on 5 × 10 4 roundtrips: here we can observe that the spectrum in the case θ = -3.04 results more continuous but not particularly broader with respect to the roll regime. For this reason we do not make any claim on the turbulent nature of the unstable roll pattern here studied.

High-peak events and statistical analysis

In Fig. 2.39(a) we depicted a zoom of the spatiotemporal diagram on 2000 roundtrips in the middle of a simulation, run for the choice of parameters explained in the previous Subsection (in particular pump parameter µ = 1.1, detuning θ = -3.04 and injection amplitude Y = 0.004). The total simulation covers a time span of 3 × 10 5 roundtrips where the first 9 × 10 4 roundtrips were excluded from the analysis to discharge any transient effect.

In Fig. 2.39(a) an abnormally high event, highlighted by line 3, emerges above an unstable roll background. In Fig. 2.39(b) we show the four intensity temporal profiles corresponding to the horizontal dashed lines illustrated in Fig. 2.39(a). Section 1 of the spatiotemporal diagram is chosen 500 roundtrips before the event: the temporal profile displays a mostly bounded roll pattern regime, which would be stable for higher values of injection. Section 3 is chosen on the roundtrip where the maximum of the high-peak event occurs. Finally sections 2 and 4 have been chosen 40 roundtrips before the event and 160 roundtrips after it.

As in the previous Chapter, to study this kind of abnormal events, we find useful to perform a statistical analysis on the intensity values explored by the system, as well as on the heights of all the observed intensity peaks. The results of this study are illustrated in Fig. 2.40: in particular in Fig. 2.40(a) we illustrated the PDF (in black) of the peak heights, measured as the difference between each maximum and the previous minimum of the full time trace, while in Fig. 2.40(b) we depicted (in black) the PDF of all the values explored by the intensity during the simulation.

The dark yellow data in Fig. 2.40 reports the results of the same statisti- cal analysis performed on a simulation with for a different value of detuning θ = -3.048: in this case the system exhibits an unstable roll regime, which remains bounded to a specific region of the phase space and does not show big intensity excursions as illustrated in the spatiotemporal diagram in Fig. 2.42(e). We can observe that the PDF (in black) of the peak heights (a) displays an initial slope, following a Gaussian distribution (green dashed line), which is to be associated with the bounded regime. The events in the tail of the distribution are escaping such statistics and occur more frequently than expected by it. In Fig. 2.40(b), for the statistics of the total intensity (in black), we can notice an initial almost flat distribution (in logarithmic scale), up to I/ I ≈ 3, followed at first by an exponential decay (fitted by the green dashed line): such behavior is closely followed by the dark yellow PDF of the simulation without high peaks. The data in the tail of the black PDF seems once again to escape from the negative exponential decay, displaying a higher probability than expected by the bounded regime. In Fig. 2.41 we show the same data analysis performed in the experimental case on a 27000 τ c -long time trace of the data displayed in Fig. 2.36: these results are in very good agreement with the numerical ones, showing once again the similarities between the experimental regime and the simulation. In the experimental case, for the PDF of the peak heights in Fig. 2.41(a) it was necessary to disregard all the events of height below 5, since for low values of intensity there was no way to distinguish the field dynamics from the detection noise: this choice avoids the use of any spectral filtering on the data which could easily modify the shape of the distribution. For this reason also in the numerical case we decided to put a threshold on the data with low values of intensity in Fig. 2.40(a), to grant a better comparison with the experiment.

From the (nontrivial) shape of the PDFs here shown, we can easily understand that the events in the tail of the distributions are generated by a different mechanism with respect to the rest of the data. For the sake of completeness in Fig. 2.42 we show the spatiotemporal diagram and PDF of peak heights for three other values of detuning, in a configuration of unstable rolls.The statistics on the heights has been computed here over 3×10 5 roundtrips. The value of θ = -3.04, chosen in the analysis, results the one producing the most high-peak events, hence the most adequate for this study.

Phase dynamics

To better understand the underlying mechanism for the formation of these abnormal high-peak events we will now focus on the role of the phase dynamics. To this aim, it is useful to set E = ρ exp (iφ) and recast Eqs. (2.8) for the field where we have neglected the diffusion term d to keep the equations in a simpler form and allow some qualitative considerations (d is any way present in all the simulations here shown).

+ ∂ρ ∂t = T [(D -1) ρ + y cos φ] (2.24a) ∂φ ∂z + ∂φ ∂t = -T θ + αD + y ρ sin φ (2.24b) ∂D ∂t = bT σ µ -D 1 + ρ 2 , ( 2.24c)
In Eq. ( 2.24b), we can observe that, for fixed z or t, the sign of the right hand side determines the rotation direction of the electric field in the complex plane.

For θ + αD = 0 and constant ρ, Eq. ( 2.24b) becomes that of a pendulum, with a stable fixed point in φ = 0 and an unstable one in φ = π. Nevertheless the unstable focus C lying in the proximity of the origin, acts as a repeller and prevents the field amplitude to assume values too close to the origin, hence the ratio y/ρ in front of the sine remains usually very small. The dominant term in Eq. (2.24b) is then the first one, which forces the electric field to rotate clockwise or counterclockwise depending on its sign. The boundary between the two rotation directions is marked by the critical value of D, D c = -θ/α , in particular with our choice of the parameters D c ≈ 1.013. If the dynamics of D were fast enough to follow adiabatically that of the field intensity, we could also define a critical amplitude through the stationary equation .25) in particular with our parameters we would have ρ c ≈ 0.292. In Fig. 2.43(a) and (b) we can observe that the temporal evolution of D and its stationary homogeneous value D s = µ/(1 + I) are not comparable in amplitude, hence enslaving D to the evolution of I is here not really justifiable. In particular, using a loose terminology we could say that a "drag" linked to the medium dynamics prevents the system from the large excursions in D predicted by the instantaneous model (where D adiabatically follows the dynamics of I). Nevertheless the critical amplitude ρ c plays a role in the phase dynamics, as will be shown in the following.

D c = µ 1 + ρ 2 c , ( 2
Re(E) (a.u.)

Im(E) (a.u.) events In Fig. 2.44 we draw a comparison between the histogram of the values assumed by the electric field in the complex plane (Re(E), Im(E)), both in the experimental (a) and numerical case (b). There are some clear similarities between the two figures, especially regarding the presence of an unstable bounded state (yellow area) from which rarer events get detached, exploring other areas of the histogram (depicted in red). As already explained in Section 2.5, the different location of the bounded state in the experimental case is due to the fact that the relative phase is measured up to a constant.

Chirality plays an important role in the formation of high-peak events both in the experimental and in the numerical case. In Fig. 2.45 we show the evolution of the unfolded phase in the experiment around the formation of an abnormal event: in this description each step depicted in Fig. 2.45 represents a 2π rotation of the system trajectory around the origin. Here we can observe that the phase does not perform any jump up to 20 ns before the abnormal event (arbitrary placed at t = 0), then a change in the phase slope identifies the acquisition by the system of a positive chiral charge and shortly after a high-peak event occurs.

In Fig. 2.46 we illustrated a zoom of the spatiotemporal diagram on an abnormally high event in the experiment (a), and four different intensity temporal profiles (b) acquired at fixed roundtrip, corresponding to the cuts highlighted in (a). Temporal profile II is acquired when the abnormal event assumes its highest value in intensity, temporal profile I corresponds to 43 roundtrips before the event, while temporal profiles III and IV are collected respectively 12 and 37 roundtrips after the event. The main experimental result about the phase dynamics is summarized in a striking topological change for the trajectory of the system in the phase space: as we can observe in Fig. 2.47, before the abnormal event, on roundtrip I, the system trajectory in the complex plane remains bounded to a specific region while when the high-peak event occurs, on roundtrip II, a 2π excursion of the phase takes place with the trajectory of the system describing an excursion around the origin in the complex plane.

In Fig. 2.48(a) we illustrated the evolution of the unfolded phase when approaching an abnormally high event in the simulations, in a situation similar to the experimental case. Here we can observe that an initial balance between the number of positive and negative chiral charges is broken by the loss of a negative charge, which gives rise to a change in the slope of the phase: then an abnormal event, highlighted by the vertical red dashed line, occurs. This description finds confirmation in the zoom of Fig. 2.48(b) around the slope change: here we can distinguish the presence of both a positive and a negative chiral charge in the initial plateau, which is then reduced to a single positive charge with the change of slope.

Figs. 2.45 and 2.48 present similar results especially when discussing the occurrence of the high-peak events in terms of a change of slope in the phase evolution from an initially balanced (flat) situation. About this point we would like to stress out the limits of the experimental results regarding the phase rotation detection: as we will see in the following, both positive and negative chiral charges (i.e. respectively counterclockwise and clockwise phase rotations) are necessary to well describe the dynamics of the high-peak events observed and eventually allow us to suggest a predictor for these abnormal events. Unfortunately the only clear observations in the experiment always regard the positive chiral charges: this is because clockwise rotations are instead bound to occur for low intensity values, beyond the experimental noise level. Furthermore the array of possible configurations regarding the evolution of the phase before or after the occurrence of an abnormal high-peak event is more complex as well illustrated in Fig. 2.49, where we depicted (a) the spatiotemporal diagram and (b) the corresponding unfolded phase time trace during 2×10 4 τ c of the simulation. The five abnormal events occurring in this time span are highlighted with red crosses in (a) and vertical red dashed lines in (b).

In general we can argue that, as a common thread, the occurrence of abnormal events may be associated with a global change in the slope of the phase: however, as it has been observed in the experiment, this change is not to be considered as a precursor since it can occur either before or after the abnormal event itself. In Fig. 2.49(b) we can also notice that, even if sometimes in a single roundtrip, the number of clockwise rotations exceeds the number of counterclockwise rotations, leading to a negative slope for the phase, the number of positive chiral charges largely exceeds the number of negative ones: this strongly asymmetric behavior, observed also experimentally [START_REF][END_REF]], is to be related to the propagative nature of our system.

In order to go more into details in the phase description of the system, it is interesting to analyze the trajectories in the phase space projections (Re(E), Im(E)) and (D, I) to point out any phase rotation peculiarities. In Fig. 2.50 on columns (a) and (b) we show the trajectories of the four horizontal sections highlighted in Fig. 2.39, where the line color changes from dark green to yellow for increasing time (from the left to the right side of each section). The (unstable) locked stationary states A, B and C of Fig. Im(E) (a.u.)

Re(E) (a.u.)

-1.5 -1 -0.5 0 0.5 1 1.5 ×10 space, where the red plane crossing the phase space identifies the critical D c . The fixed points A, B and C have been drawn every 2π in φ since their phase values are defined unless a 2π-modulus.

In Fig. 2.50 at first we can observe in (a.1)-(c.1) a mostly bounded regime, where the trajectory of the system remains in a specific area, also highlighted in Fig. 2.44, of the complex plane (a.1) without performing any full phase rotations. Then, 40 roundtrips before the event, in Fig. 2.50(a.2)-(c.2) as the trajectory of the system starts to move towards higher values of the population D, we can observe the occurrence of a first clockwise phase rotation (blue circle → blue cross) and a second counterclockwise rotation (red circle → red cross). Here we can notice that:

• the clockwise rotation is triggered for very low (almost zero) values of intensity in the proximity of the condition D = D c ;

• after performing a full -2π rotation (ending on the blue cross), the system reaches a maximum value of D and is then repelled to high values of the intensity I due to the repulsive nature of the fixed point C;

• due to a favorable φ the term in sin φ in Eq. ( 2.24b) becomes more and more important in the determination of the sign and, at the point highlighted by the red circle in (b.2), a counterclockwise phase rotation starts forming;

• such rotation is then strengthen when the trajectory approaches the red vertical dashed line, below which the first term of the phase equation becomes positive as well, reaching 2π at the point highlighted by the red cross.

The trajectory excursion toward high values of intensity, due to the interplay between these two opposite phase rotations, reaches its maximum in Fig. 2.50(a.3)-(c.3) and, after that, it starts decreasing, as depicted in Fig. 2.50(a.4)-(c.4), where the two phase rotations are still present (as well indicated in (c)) but the trajectory is slowly moving towards lower values of population, eventually evolving into a bounded state as in Fig. 2.50(a.1)-(c.1), where phase rotations will disappear (usually not at the same time).

In Fig. 2.50(c) the positive and negative phase rotations are particularly highlighted and we can well distinguish the difference between the mostly bounded state of the first frame with the rest of the dynamics depicted. Please note that in Fig. 2.50(c) the relative phase variable φ has been arbitrarily set to zero at the beginning of the trajectory in the first frame, in order to better identify the two main phase jumps. Furthermore we would like to point out that the phase rotation starting points (red and blue circles) are completely determined since the algorithm used is set to maximize the number of phase rotations encountered. In particular, analyzing the unfolded phase time trace (as in Fig. 2.49(b)), we identify all the maxima and all the minima. Starting, for the counterclockwise (clockwise) phase rotations, from the first minimum (maximum) of the phase time trace, we search for the first point with a phase value exceeding of 2π the phase value of the first minimum (maximum), which identifies the first positive (negative) phase rotation. For the next phase rotation we start from the last point of the previous rotation but if, along the phase time trace, we encounter a new lower minimum (higher maximum) we then set the starting point for the phase rotation on such spot.

The phase and amplitude dynamics preceding the emergence of these highpeak events are reminiscent of the defect-mediated turbulence discussed theo-retically in the context of forced oscillatory media [Coullet 1989, Gibson 2016].

The mechanism illustrated in Fig. 2.50 appears common to all the high-peak events in the tail of the distribution in Fig. 2.40(a) and could serve as a predictor for abnormal events. However a full analysis of the simulation results for other parameter choices remains to be performed in order to confirm this statement. In Fig. 2.51 we consider a zoom of the spatiotemporal diagram for the intensity The mechanism already illustrated in Fig. 2.50 is here even clearer: at first a negative chiral charge (in light blue) appears; then, around 100 roundtrips before the event, a positive chiral charge (highlighted in orange) is formed. The simultaneous presence of a clockwise and a counterclockwise phase rotations is needed to observe the abnormal event. Finally, both charges disappear, at different times. We would like to highlight, from Figs. 2.50(b) and 2.51(b), that clockwise rotations always occur only for high values of population, when the intensity of the electric field approaches zero, and it is indeed the interplay of clockwise and counterclockwise phase rotations that gives rise to the observed abnormal events and could serve as a predictor for such abnormal chiral events. In fact, it may happen that clockwise rotations occur without any positive chiral charge in their proximity: such a configuration does not give rise to high-peak events.

This interplay between ±2π phase rotations results completely absent in other parameter choices, for simulations that do not exhibit big excursions in intensity, as the one of θ = -3.048 illustrated in Fig. 2.40.

As a last remark we would like to point out that the observed phase rotations may be associated with a transient exploration of the PS regime, unstable for this choice of parameters since, in a stable PS regime, counterclockwise rotations would correspond to stable PSs. However clockwise rotations cannot be associated with unstable PSs due to the asymmetric nature of this propagative system [START_REF][END_REF]], since, as anticipated in Section 2.2.5, the non-instantaneous medium dynamics heavily breaks the η → -η symmetry of the system.

Conclusions

In this Chapter we have studied the dynamics of a semiconductor laser with injection, spatially extended along the propagation direction. The model used to describe this kind of laser has been carefully derived and compared with the different models also used for the description of this system. The three main regimes analyzed highlighted different aspects of the possible interactions between structures in this system. In particular we studied into details the attractive interaction between PSs, which leads to the possibility of phase soliton complexes in the PS regime, observed also experimentally. Interestingly the merging time of two separate structures appears in this case to depend exponentially from their initial distance, a fact that will be recalled also in the next Chapter. We then studied high-peak events in this model for two different configurations. In a first configuration, for high enough pump values in the PS regime, collisions between phase solitons and some transient structures leads to the emergence of extreme events in the system. In a second configuration abnormally high peaks are developed in an unstable roll pattern regime, where PSs are not a stable solution. Both these configurations present an excellent agreement with the experimental data, which has been corroborated, in particular for the last case by a detailed statistical analysis: here we observed that the abnormally high events in the tail of the distribution are clearly distinguishable from the mostly phase-bounded background, both by their size and their nature.

In all three regimes the chirality of the system plays a fundamental role. The observed phase and amplitude dynamics is strongly reminiscent of the regime of defect mediated turbulence [Coullet 1989] or vortex turbulence that has been recently associated with optical rogue waves in an oscillatory two-dimensional transverse spatial system with coherent forcing [Gibson 2016]. In our case, while counterclockwise phase rotations are intrinsically associated to phase solitons or to transient explorations of such regime, clockwise phase rotations, always unstable due to the propagative nature of the system, seem to be a necessary tool to observe the emergence of high-peak events, when combined with positive chiral charges.

Chapter 3

Broad-area semiconductor laser with injection Preface

This Chapter is devoted to the study of a broad-area vertical-cavity semiconductor laser (VCSEL) with optical injection. After a brief introduction on the relevance of this system in the context of cavity solitons studies, we are going to compare two different models for the description of the system. We will then focus on the study of the interaction between two cavity solitons in the transverse plane. In particular we will observe that cavity solitons merge for any given initial distance with a merging time displaying an exponential dependence on such distance: this allowed us to associate an exponentially decaying interaction potential to cavity solitons and permits for a possible analogy with hydrophobic materials. The results here presented have been submitted to Physical Review E [Anbardan 2017]. Furthermore we report a preliminary analysis on the presence of extreme events in this system and compare the results with [Gibson 2016].

Introduction

Broad-area semiconductor lasers, as already explained in the Introduction of this Thesis, represent a privileged framework for the study of cavity solitons, given their fast response and compactness, which makes them good candidates in the context of optical information processing [Firth 1996, Brambilla 1997[START_REF] Ackemann | [END_REF]. Spatial cavity solitons were first predicted in the transverse plane of semiconductor cavities with optical injection in [Brambilla 1997[START_REF] Michaelis | [END_REF], Spinelli 1998]. VCSELs with optical injection below threshold were the first experimental setup where cavity solitons in semiconductor microcavities were observed [Barland 2002]. In this kind of setup CSs can be independently manipulated, that is written and erased, through the use of an external optical perturbation.

Interaction between CSs is a relevant subject in the context of optical information encoding: in fact the interaction range of these structures defines the maximum information density of an optical memory array (as for example the one in [START_REF] Pedaci | [END_REF]).

In the case of temporal cavity solitons, a critical distance, above which CSs stop interacting, has been found in [Leo 2010], corresponding to 40 ps, leading to a maximum information storage density of 125 bit/m. Regarding spatial localized structures, in the early studies of [Rosanov 1990, Brambilla 1996], about nonlinear media as bistable interferometers or collections of two-level atoms, two critical distances have been observed: a first one, d 1 , below which CSs would merge or annihilate, and a second one, d 2 , above which CSs would not interact. For an initial distance d in between d 1 and d 2 the two CSs have been found to repel until the distance d 2 was reached. More specifically in the case of a VCSEL with optical injection pumped below transparency, the presence of the same critical distances was observed in [Tissoni 1999]. Furthermore some equilibrium distances d * were found, such that the interaction between the two CSs would become either attractive or repulsive in order to reach a distance d = d * between the two structures.

Model

Let us consider a broad-area VCSEL with optical injection. This kind of system is well described by the effective Maxwell Bloch equations [Hachair 2006]

Ė = σ E I -(1 + iθ) E + P + i∇ 2 ⊥ E , (3.1a) Ṗ = ξ(D) [(1 -iα) f (D)E -P] , (3.1b) Ḋ = µ -D -(E * P + EP * ) /2 + d∇ 2 ⊥ D , (3.1c)
where E and P are, respectively, the slowly varying electric field and effective material polarization and D is the carrier density. E I is the injected field amplitude, θ represents the detuning between cavity and optical injection, α is the linewidth enhancement factor, µ represents the injection current and σ is the ratio of carriers lifetime τ c to photon lifetime τ ph . Differently from the model introduced in the previous Chapter a diffusive term is introduced in the equation for the carrier density with d. Diffraction is instead described through the complex term in the transverse Laplacian ∇ 2 ⊥ in the equation for E. Time is scaled to τ c (≈ 1 ns), while the spatial scale corresponds to the square root of the diffraction parameter (here equal to 1). The gain nonlinearity, typical for quantum wells [Eslami 2014], is taken into account by the real function f (D). In particular, assuming a quadratic nonlinearity, we have f (D) = (1 -βD)D with the value of β = 0.125 resulting from the best fit of the gain in the microscopic model [Eslami 2014]. The free running laser threshold corresponds then to µ thr = (1 -1 -4β)/(2β) and specifically for β = 0.125 the threshold current results µ thr = 1.17. The complex function ξ(D) in the equation for the polarization P accounts for the dependence of the gain linewidth and the maximum gain frequency on the carrier density [Hachair 2006] similarly to the case of Chapter 2. When θ + α is different from zero, for E I smaller than the injection locking point, there are no stable solutions, the laser output is nonstationary and many transverse modes can be excited. In these conditions the polarization P acts as a spectral filter, since it includes the effects of the finite gain linewidth, thus avoiding nonphysical short wavelength instabilities. Since in the present case we are focused on the interaction of CSs occurring beyond the injection locking point, it is possible to reduce Eqs. (3.1) through the adiabatic elimination of P [Prati 2010b] without the additional introduction of a diffusive term and still obtain a model well describing the system under study, given by the following set of equations

Ė = σ E I -(1 + iθ) E + (1 -iα) f (D)E + i∇ 2 ⊥ E , (3.2a) Ḋ = µ -D -f (D)|E| 2 + d∇ 2 ⊥ D . (3.2b)
When describing the system solutions in other parameter configurations, such as below the injection locking point (for example in the study of extreme events), it will instead be necessary to consider a reduced model of the following kind

Ė = σ E I -(1 + iθ) E + (1 -iα) DE + (d + i) ∇ 2 ⊥ E (3.3a) Ḋ = µ -D 1 + |E| 2 (3.3b)
where d is the usual diffusive term, as the one introduced in Chapter 1.

HSS and linear stability analysis

The linear stability analysis of the complete model in Eqs. (3.1) for β = 0 is reported in [Prati 2010b].

The reduced model in Eq. (3.3) admits the following homogeneous stationary solution (HSS)

E s = |E s |e iφ s , D s = µ 1 + |E s | 2 with E 2 I = |E s | 2 (1 -D s ) 2 + (θ + αD s ) 2 (3.4) φ s = arctan θ + αD s D s -1 ,
similarly to what was found for the semiconductor ring laser with optical injection in Chapter 2.

The model in Eqs. (3.2) admits instead the following HSS

E s = |E s |e iφ s , D s = 1 + |E s | 2 -(1 + |E s | 2 ) 2 -4βµ|E s | 2 2β|E s | 2 with E 2 I = |E s | 2 [1 -f (D s )] 2 + [θ + α f (D s )] 2 (3.5) φ s = arctan θ + α f (D s ) f (D s ) -1 .
In particular we can observe that for a Taylor expansion at the first order in β the stationary homogeneous solution for the carrier density is approximated by

D s ≈ µ 1 + |E s | 2 1 + βµ|E s | 2 (1 + |E s | 2 ) 2 .
In Fig. 3.1 we illustrated the homogeneous stationary curve as in Eq. (3.5) (blue solid line) and in Eq. (3.4) (red solid line) for the parameters α = 4, β = 0.125, θ = -2 and µ = 1.2µ thr (where µ thr = 1 for the model in Eqs. (3.3)). The vertical black dashed line indicates a chosen injection amplitude E I = 1 for the CS interaction study. Although some differences in the shape of the two curves can be observed we can notice that the lower branches of the two solutions are almost identical.

We now study the stability of the HSS of the two reduced models for small spatiotemporal perturbations from the stationary values of E s , E * s and D, modulated in the transverse plane by a wavevector k = (k x , k y ) and exponentially growing (decaying) in time so that: where J (3.2) is the Jacobian matrix

  E E * D   =   E s E * s D s   + e λt+i(k x x+k y y)   δE δE * δD   . ( 3
  -σ 1 -α f (D s ) + i k 2 + θ 0 σ α(1 -2βD s )E s 0 -σ 1 -α * f (D s ) -i k 2 + θ σ α * (1 -2βD s )E * s -f (D s )E * s -f (D s )E s -(1 -2βD s )|E s | 2 -(1 + dk 2 )   ,
with α = 1iα (and α * = 1 + iα).

Similarly, for the model in Eqs. (3.3), the following linearized equations for the perturbations are obtained

λ   δE δE * δD   = J (3.3)   δE δE * δD   where J (3.3) is the Jacobian matrix   -σ 1 -αD s + i k 2 + θ + dk 2 0 σ αE s 0 -σ 1 -α * D s -i k 2 + θ + dk 2 σ α * E * s -D s E * s -D s E s -1 -|E s | 2   .
The From Fig. 3.2 we can observe that the stationary instability domain identifies the negative slope and part of the upper branch of the homogeneous stationary curve as stationary unstable. The Hopf instability domain in Fig. 3.3 displays a peculiar behavior for low values of intensity, due to the absence of a diffusion coefficient in the equation for the electric field: as we have already observed in Chapter 2 (especially in Fig. 2.4) when the adiabatic elimination of the polarization variable is not compensated by a diffusive term in the second spatial derivative of the electric field, no limitation on the number of sidemodes is introduced, leading to nonphysical short wavelength instabilities. Nevertheless, when avoiding the parameter regions where the homogeneous stationary curve presents low values of intensity I 0.23, the model in Eqs. (3.2) does not show any self-collapsing problem. From Figs. 3.2, 3.3 we can observe that, for the chosen value of injection amplitude E I = 1 (which corresponds to a value of intensity I = 0.2825 on the homogeneous stationary curve for Eq. (3.5)), the lower branch of the homogeneous stationary curve of the model in Eqs. (3.2) results stable, while the negative slope and the upper branch are stationary unstable.

With the chosen parameters, the model in Eqs. (3.2) presents stable cavity solitons on a stationary homogeneous background for E I = 1. In particular the radius (HWHM) of a CS in these conditions is around one space unit.

In Fig. 3.4 we depicted the stationary curve (a) and the stationary (b) and Hopf (c) instability domains for the model in Eqs. (3.3) and the parameters α = 4, θ = -2, µ = 6, d = 0.01 and σ = 400. These results will be useful for the investigation of extreme events in the system, to which we dedicated the second part of the present Chapter. In particular we can notice that for this choice of parameters the model in Eqs. (3.3) results stationary unstable in the negative slope branch and part of the upper branch of the stationary curve and Hopf unstable in the lower branch.

CS interaction

Interactions between CSs in a VCSEL with optical injection below transparency and above transparency (but below threshold) have been studied extensively in [Tissoni 1999, Barland 2002]. In the present Section we will show that the interaction of CSs in a VCSEL with optical injection above threshold presents some fundamental differences with the previous case and displays some interesting analogies with hydrodynamics.

The simulations were performed on the model in Eqs. (3.2) as follows: two CSs were switched on at the same time in the transverse plane of the electric field intensity at an initial distance d 0 through the use of two Gaussian input beams, of width (FWHM) equal to one spatial unit, amplitude 1, for the initial 3 time units. After this initial time the two Gaussian inputs are turned off and two stable CSs persist in the transverse plane of the electric field intensity. At this point the attractive interaction between the two CSs starts to come into play as illustrated for example in Fig. 3 periodic boundary conditions. Although the spatial grid dimension of 64 × 64 pixels results slightly small to well describe the spatial structure of CSs, the distance between two CSs is very well evaluated through the use of a simple data interpolation of the spatial profiles as in Fig. 3.5(b,d,f,h). More details about the integration method of the model in Eqs. (3.2) are given in Appendix G.

Even for an initial distance of ≈ 20 spatial units, which is much larger than the CS diameter, the two CSs still experience an attractive force, leading to the merging into a single structure. It is important to observe from Fig. 3.5 that for the first 10 3 ns the two CSs approach each other very slowly, while the final merging occurs in about 6 ns (for these value of σ and d 0 ).

In Fig. 3.6 we displayed the time evolution of the cavity soliton distances for different values of r 0 = d 0 /2: also in this case we can observe that the interaction results very weak in the beginning, in particular for large distances while, in comparison, the final steps show a very strong attraction, leading to a fast merging. From Fig. 3.6 we can also observe that the merging time separation between two consecutive initial distances results almost constant in a logarithmic scale, which suggests that the merging time might depend exponentially on the initial distance. In Fig. 3.7 we plotted the CSs merging times as a function of the initial distance, in a logarithmic scale: although the data relative to small distances slightly deviates from the exponential fit for higher distances, the assumption of an exponential dependence of the merging time on the initial distance results here essentially confirmed. This result is similar to what we have already observed in the case of phase soliton (PS) complexes in Fig. 2.22 in Chapter 2, for the merging of two single-charge PSs, where the dependence of the merging times on the initial distance was also exponential [Gustave 2017]. The main difference between the two cases consists in the fact that, for phase soliton complexes, the merging leads to a different structure where both chiral charges survive in a bounded state, while, in the present case, the CS after the merging is identical to the two previously-existing CSs.

Interaction potential and merging time

The particle-like interaction of dissipative solitons is often described through the use of effective potentials, especially to the aim of finding the equilibrium distance in bound states [Malomed 1993, Malomed 1994, Vladimirov 2002]. Here we aim to extend this idea to the free fall motion of one soliton toward the other. In general the merging time of two masses under gravitational potential U(r) ∝ 1/r scales as r 3/2 0 [Foong 2008]. In the present case, given a merging time that increases exponentially in time with the initial distance, we make the assumption that the two CSs are subject to an interaction potential V(r) exponentially decreasing with the distance between the two objects.

V(r) = -K 2 e -r/R .
where K and R represent respectively the strength and the range of the potential. We assume that the motion of the two CSs, considered as particles under a potential, is conservative. Hence, given an initial distance r 0 between the two bodies, initially at rest, the total system energy is

E = V(r 0 ) = -K 2 e -r 0 /R .
We obtain then for energy conservation, for a particle mass equal to 2, that at a given distance 0

≤ x ≤ r 0 mv 2 (x) 2 -K 2 e -r/R = -K 2 e -r 0 /R v 2 (x) -K 2 e -r/R = -K 2 e -r 0 /R hence v(x) = -K e -x/R -e -r 0 /R ,
where we can observe that the potential strength K represents the impact velocity in the limit r 0 R. Finally the time needed to reach a distance d corresponds to

t(r) = 1 K r 0 r dx v(x) = 2R K e r 0 /(2R) arctan e (r 0 -r)/R -1 . (3.7)
and the merging time t m = t(0) is t m = 2R K e r 0 /(2R) arctan e r 0 /R -1 .

For an initial distance r 0 even just a few times larger than the interaction range R, the arctan term can be replaced by π/2 thus leading to the following approximated exponential law for the merging time t m,ap = π R K e r 0 /(2R) .

(3.8)

In Fig. 3.7 we plotted, in black solid line, a fit of the numerical data obtained discarding the distances r 0 < 10, since Eq. (3.8) results valid only in the limit of large initial distances. The slope of this linear fit (in a logarithmic scale) allows us to obtain a high precision estimation for the potential range R as R = 0.474: this value corresponds approximately to one half the CS radius, hence we can conclude that the merging dynamics of the two CSs can be described as two particles under a short-range potential.

A complete time evolution of the distance d could be obtained by graphically inverting Eq. (3.7) if both R and K were known. Unfortunately the uncertainty on the intercept of the linear fit results much larger than that of the slope, thus precluding a precise estimation of the interaction strength. The analytic curve for the initial distance r 0 = 10 in Fig. 3.8 is plotted for a value K = 48.77, which corresponds to the intercept a = -3.382, a value compatible with the uncertainty of the fit, chosen so that the merging time results the same for the numerical and analytic curves. The analytic curve seems to well describe the exact time dependence of the distance, since the larger deviations observed (in the central stage of the motion) remain smaller than 3%. For instance, at t = 700 t.u., the distance in the numerical simulation is r = 9.493 and the theoretical one is r = 9.235.

In the simulations results illustrated in the previous figures the relative decay rate σ = τ c /τ ph was kept fixed as σ = 400: this corresponds to assuming photons to be 400 times faster than the carriers, which is a typical value for semiconductor microresonators. We studied the influence of the relative decay rate σ on the interaction strength, running simulations with σ = 40, 4 and 1. A smaller value of σ corresponds to a longer (external) resonator. The linear stability analysis results for these values of σ are illustrated in Figs. 3.2 and 3.3, where we can observe that, for the fixed value of injection E I = 1, the fixed point on the lower branch of the homogeneous stationary curve is stable. In Fig 3.9 we illustrated the CSs merging time (in logarithmic scale) as a function of their initial distance for σ = 400, 4 and 1. The results for σ = 40 differ very little from those for σ = 400 and were, for the sake of clarity, excluded from this plot. We can observe that the three linear fits (in black solid line) displayed in Fig. 3.9 present the same slope: the only difference consists in an upward shift of the lines for decreasing values of σ, implying that the CSs merging time increases with the photon lifetime. In conclusion the range R of the interaction potential, related to the slope of the data, seems independent from the choice of σ while the interaction strength K decreases with σ. Nevertheless, given the high uncertainty on the values of K it is not possible to assume a particular functional dependence between K and σ.

In Fig. 3.10 we plotted the CSs merging time dependence on σ -1 for different values of initial distances, which clearly results linear. In general we observed that

t m (r 0 ) = f (r 0 ) 1 + η σ . (3.9)
where f (r 0 ) represents the limit value of the merging time when the photon life- time approaches zero and can be estimated as f (r 0 ) = 0.0321e 1.053r 0 . Comparing with Eq. (3.9) we obtain R ≈ 0.474, as expected, and K ≈ 51.5, which is the maximal interaction strength and also the maximal possible impact velocity. Assuming that one space unit corresponds to 4 µm and one time unit corresponds to 1 ns, the maximal impact velocity amounts to 200 µm/ns = 200 km/s, which is at least one order of magnitude larger than the average velocities measured for a CS set in motion by a phase gradient [Pedaci 2008]. The parameter η in Eq. (3.9) does not seem to depend on r 0 , in particular, from the four fits in Fig. 3.10, we obtain η = 1.7 ± 0.008. The time variable t is here scaled as t = τ/τ c where τ is the physical time and τ c the carrier relaxation rate.

Considering that σ = τ ph /τ c , we can rewrite Eq. (3.9) as follows

τ m (r 0 ) = f (r 0 ) τ c + 1.7τ ph . (3.10)
In conclusion in an optically injected VCSEL above threshold, cavity solitons experience an attractive force even at distances much larger than their diameter, in contrast with what reported for passive semiconductor cavities [Tissoni 1999].

The merging time with good approximation depends exponentially on the initial distance and the CSs merging dynamics can be interpreted as the conservative motion of two particles under an exponentially-decaying potential.

Possible analogy with hydrophobic materials

High-dimensional optical dynamics, arising from the competition of a large number of spatial and/or temporal degrees of freedom, often present several analogies with hydrodynamics [Brambilla 1991, Staliunas 1993]: among these we can list optical vortices [Coullet 1989], turbulence [Arecchi 1991], photon flux along channels [Vaupel 1996], and, as extensively explained in the Introduction of this Thesis, rogue waves [Solli 2007, Onorato 2013, Dudley 2014].

In the case of localized structures there are some striking similarities between the trajectories of confined self-propelled cavity solitons in a laser with saturable absorber [Prati 2011] and those of the so-called "walkers", which are droplets bouncing over a vessel containing the same liquid, vertically vibrating close to the Faraday instability threshold [START_REF] Protière | Particle-wave association on a fluid interface[END_REF]].

An exponentially decaying interaction potential is usually associated to the hydrophobic force [Israelachvili 1982], which is the strong attraction (stronger than the van der Waals force) experienced by nonpolar molecules and surfaces in water [Donaldson 2015]. We can then suggest that two CSs in an optically driven laser above threshold interact in a way analog to hydrophobic materials. These results are once again in contrast with the overdamped dynamics of CSs in driven optical cavities pumped below transparency, where it was shown that the interaction of two nearby solitons causes a non-Newtonian motion where the velocity, rather than the acceleration, is proportional to the perturbation generated by one soliton onto the other [Vladimirov 2002].

Extreme event investigation

In this Section we will illustrate some preliminary results regarding the investigation of extreme events in the transverse plane of a broad-area VCSEL with optical injection.

The chaotic regime where we searched for extreme events corresponds to low values of optical injection. From the homogeneous stationary curve (blue solid line), depicted in Fig. 3.1, and the linear stability analysis results, illustrated in Figs. 3.2 and 3.3, it seems clear that the model in Eqs. (3.2) is unfit to this study, given the absence of a diffusive term in the equation for the electric field, which induces nonphysical short wavelength instabilities for intensity values I 0.23. Simulations of the model in Eqs. (3.1) require a much larger computational time, due to the equations stiffness, furthermore the large sampling required for the investigation of extreme events tends to slow even further the time required for each simulation, making the complete model unfit to this study as well. In order to investigate extreme events in this system we will then adopt the model described in Eqs. (3.3), where the diffusive term in the electric field equation sets a limit to the number of Hopf unstable modes. The chaotic nature of the regime observed has not yet been rigorously characterized in 2D+time, but we would like to highlight that in [Zamora-Munt 2013] for the purely temporal description of the same kind of laser, a chaotic regime was observed and characterized through Lyapunov exponent analysis. Furthermore the observed regime results visually similar (see for instance Fig. 3.11) to the extended spatiotemporal chaos observed in Chapter 1 and characterized in [Selmi 2016] in the 1D+time case through Lyapunov exponent analysis and the study of the Kaplan-Yorke dimension.

As illustrated in the Introduction of this Thesis, extreme events in a class-A laser with optical injection have been already observed numerically in [Gibson 2016], with the model of a forced Ginzburg-Landau equation. In the present case we are considering a class-B laser, where the carrier density dynamics is not fast enough for D to be adiabatically eliminated. Furthermore the laser considered is a VCSEL where the semiconductor nature of the system implies the addition of a linewidth enhancement factor in the equation for the electric field. The initial idea of the study consisted in searching the validity limits of the re- sults observed in [Gibson 2016] for vortex turbulence. Unfortunately inserting a linewidth enhancement factor α, with a reasonable value, appears to change the complexity of the turbulence observed.

We then studied the system dynamics in the chaotic regime occurring for low values of injection and usual values of linewidth enhancement factor, setting aside, for the time being, the search for the validity limits of the results of [Gibson 2016]. Similarly to what was done in Chapter 1 we focused on the statistics of the spatiotemporal maxima of the transverse field intensity, detected through the new method already illustrated in Section 1.5. The simulations here displayed were performed on a spatial grid of 256 × 256 pixels (corresponding roughly to 256 × 256 µm 2 ). As in the case of Chapter 1, the sampling of the transverse plane occurs every 1 ps and the integration time unit corresponds to 100 fs. The statistical analysis has been performed on a window of 25 ns after the system had been initialized for a time of 10 ns.

In Fig. 3.12 we have illustrated the probability density functions (PDFs) of the total intensity (a) and the spatiotemporal maxima (b) for the parameters σ = 400, α = 4, d = 0.01, µ = 6 and E I = 0.5. In particular, with respect to the parameters considered for the study of CS interaction we have increased the value of the pump µ and reduced the injection amplitude y so that the system results below the injection locking point, where it is Hopf unstable as already illustrated in Figs. 3.4. Even if the set of parameters is not optimized to obtain the maximum number of extreme events, the data shows a considerable deviation in (a) from the negative exponential tail, expected in case of Gaussian statistics, and a fair amount of events above the third extreme event threshold in (b). From these preliminary results we can conclude that extreme events are present in the system for this specific choice of parameters.

In Fig. 3.13 we illustrated an example of the kind of extreme events observed in this system. We can observe that a series of ripples of the electric field inten- sity is developed, starting from the extreme event spot: this appears even clearer in Figs. 3.14 where we plotted a color map of (a) the same electric field intensity, (b) the corresponding phase of the electric field and (c) the carrier density D at the instant where the extreme event (highlighted by the dark red cross) occurs.

In general, for this specific set of parameters, extreme events appear to occur in spots where the carrier density presents a higher value with respect to the background, similarly to what we have observed in Chapter 2. Given the high value of the relative decay rate σ, the evolution of the carrier density results very slow when compared to that of the electric field: this allows for an easier localization of the extreme event and its evolution in time, unlike the case of Chapter 1 where the structures observed in the intensity would move around in the transverse plane, making tracking them a difficult task.

In Fig. 3.15 we illustrated the trajectory (in time) of the system in the (D, φ, I) phase space in the spot where this specific extreme event is observed: as in the case of Chapter 2, the trajectory of the system towards a high-peak event is preceded by an evolution down to very low values of intensity and high values of carrier density, in the direction of the fixed point C (an unstable focus) on the lower branch of the homogeneous stationary curve. Thanks to the repulsive nature of the point C, the system is then expelled towards high values of intensity and finally ends up in a more bounded state for low values of carrier density. Even if some similarities can be drawn with the case of high-peak events observed in Chapter 2 in the unstable roll regime, one main difference consists in the role of the phase: in fact, as illustrated in Fig. 3.15, in the temporal evolution of the spot where the extreme event occurs, the phase continues to decrease and does not exhibit any counterclockwise rotation. Instead, as illustrated in the zoom in Fig. 3.16, around the extreme event plenty of vortices occur: these objects are clearly identified in (b) by the crossing of the zero value isoline for the real and the imaginary part of the electric field.

The results here shown consist just in a preliminary analysis of the system dynamics. A full study of the extreme event dependence on the parameters considered (especially the injection y and the pumping current µ) as well as a deeper investigation of the extreme events generating mechanism and its similarities with the vortex turbulence observed in [Gibson 2016], remain to be performed and will be addressed in a future work.

Conclusions

In this Chapter we have studied the model of a broad-area VCSEL with optical injection, described as a class-B laser. The main focus of the Chapter consisted in the study of the interaction between two cavity solitons in the transverse plane of the electric field intensity for values of pumping current above the free running laser threshold. We observed that CSs can be described as two particles subject to an interaction potential exponential decreasing with the distance between the two objects. This assumption lead us to an analytic expression for the cavity soliton merging time, which has, in good approximation, an exponential dependence on the initial cavity soliton distance. The numerical results confirmed this analytic description and gave us some specific values for the interaction range and strength of the potential: the attractive force is experienced by cavity solitons at distances much larger than the structures diameters. The exponential dependence of the cavity soliton merging time on the initial distance presents some clear differences with the similar study of cavity soliton interaction for the same system, not pumped and below transparency [Tissoni 1999], where some equilibrium distances and a particular maximal distance of interaction were observed. Furthermore this potential allows for a possible analogy with the attraction experienced by hydrophobic materials, that is nonpolar molecules and surfaces in water [Israelachvili 1982].

A preliminary analysis of the chaotic regime present in the system for low values of injection (below the injection locking point) lead to the observation of extreme events for a specific choice of parameters, in particular for high values of pumping current. Extreme events seem in this case to occur close to large values of carrier density, in spots where the phase of the electric field continuously decrease and generates some sort of ripples in the electric field intensity through the transverse plane. The dependence of extreme events on the optical injection amplitude and the values of pumping current considered will be addressed in a future work, as well as the investigation of the physical mechanism behind the formation of these objects and its connection with the vortex turbulence identified in [Gibson 2016] and the dynamics observed in [Rimoldi 2017a].

Conclusions

In this Thesis we presented the study of extreme and abnormal events in three different active optical systems. Localized structures represent an important possible solution in all these models, hence their interaction and the role played in the formation of extreme events have been also investigated into details. The main focus in all the work here presented consisted in searching for the physical and dynamical causes behind the formation of extreme events. Differently from most of the work developed in fibers and other optical systems, all the three types of laser considered in this Thesis present spatial degrees of freedom, either in the transverse plane (Chapters 1 and 3), perpendicular to the direction of propagation of light, or in the propagation direction (Chapter 2).

In Chapter 1 we studied the model of a monolithic broad-area semiconductor laser with an intracavity saturable absorber. A detailed analysis of all the possible solutions of the system was performed, with a particular attention to cavity solitons and extended spatiotemporal chaos. Extreme events were investigated in the latter and analyzed in their statistics and parameter dependence, showing a larger amount of extreme events in the case of fast absorbers and low pump values. This study required the development, for the first time, of an algorithm for the detection of the spatiotemporal maxima of the electric field intensity in three dimensions (2D+time). The comparison with other works in the community [Selmi 2016, Coulibaly 2017], developed in 1D+time, showed that the inclusion of 2D spatial effects may affect the way extreme events manifest themselves in the system. The striking similar shape, in space and time, between these events and cavity solitons, both stationary and self-pulsing suggests an enhancement of extreme events when the system results close to the dissipative soliton attractor [Rimoldi 2017b]. Furthermore, these similarities suggest that the same generating mechanism could lie at the formation of both structures: that is modulational instability of the homogeneous solution, in the spatial case, and Hopf instability of the homogeneous and cavity soliton solutions, leading to Q-switching, in the temporal case. The investigation of predictability in the system as well as the study of its conservative limit are among the future development of this work.

Chapter 2 is dedicated to the study of a semiconductor laser with optical injection, spatially extended along the propagation direction. After a careful derivation of the model and its comparison with other models, also used to describe this system, we analyzed into details three different regimes, each highlighting a different aspect of possible structure interactions. The main kind of localized structures observed in this system are phase solitons [START_REF][END_REF]], that is objects propagating inside the cavity maintaining their shape in time and carrying a 2π chiral charge. The first study focused on the attractive interaction between phase solitons, leading to phase solitons complexes [Gustave 2017], carrying multiple chiral charges. One interesting result consisted in the observation of an exponential dependence of the merging time between two separate structures from their initial distance [Gustave 2017]. The second and the third studies focused more on the search for extreme events. In particular in a first configuration, for high enough pump values, extreme events have been observed due to the collision of phase solitons with another kind of transient structure carrying a negative chiral charge [Walczak 2017]. In a second configuration instead, high-peak events have been observed in an unstable roll pattern regime where phase solitons are not a stable solution. In this case the interplay of positive and negative chiral charges can be suggested as a predictor for the formation of high-peak events in the system [Rimoldi 2017b]. The three different regimes studied numerically, both in the context of phase soliton merging and extreme events, present excellent agreement with the experimental results corroborated, in the last case, by a detailed statistical analysis. There are some clear similarities in the dynamical description of abnormal event appearance in the regime of unstable rolls and in the case of a system without space of [Bonatto 2011, Zamora-Munt 2013]. In particular, in both cases, these events seem associated to an evolution of the trajectory toward an unstable focus, presenting a high value of carrier density. As the system approaches the fixed point, the electric field intensity becomes almost zero, thus increasing the value of carrier density, then due to the fixed point repulsive nature the system trajectory is expelled toward lower values of carrier density, performing a large excursion in intensity. Differently from [Zamora-Munt 2013] the system appears to feel the repulsive effects of the unstable focus at a larger distance. The collisional regime highlights instead another mechanism for extreme event formation where space is a necessary requirement. In this system chirality plays a fundamental role. In particular the observed phase and amplitude dynamics appear reminiscent of the regime of defect mediated turbulence observed in [Coullet 1989] and vortex turbulence, which has been recently associated to rogue wave formation in oscillatory two-dimensional transverse spatial systems with coherent forcing [Gibson 2016].

Chapter 3 is focused on the study of cavity soliton interaction in a monolithic broad-area VCSEL with optical injection, pumped above threshold. Some preliminary results on the study of extreme events in this system are also given. Differently from the model described in [Gibson 2016], the present system is a class-B laser, where the dynamics of the carrier density is not fast enough to justify its adiabatic elimination. We observed that cavity solitons can be described as two particles subjected to an interaction potential exponentially decreasing with the distance between the two objects. As a consequence the merging time between two localized structures presents an exponential dependence on their initial distance, a result confirmed both analytically and numerically [Anbardan 2017]. Furthermore this observation allowed for a possible analogy with the attraction experienced by hydrophobic materials. In contrast with previous work [Tissoni 1999], developed for the same system below transparency, no critical distance, above which cavity solitons stop interacting, has been found. From the point of view of information processing, cavity solitons in semicon-ductor lasers have been long considered good candidates for optical information encoding, due to the compactness and fast response of the device and to the plasticity of cavity solitons themselves [Firth 1996, Brambilla 1997]. In this context cavity soliton interaction can be considered as a limiting aspect of optical information storage since, in order to treat cavity solitons as independent bits of information, ideally no interaction effect should be exhibited: hence the presence of a critical distance, defines the maximum information density of a possible optical memory array. Nevertheless we would like to point out that cavity soliton interaction could be positively exploited in order to develop logic gates [Columbo 2014]. Furthermore, from the observations of Chapter 2, it is clear that information could be encoded as well in the phase of localized structures (phase soliton complexes). Extreme events have been preliminary investigated in this system in a chaotic regime for low values of injection. In particular we observed extreme events for a specific parameter choice, with high values of the pumping current. Extreme events appear, in this case, to occur through a mechanism similar to [Rimoldi 2017a], in particular, once the trajectory of the system in a single spot of the transverse plane approaches low values of intensity (and high values of carrier density), getting closer to an unstable focus, it is repelled to low value of carrier density, through a large excursion in intensity. Differently from [Rimoldi 2017a], the phase of the electric field appears to monotonically decrease in time in the spot of the transverse plane where the extreme event is about to occur. Further investigation on the extreme event parameter dependence and their formation mechanisms, in connection with the vortex turbulence observed in [Gibson 2016] and the dynamics observed in [Rimoldi 2017a] are to be addressed in future work.

From the results here presented it is clear that, even if extreme events can be observed in system without space, describable through purely temporal models [Bonatto 2011, Zamora-Munt 2013], spatial effects may play an important role in the formation of extreme events [Bonazzola 2013], either enabling or disrupting some of the dynamical and physical mechanisms for their formation and preserving others, when the spatial dimensionality of the system is altered.

As stated before, the general study of extreme events in active optical systems allows to address the dynamical and physical mechanisms behind the formation of these objects, in systems that are mainly deterministic. Apart from the obvious positive effects of a deeper understanding for these disruptive phenomena, the study of extreme events will also have as a consequence, a general deeper knowledge of the systems considered. The connection of this field to the wider topic of extreme events in physics will allow for interesting developments in possible analogies with other physical systems, especially considering the high capability to observe, control and study these events in optics.

LSA equations in physical variables

The variables in Eqs. (1.1) are rescaled as follows [Bache 2005]:

D = η 1 N 1 N 0,1 -1 , D = η 2 N 2 N 0,2 -1 , F = 0 v g n 2 a 1 τ 1 2hω E .
The pump and absorption parameters are given as

µ = η 1 (1 -Bη 1 ) I 1 I 0,1 -1 , γ = η 2 (1 -Bη 2 ) 1 - I 2 I 0,
∂ F(k x , k y , t) F(k x , k y , t) = -k 2 (d + i)∂t log F(k x , k y , t) = -k 2 (d + i)t F(k x , k y , t) = e -k 2 (d+i)t .
Then at a following time t + ∆t F(k x , k y , t + ∆t) = e -k 2 (d+i)(t+∆t)

F(k x , k y , t + ∆t) = F(k x , k y , t)e -k 2 (d+i)∆t .

We can in this way appreciate the advantage of solving the spatial part of the equations evolution in the Fourier space: any spatial derivative leads to a simple multiplication in the Fourier space.

In our programs we implement the split-step method executing first a half step ∆t/2 in the Fourier space, then a full step ∆t of the purely temporal part of the equations through a Runge-Kutta method at second order and then completing with another half step in the Fourier space.

where we now have to normalize with all the possible changes in the order of the k -1 events. Setting q i = m i -1 we obtain, in the limit 1 k t:

P k (t) = 1 (k -1)! t-(k-1) ∑ q 1 =1 t-(k-2) ∑ q 2 =q 1 +1 . . . t-1 ∑ q k-1 =q k-2 +1
P k|1,q 1 ,...,q k-1 (t)

= 1 (k -1)! 1 t t-(k-1) ∑ q 1 =1 1 q 1 t-(k-2) ∑ q 2 =q 1 +1 1 q 2 . . . t-1 ∑ q k-1 =q k-2 +1 1 q k-1
Passing to the integral

P k (t) ∼ 1 (k -1)! 1 t t-(k-1) 1 dq 1 q 1 . . . t-1 q k-2 +1 dq k-1 q k-1 = 1 (k -1)! 1 t t-(k-1) 1 dq 1 q 1 . . . t-2 q k-3 dq k-2 q k-2 log t -1 q k-2 + 1 ∼ 1 (k -1)! 1 t t-(k-1) 1 dq 1 q 1 . . . t-2 q k-3 dq k-2 q k-2 log (t -1) = 1 (k -1)! log (t -1) t t-(k-1) 1 dq 1 q 1 . . . t-3 q k-4 dq k-3 q k-3 log t -2 q k-3 + 1 ∼ 1 (k -1)! log (t -1) log (t -2) t t-(k-1) 1 dq 1 q 1 . . . t-3 q k-4 dq k-3 q k-3 ∼ 1 (k -1)! log (t -1) log (t -2) • • • log (t -(k -1)) t ∼ 1 (k -1)! (log t) k-1 t .
Considering then a more generic time span between t w and t w + t the probability of the occurrence of k events in such a time span is the following Eq. (D.1) can be recognized as a Poisson distribution for the variable log (t w + t)/t w also called log-Poisson. We can also observe that the logarithmic rate of events approaches the constant value α for t t w . The total number of records produced by α independent sequences is then a sum of log-Poisson variables, therefore it is also log-Poisson distributed [Sibani 2003]. Finally from Eq. (D.1) we can observe that in the limit 1 k t the probability goes as Hence once solved Eq. (E.2) for x, we find also the solutions for n and n.

P k (
In the limit of large energies x 0 V(x) ≈ e x (E.4)

x 0 V(x) ≈ -x
hence the energy of the system E E = ω 2 LSA V(x max ) ≈ ω 2 LSA e x max ≈ -ω 2 LSA x min (E.5) where x max and x min are as indicated in Fig. 1.32. It is then possible to obtain an analytic expression for x, n and n for large energies in the conservative limit.

Setting as initial condition in x(t = 0) = x min , we obtain the following expressions for the first period of the particle, which is the time that takes to the particle to pass from x min to x max and then get back. where the first solution is related to times -t 0 + 2mt max < t < t 0 + 2mt max and the second solution is related to times t 0 + 2mt max < t < -t 0 + 2(m + 1)t max . We verified that this analytic approximation well describes the behavior of the conservative system in the limit of high energies. About the expression we found for t max , it is important to draw a quick comparison with what is shown in [START_REF] Oppo | Toda potential in laser equations[END_REF]] regarding the equation relating the period of oscillation to the energy in the system (there multiplied by the frequency, due to a renormalization): in fact the same expression holds in our system. Furthermore, in the limit of large energies H = e x max , if we define E + 1 = H due to a different definition of the Toda potential in the article considered, we obtain

ω LSA T = 2ω LSA t max = 2 √ 2e x max + 2 e x max arctanh √ 1 -e -x max = √ 2 2e
x max /2 + e -x max /2 log e x max /2 + (e x max -1) 1/2 e x max /2 -(e x max -1) 1/2 = √ 2 2(E + 1) 1/2 + (E + 1) -1/2 log (E + 1) 1/2 + E 1/2 (E + 1) 1/2 -E 1/2 .

Hence we conclude that the only actual difference between the period of oscillation in a class B laser and in a LSA consists in a different frequency of relaxation (ω for the class B laser or ω LSA for the LSA).

Given the analytic expressions for x, n and n in case of conservative large energies we can obtain an analytic approximation also for the variation of energy along a period. In particular the time derivative of the total energy of the system is, according to the complete model in Eq. ( 1 We can notice that this equation can be zero (since D < 0), differently from the case of a class B laser [START_REF] Oppo | Toda potential in laser equations[END_REF]], where it is always negative. Inserting the solution x(t) in Eq. (E.6) we obtain a complicated expression in the solution of the integrals, depending on x max and all the system parameters, which we are not going to report in this Appendix. The value of x max for which Eq. (E.6) equals zero should allow us to find the Q-switching solution for any set of parameters. The details of this study will be addressed in some future work.

Period and width of the pulses

In [START_REF] Oppo | Toda potential in laser equations[END_REF]] and [Cialdi 2013] two particular temporal intervals are considered to verify the validity of the class B laser description through the Toda potential: these are T 1 and T 2 (which correspond to 2t max and 2t 0 ). T 1 represents the laser period of oscillation, that is the time that takes to the particle to pass from x min to x max and get back. On the other hand T 2 is the time that takes to the particle to get from 0 to x max and come back. These two intervals have also a clear description on terms of usual laser dynamics: in fact T 1 corresponds to the temporal interval between two consecutive intensity peaks and T 2 is the pulse width, which corresponds to the temporal interval where the laser has an intensity higher than its stationary value, as it has been illustrated in Fig. E.1(a). We can apply the same kind of analysis to our system, in particular we have

T 1,2 = 2 x max
x min ,0 dx p = T 0 √ 2π

x max

x min ,0 dx E -V LSA (x) . (E.7)

with T 0 = 2π period of the small oscillations with time determined by the frequency ω LSA of the relaxation oscillations.

Using the semi-open extended Simpson formulas we can numerically evaluate the two integrals and hence make a comparison between the two temporal intervals below and above the the Hopf bifurcation.

For µ > µ h the system exhibits relaxation oscillations: in Fig. E.2(a) we have plotted in solid line the curves predicted from the theoretical model (in blue for T 1 /T 0 and red for T 2 /T 0 ) and compared them with the data obtained from the measured periods during the full system evolution. For µ < µ h the system, after a transient, reaches a pulsed solution of fixed amplitude, typical of the Q-switching regime. In Fig. E.2(b) we can observe that also the transient phase is well described by the Toda potential and in particular the periods of the final pulsed solution are well aligned with the theoretical curve. Finally we did different simulation where we kept fixed the parameters s = 3, b = 10 -4 , r = 1 and the initial conditions (I 0 , n 0 , n0 ) = (0.001, 0, 0) and meanwhile changed the pump parameter µ, always below the Hopf bifurcation. In Fig. E.2(c) we reported the data of the final pulsed solution for 30 different values of µ thr < µ < µ h : we can notice a good agreement with the theoretical model. We hence conclude that the Toda oscillator model well describes both the nonlinear damped oscillation above the Hopf bifurcation and the passive Q-switching regime below the bifurcation.

In the large energy limit we can obtain an analytic expression for the two in- we can observe that the curves given by the analytic approximation of these two temporal intervals seem almost the same as the theoretical curves given by the numerical integration of Eq. (E.7).

Appendix F

Split-step method for Eqs. (2.7), (2.8) Let us consider the rate-equation model with a term of diffusion, same as in Eqs.

(2.8),

∂E ∂η

+ ∂E ∂τ -d ∂ 2 E ∂η 2 = T [y -(1 + iθ) E + (1 -iα)DE] , ∂D ∂τ = bT σ µ -D 1 + |E| 2 .
As in Appendix B we separate the purely temporal part of the evolution of the dynamical equations from the spatial part, which will be computed in the Fourier space.

In this case the temporal part is the following

∂E ∂τ = T [y -(1 + iθ) E + (1 -iα)DE] , ∂D ∂τ = bT σ µ -D 1 + |E| 2 .
which is integrated through a Runge-Kutta method at the second order, and the spatial part is contained just in the differential equation for E ∂E ∂η

+ ∂E ∂τ -d ∂ 2 E ∂η 2 = 0 , (F.1)
The Fourier transform of E is defined as follows

F [E(η, τ)] = Ê(k η , τ) = +∞ -∞
E(η, τ)e -ik η η dη with F Fourier transform operator, and its inverse transform is

E(η, τ) = +∞ -∞
Ê(k η , τ)e ik η η dk η .

Applying the operator F to Eq. (F.1), we obtain Then at a following time τ + ∆τ Ê(k η , τ + ∆τ) = e -k η (i+kηd)(τ+∆τ)

Ê(k η , τ + ∆τ) = Ê(k η , τ)e -k η (i+kηd)∆τ .

Similarly, for the pure rate-equation model in Eq. (2.7), we have Ê(k η , τ + ∆τ) = Ê(k η , τ)e -ik η ∆τ .

In our programs we implement the split-step method executing first a half step ∆t/2 in the Fourier space, then a full step ∆t of the purely temporal part of the equations through a Runge-Kutta method at second order and then completing with another half step in the Fourier space. Then at a following time t + ∆t Ê(k x , k y , t + ∆t) = e -iσk 2 (t+∆t) , D(k x , k y , t + ∆t) = e -dk 2 (t+∆t) ,

  ] . . . . . . . . . . . . . . . . . 1.1 Example of a VCSEL with intracavity SA . . . . . . . . . . . . . . 1.2 Scheme of the LSA system . . . . . . . . . . . . . . . . . . . . . . . 1.3 HSSs and ω s for LSA . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Stationary and Hopf instability domains . . . . . . . . . . . . . . . 1.5 HSSs and other solutions for the LSA . . . . . . . . . . . . . . . . . 1.6 Stationary CS in the intensity transverse plane and its Fourier optical spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Extended spatiotemporal chaos in the intensity transverse plane and its phase . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.8 Fourier optical spectrum for the extended spatiotemporal chaos both in linear and logarithmic scale . . . . . . . . . . . . . . . . . . 1.9 Fourier spectrum profiles for CS, CCS and STC . . . . . . . . . . . 1.10 CCS in the intensity transverse plane and its Fourier optical spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.11 Phase of the electric field for the CCS . . . . . . . . . . . . . . . . . 1.12 OCS in the intensity transverse plane and its Fourier optical spectrum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.13 Fourier spectrum profile and time trace for a OCS and plot of the pulsation period as function of µ . . . . . . . . . . . . . . . . . . . 1.14 Fourier spectrum for different values of diffusion coefficient . . . 1.15 Spatiotemporal maxima selection . . . . . . . . . . . . . . . . . . . 1.16 3D spatiotemporal grid representation . . . . . . . . . . . . . . . . 1.17 PDFs of the spatiotemporal maxima . . . . . . . . . . . . . . . . . 1.18 PDFs of the total intensity . . . . . . . . . . . . . . . . . . . . . . . 1.19 Density plots of RW ratio and kurtosis for the PDFs of spatiotemporal maxima . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.20 Density plots of RW ratio and kurtosis for the PDFs of total intensity 1.21 Time trace of the spatially averaged intensity for µ = 5 and r = 2.4 1.22 Kurtosis and RW ratio for fixed r and µ . . . . . . . . . . . . . . . 1.23 HSSs and extended spatiotemporal chaos branch for different parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.24 Kurtosis and RW ratio for fixed r for different parameters . . . . . iii 1.25 Results from [Selmi 2016] . . . . . . . . . . . . . . . . . . . . . . . 1.26 Time trace of the intensity averaged on a small area . . . . . . . . 1.27 PDFs of the spatiotemporal maxima and of the temporal maxima of a time trace spatially averaged on a small area . . . . . . . . . . 1.28 PDF of the logarithmic waiting times between two consecutive extreme events . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.29 Extreme event in the intensity transverse plane and PDFs of spatiotemporal maxima and total intensity . . . . . . . . . . . . . . . 1.30 Temporal and spatial profiles for I, D and D of the extreme event in Fig. 1.29 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.31 Comparison between the spatial and temporal profiles of cavity solitons and extreme events . . . . . . . . . . . . . . . . . . . . . . 1.32 Toda potential . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

  v

event120 2 .

 2 48 Numerical phase time trace . . . . . . . . . . . . . . . . . . . . . . 2.49 Spatiotemporal diagram and phase time trace with highlighted 5 abnormal events . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.50 Phase portrait in the (Re(E), Im(E)) and (D, I) plane around an abnormal event . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.51 Spatiotemporal diagram for intensity and carrier density around an abnormal event . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 HSS for the semiconductor laser with optical injection. . . . . . . 3.2 Stationary instability domain of the model in Eqs. (3.2) for a specific choice of parameters. . . . . . . . . . . . . . . . . . . . . . . . 3.3 Hopf instability domains of the model in Eqs. (3.2) for specific choices of parameters. . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Stationary and Hopf instability domains of the model in Eqs. (3.3) for a specific choices of parameters. . . . . . . . . . . . . . . . . . . 3.5 CSs merging process for σ = 400 and d 0 = 10 spatial units. . . . . 3.6 Time evolution of cavity soliton distances for different initial values and σ = 400. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7 Merging time of two CSs as a function of their initial distance. . . 3.8 Temporal evolution of the distance between two CSs for σ = 400 and r 0 = 10. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 CSs merging time as a function of their initial distance for different values of σ. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.10 CSs merging time as a function of 1/σ for different initial distances.139 3.11 Example of the chaotic regime observed for low values of optical injection. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.12 PDFs of total intensity and spatiotemporal maxima. . . . . . . . . 3.13 Example of extreme event in the electric field intensity. . . . . . . 3.14 Color maps of the electric field intensity, the phase of the electric field and the carrier density in presence of an extreme event. . . . 143 3.15 System trajectory in the (D, φ, I) phase space in the spot of the extreme event. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 144 3.16 Phase and zero isolines of the real and imaginary part of the electric field around the extreme event. . . . . . . . . . . . . . . . . . . 144 E.1 Relaxation oscillations in a class B laser and comparison between the numerically integrated and theoretical periods T 1 and T 2 . . . 164 E.2 Comparison between the predicted periods and the data measured during the simulation . . . . . . . . . . . . . . . . . . . . . . 165

  Fig. 1 -Examples of extreme events in nature and society. (a) Hurricane Hugo occurring in 1989. Photo by NOAA / Satellite and Information Service, from Wikimedia Commons -Public Domain. (b) Crowd gathering on Wall Street after the 1929 Wall Street Crash which many believe lead to the following 12-year Great Depression affecting all Western industrialized countries. Photo from Wikimedia Commons -Public Domain.

Fig. 4 -

 4 Fig. 4 -Reprinted with permission from [Liu 2015]. Rogue wave formation in the spatial distribution of the electromagnetic energy. The color bar values are rescaled on the significant wave height of the field intensity. The white circle highlights the ultrafast, subwavelength rogue wave triggered by spontaneous synchronization.

Fig. 5 -

 5 Fig.5-Reprinted with permission from[Zamora-Munt 2013]. System dynamics in the phase space when approaching and extreme event. S1, S2, S3 are three fixed points of the system, in particular S1 is a saddle, S2 and S3 are unstable foci. In (a) the blue dotted and green solid lines indicate, respectively, the stable and unstable manifolds of S2. (b) illustrates the system trajectory (highlighted in red) when approaching an extreme event, zoomed in (c). (d) Shows the intensity time trace.

Fig. 6 -

 6 Fig. 6 -Reprinted with permission from [Gibson 2016]. Vortices and extreme event formation in the transverse plane of the electric field intensity. (a) hexagonal Turing pattern whose breaking leads to (b) optical vortex-mediated turbulence. (c) is a zoom on the extreme event highlighted in (b) where we can distinguish the presence of two optical vortices where the phase (d) in the transverse plane performs full 2π rotations in the opposite directions.

Fig. 1

 1 Fig. 1.1 -Example of experimental setup used at the laboratory of C2N in Marcoussis, Paris. Courtesy of C2N -Paris

2 . 2 -

 22 Fig. 1.2 -Scheme of the system studied, a broad-area monolithic VCSEL with an intracavity saturable absorber. Light is emitted perpendicularly to the mirrors in the Laser. The transverse profile presents for certain parameter conditions, chaotic patterns, as the one here represented.

µ

  = (1 -Ds,nt ) [B(1 -Ds,nt ) + 1 + I s ] (1.6) which is plotted in Fig. 1.3(a) together with the trivial solution for a specific set of parameters. For completeness we also plot in Fig. 1.3(b) the laser frequency ω s as function of I for the same choice of parameters.

( a )

 a Trivial and nontrivial HSSs. (b) ω s as in Eq. (1.5).

Fig. 1

 1 Fig. 1.3 -(a) Trivial I 0 = |F s,0 | 2 and nontrivial I s HSSs for the electric field intensity in Eq. (1.6) as functions of the pump parameter µ. (b) Laser frequency ω s as in Eq. (1.5) as function of I. The parameters used in both plots are α = 2, β = 1, γ = 2, s = 1 and B = 0.1. µ th = 5.18 is the laser threshold. Dashed lines highlight the unstable branches, as illustrated in the following stability analysis.

Fig. 1

 1 Fig. 1.4 -Pattern instability domains plotted for the square modulus of the wave vector k 2 as a function of the laser intensity I. The common parameters considered are: s = 1, B = 0.1, α = 2, β = 1, γ = 2 and d = 0.01. In (b) a specific value of r = 1.75 was considered. The grey-shaded area is (a) stationary or (b) Hopf unstable while the red dashed (and solid) line corresponds to the instability boundary for d = 0.

  Fig. 1.5 -(a) Homogeneous stationary solution (in black), as given in Eqs. (1.2) and (1.3) plotted together with the stationary soliton branch (blue line with triangle markers) and the time-averaged maximum intensity of the chaotic solitons (green line with diamond markers) and of the extended spatiotemporal chaos (orange line with square markers) for r = 1. (b) Homogeneous stationary solution (in black) plotted together with the oscillatory cavity soliton branch for r = 1.75. All the data here plotted are relative to simulations for b = 0.01, s = 1, B = 0.1, α = 2, β = 1, γ = 2 and variable µ. µ th = 5.18 is the laser threshold.

  (a) Stationary CS in the transverse plane of the electric field intensity. (b) Transverse optical spectrum of the stationary CS (in logarithmic scale).

Fig. 1

 1 Fig. 1.6 -Stationary cavity soliton (CS) simulated for the following parameters b = 0.01, r = 1, s = 1, B = 0.1, α = 2, β = 1, γ = 2, δ = 0.01, µ = 4.8 and represented (a) in the transverse plane of electric field intensity and (b) in the Fourier space.

( a )

 a Extended spatiotemporal chaos in the transverse plane of the electric field intensity. (b) Phase of the electric field in the extended spatiotemporal chaotic regime.

Fig. 1

 1 Fig. 1.7 -Intensity (a) and phase (b) of the electric field in the transverse plane in the regime of extended spatiotemporal chaos simulated for the following parameters b = 0.01, r = 1, s = 1, B = 0.1, α = 2, β = 1, γ = 2, δ = 0.01, µ = 4.8.

Fig. 1

 1 Fig. 1.8 -Transverse optical spectrum in linear (a) and logarithmic (b) scale, in the regime of extended spatiotemporal chaos, simulated for the following parameters b = 0.01, r = 1, s = 1, B = 0.1, α = 2, β = 1, γ = 2, δ = 0.01, µ = 4.8 (same as in Fig. 1.7).

Fig. 1

 1 Fig.1.9 -Averaged spectrum profiles of the CS in Fig.1.6 (in blue), the chaotic cavity soliton (CCS) in Fig.1.10 (in red) and the extended spatiotemporal chaotic regime (STC) in Fig.1.7 and 1.8 (in yellow). In (a) no normalization has been applied to the data, instead in (b) and (c) all the spectrum profiles are normalized to their maximum value. The black dashed line represent a possible spectrum profile with a power law decay that best fits the spectrum profile of the extended chaotic regime. We highlight that all these spectra have been acquired for the same value of pumping current µ = 4.8.

  (a) CCS in the transverse plane of the electric field intensity. (b) CCS in the transverse plane of the electric field intensity. (c) Transverse optical spectrum of the CCS. (d) Transverse optical spectrum of the CCS.

Fig. 1 .

 1 Fig. 1.10 -Chaotic cavity soliton (CCS) simulated for the following parameters b = 0.01, r = 1, s = 1, B = 0.1, α = 2, β = 1, γ = 2, δ = 0.01, µ = 5 and represented (a-b) in the transverse plane of electric field intensity and (c-d) in the Fourier space at two different times. The white circles in (a-b) correspond to the same spots highlighted in Fig. 1.11 and represent optical vortices.

( a )

 a Phase of the electric field for the CCS. (b) Phase of the electric field for the CCS.

Fig. 1 .

 1 Fig. 1.11 -Phase of the electric field for the CCS, simulated for the following parameters b = 0.01, r = 1, s = 1, B = 0.1, α = 2, β = 1, γ = 2, δ = 0.01, µ = 5, the white circles denote the presence of two optical vortices in the transverse plane. (a) has been recorded at the same time as Fig. 1.10 (a) and (c), (b) has been recorded at the same time as Fig. 1.10 (b) and (d).

Fig. 1 .

 1 10 we present (a-b) a CCS in the transverse plane of the electric field intensity at two different times and (c-d) the corresponding Fourier spectrum. The structure presented in Fig.1.10 continues to evolve and change in time but it remains localized in space on a specific area of the transverse plane. The transverse optical spectrum reflects this evolution of the intensity and, in time, its possible to observe that the modes interested by the CCS dynamics either change (aperiodically) in a rotation around the central (k x , k y ) = (0, 0) mode, as visible in Fig.1.10(c), or the energy associated is split in two lobes on the opposite sides of the central mode, like in the case ofFig. 1.10(d).

  (a) OCS in the transverse plane of the electric field intensity. (b) OCS in the transverse plane of the electric field intensity. (c) Transverse optical spectrum of the OCS. (d) Transverse optical spectrum of the OCS.

Fig. 1 .

 1 Fig. 1.12 -Oscillatory cavity soliton (OCS) simulated for the following parameters b = 0.01, r = 1.75, s = 1, B = 0.1, α = 2, β = 1, γ = 2, δ = 0.01, µ = 5 and represented (a-b) in the transverse plane of electric field intensity and (c-d) in the Fourier space at two different times.

  Fig. 1.12 we illustrated the behavior of the intensity in the transverse plane in presence of an OCS at two different times (a-b) as well as the corresponding Fourier transverse optical spectrum (c-d). In Fig. 1.13(a) we compared the two (a) Fourier spectrum profiles at peak and baseline points of the time trace. (b) Time trace of the OCS peak. (c) Period of OCSs as function of µ.

Fig. 1 .

 1 Fig. 1.13 -(a) Comparison of the spectrum profiles at the same times depicted in Fig. 1.12(c) and (d), the data plotted in blue refers to a case of very small peak intensity, and the data plotted in red refers to a high peak intensity. (b) Time trace of intensity in the OCS peak during the same simulation as in Fig. 1.12. (c) Period of the OCS solutions illustrated in Fig. 1.5(b) as a function of µ.

Fig. 1 .

 1 Fig. 1.14 -Fourier spectrum profile of the chaotic regime for different values of the diffusion coefficient d. The violet profile (d = 0.01) corresponds to the actual value chosen in the simulations.

  (a) t = 0 ps. (b) t = 9 ps. (c) t = 18 ps.

Fig. 1 .

 1 Fig. 1.15 -Snapshots of a zoom on the intensity transverse plane, separated 9 ps from each other to exemplify the spatiotemporal maxima method of selection. The peak highlighted by the arrows is not selected as long as it is increasing (a) or decreasing (c) in time, but only at its maximum value (b). The red circles indicate other spatiotemporal maxima occurring at the same time of the snapshots.

Fig. 1 .

 1 Fig. 1.16 -3D spatiotemporal grid representation illustrating the numerical procedure of selection for the maxima. Each point of the grid (central point) is compared, in intensity values, with all its 26 neighbors.

  The first threshold, mostly used in spatially extended systems[Arecchi 2011, (a) PDF of the spatiotemporal maxima for µ = 7 and r = 2.5. (b) PDF of the spatiotemporal maxima for µ = 5 and r = 2.4. (c) PDF of the spatiotemporal maxima for µ = 5 and r = 1.

  (a) PDF of the total intensity for µ = 7 and r = 2.5. (b) PDF of the total intensity for µ = 5 and r = 2.4. (c) PDF of the total intensity for µ = 5 and r = 1.

Fig. 1

 1 Fig. 1.18 -PDFs of the total intensity. (a) corresponds to a simulation for µ = 7 and r = 2.5, (b) corresponds to a simulation for µ = 5 and r = 2.4 and (c) corresponds to a simulation for µ = 5 and r = 1. All simulations lasted 25 ns. The vertical black dashed line corresponds to the first threshold described in the text, computed on the average and standard deviation of the total intensity. The light blue line is a negative exponential PDF computed on the average of the data as exp(-I/ I )/ I .

  (a) RW ratio according to threshold 1. (b) RW ratio according to threshold 2. (c) RW ratio according to threshold 3. (d) Kurtosis.

Fig. 1 .

 1 Fig. 1.19 -RW ratios of the PDF of spatiotemporal maxima according to the first (a), second (b) and third (c) RW threshold definitions introduced in the text. (d) Kurtosis of the different spatiotemporal maxima PDFs.

  (c) and (d) illustrate a slightly different situation, in comparison with (a) and (b), showing as the most favorable regime for extreme events an area for µ close to 5 and r close to 2.4, corresponding to the PDFs plotted in Figs. 1.17(b) and 1.18(b). Such favorable situation seems to disappear for values of r higher than 2.7, at difference with what occurs in (a) and (b). In its peak the kurtosis reaches a value of K ≈ 90 and the maximum ratio of extreme events associated to the same simulation, according to the third threshold, is RW ratio ≈ 0.003. (a) RW ratio according to threshold 1. (b) Kurtosis -K exp .

Fig. 1 .

 1 Fig. 1.20 -RW ratio of the PDF of total intensity according to the first RW threshold definition (a) introduced in the text. (b) Kurtosis of the total intensity PDFs, relative to the negative exponential (see text).

  Fig. 1.20(a) and (b) present more similarities with Fig. 1.19(a) and (b), especially showing a maximum value of RW ratio and kurtosis for higher values of r 2.5, instead than for r ≈ 2.4 and µ ≈ 5 as Fig. 1.19(c) and (d).

Fig. 1

 1 Fig. 1.21 -Time trace of the spatially averaged intensity for a simulation with µ = 5 and r = 2.4.

  (a) Kurtosis and RW ratio for r = 2.2. (b) Kurtosis and RW ratio for µ = 5.

Fig. 1 .

 1 Fig. 1.22 -Values of kurtosis, in green solid line with circle markers, and (third threshold) RW ratio, in blue dashed line with diamond markers, for (a) r = 2.2 and variable µ and for (b) µ = 5 and variable r.

Fig

  Fig.1.23 -Homogeneous stationary solution (in black), as given in Eq. (1.6) plotted together with the extended spatiotemporal chaos (orange line with square markers) branch. All the data here plotted is relative to simulations for b = 0.005, r = 1, s = 10, B = 0, α = 2, β = 0, γ = 0.5[Bache 2005] and variable µ. µ th = 1.5 is the laser threshold.

Fig. 1

 1 Fig. 1.24 -Values of kurtosis, in green solid line with circle markers, and (third threshold) RW ratio, in blue dashed line with diamond markers, for (a) r = 1 and variable µ. The simulations were run for fixed parameters b = 0.005, r = 1, s = 10, B = 0, α = 2, β = 0, γ = 0.5 [Selmi 2016].

Fig. 1

 1 Fig. 1.25 -Reprinted with permission from [Selmi 2016]. On the left side proportion p EE of extreme events (blue circles) and kurtosis γ 2 (red squares) as functions of µ. On the right side the Kaplan-Yorke D KY dimension as function of µ.

Fig. 1

 1 Fig. 1.26 -Time trace of the intensity, averaged on a central square of 5 × 5 µm 2 . The simulation was run for b = 0.005, r = 1, s = 10, B = 0, α = 2, β = 0, γ = 0.5 and µ = 1.49 (left boundary of the chaotic branch in Fig. 1.23).

  (a) PDF of the spatiotemporal maxima µ = 1.49.(b) PDF of the spatiotemporal maxima µ = 3.(c) PDF of the temporal maxima on a small area µ = 1.49.(d) PDF of the temporal maxima on a small area µ = 3.

Fig. 1 .

 1 Fig. 1.27 -(a) PDF of the spatiotemporal maxima, detected during the simulation with our usual method. (b) PDF of the temporal maxima, individuated on the time trace of the intensity, spatially averaged on a small 5 × 5 µm 2 area as in the experiments [Selmi 2016, Coulibaly 2017]. In (d) the Weibull PDF is not depicted because it does not describe well the data.

  (a) Extreme event in the transverse plane. (b) Spatiotemporal maxima PDF.

Fig. 1 .

 1 Fig. 1.29 -(a) Example of extreme event in the transverse plane for a simulation at µ = 4.8 and r = 2.2. (b) Spatiotemporal maxima PDF and (c) total intensity PDF for this specific simulation.

Fig. 1 .

 1 Fig. 1.30 -Temporal (a) and spatial (b) profile of the extreme event depicted in intensity I in Fig. 1.29, together with the corresponding profiles of the carrier populations D and D.

( a )

 a Comparison with the spatial profile of CS. (b) Comparison with the temporal profile of OCS.(c) Comparison with the spatial profile of OCS.

Fig. 1 .

 1 Fig. 1.31 -(a) Comparison of an extreme event spatial profile with that of a CS, both simulated for µ = 5 and r = 1. Comparison of an extreme event spatial (b) and (c) temporal profiles with that of an OCS, both simulated for µ = 5 and r = 1.75.

Fig. 1 .

 1 Fig. 1.32 -Toda potential V(x) = e xx.

  (1.24) (in the non conservative limit) become ẋ = 2 (D s n + Ds n) , (1.28a) ṅ = -b I s (e x -1) + n(1 + I s e x + 2BD s ) + n 2 BD s , (1.28b) ṅ = -rb sI s (e x -1) + n(1 + sI s e x + 2B Ds ) + n2 B Ds , (1.28c)

  Fig. 2.1 -Schematization of a semiconductor ring laser with optical injection.

Fig. 2 . 2 -

 22 Fig. 2.2 -Schematization of the semiconductor ring laser used for the model derivation.

  for the polarization at order 0 reads P = (1iα)DE .

Fig. 2 .

 2 Fig. 2.3 -(a) Stationary intensity I = ρ 2 s of the electric field as a function of the injection amplitude y, as indicated in Eq. (2.15) for θ = -3.04, and µ = 1.1. (b) Bistability boundaries as in Eq. (2.17) for the same parameters.

  Fig. 2.4 -Comparison between the instability domains of the complete model of Eqs. (2.5) (in green), the rate-equation model of Eqs. (2.7) (in blue) and the rate-equation model with diffusion of Eqs. (2.8) (in red). The parameters considered are α = 3, σ = 3 × 10 -6 , T = 0.3, µ = 1.1, θ = -3.04, b = 0.0001 and d = 10 -6 .

Fig. 2 .

 2 Fig. 2.5 -Stationary curve as in Eq. (2.15) for the parameters µ = 1.1, θ = -3.04 and α = 3. The highlighted fixed value of injection corresponds to y = 0.004 and will be relevant in the next Sections.

( a )

 a Stationary curve for µ = 1.1, θ = -2.7. (b) Stationary curve for µ = 2, θ = -2.7.(c) Instability domains for µ = 1.1, θ = -2.7, b = 0.0005. (d) Instability domains for µ = 2, θ = -2.7, b = 0.0001.

Fig. 2 .

 2 Fig. 2.6 -Stationary curves (a-b) and instability domains (c-d) for µ = 1.1, b = 0.0005 (a,c) and for µ = 2, b = 0.0001 (b,d). In (c-d) the instability domains are depicted for the rate-equation model with diffusion (in red) and for the pure rate-equation model (in blue). The other parameters are θ = -2.7, α = 3, σ = 3 × 10 -6 , T = 0.3.

  Zoom of the instability domain in (a).

Fig. 2 .

 2 Fig. 2.7 -Instability domains for the pure rate-equation model (in blue), the complete model (in green) and the rate-equation model with diffusion (in red) for a smaller value of d = 10 -12 . The other parameters are the same as in Fig. 2.4.

  Fig. 2.8 -Stationary curves for (a) µ = 1.01, θ = -3.01, (b) µ = 1.01, θ = -3.1, (c) µ = 1.1, θ = -3.01 and (d) µ = 1.1, θ = -3.1.

  Fig. 2.9 -Instability domains for (a) µ = 1.01, θ = -3.01, (b) µ = 1.01, θ = -3.1, (c) µ = 1.1, θ = -3.01 and (d) µ = 1.1, θ = -3.1. The instability domain for the rateequation model with diffusion is depicted in red and for the pure rate-equation model is depicted in blue.

Fig. 2 .

 2 Fig. 2.10 -Schematization of the main experimental setup. M are the mirrors, BS are the beam splitters, OI is an optical isolator and piezo is a piezoelectric actuator. SOA stays for semiconductor optical amplifier.

Fig. 2 .

 2 Fig. 2.11 -Schematization of a second experimental setup in a Fabry-Perot configuration. M are the mirrors, BS are the beam splitters, OI is an optical isolator.

( a )

 a Fig.2.12 -Spatiotemporal diagrams in the (a) experimental reference frame, where η is fixed at η = 1 (z = L) and (b) the numerical reference frame, where all the full spatial intensity profiles, acquired one each roundtrip time, are stacked onto each other.

  Fig. 2.13 -Homogeneous stationary solution of the model in Eq. (2.7) for µ = 1.1, α = 3 and θ = -2.7 with superimposed the stability domain of PSs. The square and triangle markers correspond respectively to the maxima and minima of the PS profile.

( a )

 a Fig. 2.14 -PSs at the boundaries of the stability domain: y = 0.0445 (left column) and y = 0.0469 (right column). (a,b) are the spatiotemporal diagrams, (c,d) represent the phase portraits and (e,f) are the temporal profiles of the PSs. The points 1 and 2 define the width of the PS. For the points A, B and C see the text.

  Fig. 2.15 -Trajectory of the system during a single roundtrip in the presence of a PS for a value of injection amplitude y = 0.0458. (a) is the trajectory in the phase space while (b) is the projection of the same trajectory in the (Re(E), D) plane.

Fig. 2 .

 2 Fig. 2.16 -Reprinted with permission from [Gustave 2016b]. On the left the authors depicted the spatiotemporal diagram in a regime where phase solitons are stable, on the top right the authors illustrated the phase portrait and, below it, the intensity temporal profile of the PS.

  (a) Velocity and width of PSs as functions of y (b) Velocity as function of the width of PSs

Fig. 2 .

 2 Fig. 2.17 -(a) Velocity and width of PSs for varying values of the injection amplitude y. (b) Dependence of the velocity on the width of the PS. The PS size is measured in percentage of the roundtrip time.

Fig. 2 .

 2 Fig. 2.18 -Spatiotemporal diagram of attractive PS interaction between two PSs carrying a single charge each, observed experimentally. This result has been obtained by P. Walczak in our team, using the Fabry-Perot experimental setup illustrated in Fig. 2.11.

Fig. 2 .

 2 Fig.2.19 -Attractive interaction between PSs of single (PS 1 ) and n -1 charges (PS n-1 ), as observed numerically. All the simulation results here depicted have been run for a value of injection y = 0.0458. The single charge PS is generated at a similar distance of approximately 20% of the roundtrip time from the other structure for all the simulations.

  Fig.2.20 -Merging of two PSs in two peculiar cases: in (a) the interaction between the two PS does not lead to a PS whose charge is the sum of the charges of the previous objects; in (b) the interaction causes the system to lose its stability. The simulations have been run for y = 0.0458 and the single-charge PS is generated at a distance of approximately 20% of the roundtrip time from the other structure.

( a )

 a Fig. 2.21 -Interaction between two single-charge phase solitons at the initial distances of 5%, 7% and 10%.

  Fig. 2.22 -Merging time between two single-charge PSs as a function of the initial distance dist ini between them, plotted in linear (a) and semi-logarithmic (b) scale.

Fig. 2 .

 2 Fig. 2.24 -Experimental PSs of charge 6 (a) and 11 (b), observed with the ring laser setup of Fig. 2.10.

Fig. 2 .

 2 Fig. 2.25 -Velocity as a function of the chiral charge of the PSs.

2 Fig. 2 .

 22 Fig. 2.26 -Experimental PS of charge 2, represented in its (a) spatiotemporal diagram, (b) trajectory in the Argand plane, (c) intensity temporal profile. The experimental setup used for these observations is the ring laser illustrated in Fig. 2.10.

Fig. 2 .

 2 Fig. 2.27 -Spatiotemporal diagrams (a,b) and phase portraits (c,d) of a phase soliton of charge 2 (a,c) and a phase soliton of charge 3 (b,d), observed experimentally with the setup illustrated in Fig. 2.11.

( a )

 a Phase space for a PS 2 (b) Same as (a), in the (Re(E), D) plane. (c) Same as (a), in the (Re(E), Im(E)) plane.

  (d) Phase space for a PS 3 (e) Same as (a), in the (Re(E), D) plane.

  Fig. 2.29 -Collisional regime as observe experimentally. (a) represents the spatiotemporal diagram where the black dashed lines highlight the extreme event with the maximum intensity observed. (b) is a zoom on such extreme event where the two arrows point out the collision between a PS of charge 2 and another structure propagating in the opposite direction. (c) is the intensity temporal profile of the extreme event of maximum intensity.

Fig. 2 .

 2 Fig. 2.30 -Experimental PDF of the intensity. The two vertical lines are two thresholds for extreme events computed as the average of the total intensity plus eight times the standard deviation computed on the full intensity time trace (black dashed line) or just on the time trace in the collisional regime (yellow dash-dotted line) of Fig. 2.29(a). For the different PDFs description see the text.

Fig. 2 .

 2 Fig.2.31 -Example of the kind of collisions observed numerically. We can distinguish PS complexes and other transient structures propagating, in this reference frame, in the opposite direction. The collisions between these two objects can give rise to multiple peaks as well illustrated here.

  (a) Spatiotemporal diagram. (b) Phase portraits.

Fig. 2 .

 2 Fig. 2.32 -(a) Zoom of the spatiotemporal diagram around a collision observed in the simulation. (b) phase portraits of the roundtrip sections highlighted in (a). On roundtrip section 3 a high intensity event occurs, preceded by a clockwise phase rotation on roundtrip section 2.

Fig. 2 .

 2 Fig. 2.33 -Example of a roll pattern in the (a) spatiotemporal diagram, (b) its trajectory in the phase space and (c) in the Argand plane. The simulation has been run for a value of injection Y=0.1.

  Fig. 2.34 -Stationary homogeneous solution with depicted (green triangles) the maximum and minimum values of the roll intensity profile for different values of injection Y. The stationary curve is depicted with a dashed (solid) line where it is unstable (stable). The blue dotted line marks the boundary of the roll solution. The orange markers correspond to the maxima and minima registered during the unstable roll regime.

  (a) Spatiotemporal diagram (b) (Re(E), Im(E)) plane (c) (I, D) plane (d) Spatiotemporal diagram (e) (Re(E), Im(E)) plane (f) (I, D) plane (g) Spatiotemporal diagram (h) (Re(E), Im(E)) plane (i) (I, D) plane (j) Spatiotemporal diagram (k) (Re(E), Im(E)) plane (l) (I, D) plane

Fig. 2 .

 2 Fig. 2.35 -Roll pattern breaking for decreasing values of injection in the spatiotemporal diagram (first column), in the Argand plane (second column) and in the (D, I) plane (third column). The simulations have been run for a value of injection Y = 0.02 (a-c), 0.006 (d-f), 0.004 (g-i) and 0.003 (j-l).

Fig. 2 .

 2 Fig. 2.36 -Spatiotemporal diagram observed experimentally, where abnormally high events emerge above an unstable roll background.

  Fig. 2.37 -Stationary homogeneous solution of the model in Eqs. (2.8) for the intensity as a function of the detuning θ. The parameters considered are α = 3, µ = 1.1 and Y = 0.004. The vertical red dash-dotted line indicates the chosen value of detuning θ = -3.04.

( a )

 a Fig. 2.38 -Fourier spectra in a zoom of the frequency-time diagram along 10 4 roundtrips for a simulation (a) at θ = -3.03, where the roll regime is stable, and (b) at θ = -3.04 where rolls are unstable. Both spectra are shown in (c), averaged along 5×10 4 roundtrips. The full frequency axis covers a range of (-50,50) GHz.

  Fig. 2.39 -Zoom (a) of the spatiotemporal diagram of a simulation for µ = 1.1, θ = -3.04 and Y = 0.004 centered on an event of high intensity (3). (b) shows the intensity time traces at fixed roundtrip corresponding to the horizontal cuts highlighted on the diagram.

  Fig. 2.40 -(a) PDF of the temporal peak heights (in black), computed over 2.1×10 6 τ c , the green dashed line highlights a Gaussian fit of the initial slope of the distribution. (b) PDF of all the values explored by the intensity (in black), computed over 2.1×10 5 roundtrips. The green dashed line corresponds to a negative exponential fit highlighting the first exponential decay of the data for high values of intensity. The data in dark yellow shows the corresponding PDFs for a simulation run with the same parameters but a different detuning θ = -3.048.

  Fig. 2.42 -Spatiotemporal diagrams and PDFs of peak heights for three different simulations run at θ = -3.042 (a-b), θ = -3.044 (c-d) and θ = -3.048 (e-f).

  Fig. 2.43 -(a) Time trace for D (blue solid line) and its adiabatic value D s = µ/(1 + I) (green solid line) over 10 4 roundtrips. (b) Zoom on the last 6 roundtrips.

Fig. 2 .

 2 Fig.2.44 -Histogram of the values assumed by the electric field in the complex plane (Re(E), Im(E)) in the experimental case (a), computed over 10 4 roundtrips, and in the numerical case (b), computed over 10 5 roundtrips, in order to maintain the same number of events associated to the bounded state.

Fig. 2 .Fig. 2 .

 22 Fig. 2.45 -Experimental phase time trace in the proximity of an abnormally high event in a time span of 100 ns (a) and in a zoom at the beginning of the slope change (b). The initial phase is arbitrarily set to zero and the high-peak event occurs at t = 0.

Fig. 2 .Fig. 2 .Fig. 2 .Fig. 2 .

 2222 Fig. 2.47 -Phase portrait of the temporal profiles I and II, highlighted in Fig. 2.46. (a) displays a bounded regime while (b) shows a 2π excursion of the phase around the origin.

  Fig. 2.51 -Spatiotemporal diagram of the intensity I of the electric field (a) and of the carrier population D (b). The light blue (orange) line corresponds to clockwise (counterclockwise) phase rotations, the abnormal event has been highlighted by the green circle.

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 127 3.2 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 128 3.3 HSS and linear stability analysis . . . . . . . . . . . . . . . . . . . . 129 3.4 CS interaction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 133 3.4.1 Interaction potential and merging time . . . . . . . . . . . . 136 3.4.2 Possible analogy with hydrophobic materials . . . . . . . .

Fig. 3

 3 Fig. 3.1 -Stationary curve as in Eq. (3.5) (blue solid line) and in Eq. (3.4) (red solid line) for the parameters α = 4, β = 0.125, θ = -2, µ = 1.2µ thr . The highlighted fixed value of injection corresponds to E I = 1, used in the simulations of the model in Eqs. (3.2).

2 -

 2 Fig. 3.2 -Stationary instability domain of the model in Eqs. (3.2) for fixed parameters α = 4, β = 0.125, θ = -2, d = 0.052, µ = 1.2µ thr = 1.404.

Fig. 3

 3 Fig. 3.3 -Hopf instability domains for fixed parameters α = 4, β = 0.125, θ = -2, d = 0.052, µ = 1.2µ thr = 1.404 and σ = 1(a), 4(b), 40(c), 400(d).

  fixed parameters in the CS study are α = 4, θ = -2, d = 0.052, E I = 1, µ = 1.2µ thr = 1.404 and β = 0.125. The diffusion coefficient d in Eqs. (3.3) is set to d = 0.01, similarly to the model described in Chapter 1. In Fig. 3.2 and 3.3 we indicated, respectively, the stationary and Hopf instability domains for σ =1(a), 4(b), 40(c) and 400(d) for the model in Eqs. (3.2).

( a )

 a Stationary curve. (b) Stationary instability domain.(c) Hopf instability domain.

Fig. 3

 3 Fig. 3.4 -Stationary and Hopf instability domains for fixed parameters α = 4, θ = -2, µ = 6, d = 0.01 and σ = 400. The value of injection amplitude highlighted by the vertical black dashed line in (a) corresponds to E I = 0.5.

  .5 for the parameters α = 4, θ = -2, d = 0.052, E I = 1, µ = 1.2µ thr = 1.404, β = 0.125 and σ = 400. As usual the dynamical equations are integrated through a Fourier split-step method with (a) t = 0 ns. (b) t = 0 ns. (c) t = 1154 ns. (d) t = 1154 ns. (e) t = 1158 ns. (f) t = 1158 ns. (g) t = 1160 ns. (h) t = 1160 ns.

Fig. 3

 3 Fig. 3.5 -CS interaction process for the parameters α = 4, θ = -2, d = 0.052, E I = 1, µ = 1.2µ thr = 1.404, β = 0.125 and σ = 400. The initial distance between the two CSs is d 0 = 20 spatial units. On the left side we illustrated the CSs in the transverse plane of the electric field intensity, on the right side we depicted the intensity profile along the y-axis at x = 0.

  Fig. 3.7 -Merging time of two CSs, in logarithmic scale, as a function of their initial distance d 0 /2 for σ = 400. The black solid line represents a fit y = a + bx, with a = -3.460 ± 0.12 and b = 1.055 ± 0.01, computed on the eight largest initial distances, from d 0 /2 = 10 to d 0 /2 = 13.5.

  Fig.3.8 -Temporal evolution of the distance between two CSs for σ = 400 and and initial distance r 0 = 10 (same simulation as in Fig.3.5). The black solid line indicates the simulation results, while the black dashed line depicts the analytic curve given by Eq. (3.7) with the parameters R = 0.474 and K = 48.77.

  Fig. 3.9 -Merging time of two CSs (in logarithmic scale) as a function of their initial distance for different values of the relative decay rate σ.

Fig. 3 .

 3 Fig. 3.10 -Merging time of two CSs as a function of the inverse of the relative decay rate σ for different values of the initial distance r 0 .

Fig. 3 .

 3 Fig. 3.11 -Example of the chaotic regime observed in the transverse plane of the electric field intensity for low values of optical injection in the simulations of Eqs. (3.3). Here in particular we considered the parameters σ = 400, θ = -2, µ = 1.404, α = 4, d = 0.01 and a value of injection E I = 0.5 below the locking point.

  (a) Total intensity PDF.(b) Spatiotemporal maxima PDF.

Fig. 3 .

 3 Fig. 3.12 -PDFs of total intensity (a) and spatiotemporal maxima (b) for a simulation with the parameters σ = 400, α = 4, d = 0.01, µ = 6 and E I = 0.5. The light blue solid line in (a) is a negative exponential PDF computed on the average of the data as exp(-I/ I )/ I . The three vertical dashed lines are three extreme event thresholds and the green and magenta solid lines in (b) are, respectively, a Gumbel and a Weibull PDF computed through the average and the standard deviation of the data (the definition of these PDFs and of the three extreme event thresholds are given in Section 1.6 in Chapter 1).

Fig. 3

 3 Fig. 3.13 -Example of an extreme event observed in the transverse plane of the electric field intensity for the parameters σ = 400, α = 4, d = 0.01, µ = 6 and E I = 0.5. The dark red cross highlights the extreme event.

  (a) Electric field intensity. (b) Phase of the electric field. (c) Carrier density D.

Fig. 3 .

 3 Fig. 3.14 -Color maps of (a) the electric field intensity, (b) the phase of the electric field and (c) the carrier density in presence of the extreme event illustrated in Fig. 3.13. The dark red cross corresponds to the spot where the extreme event occurs.

Fig. 3 .

 3 Fig.3.15 -Trajectory (in time) of the system in the (D, φ, I) phase space in the spot highlighted in Fig.3.14 by a dark red cross, where an extreme event occurs. Time evolves along the trajectory from green to yellow. The blue line highlights the fixed point C (unstable focus), whose phase has a 2π modulus.

Fig. 3 .

 3 Fig. 3.16 -(a) Zoom of the phase of the electric field around the extreme event in Figs. 3.13, 3.14. (b) zero value isolines for the real (blue solid line) and imaginary (red dotted line) part of the electric field zoomed in the same area as (a): the point of crossing between the two different isolines correspond to optical vortices. The dark red cross corresponds to the spot where the extreme event occurs.

P

  (x) = ẋ[ ẍ + ω 2 LSA (e x -1)] = ẋ ẍ + dV LSA (x) dx =ẋ [2b (D s n + Ds r n) + 2b (D s n + Ds rs n) I s e x ] Ds r 2 s ẋ2 + D s + Ds r 2 s 2 I s ẋ2 e x ,where we have used the expressions for n and n as functions of ẋ. The time derivative of the mean energy on a period 2t max is then given by Ḣ dt + D s + Ds r 2 s 2 I

  Fig. E.1 -(a) Relaxation oscillations in a class B laser; (b) comparison between the periodsT 1 e T 2 given by the numerical integration of Eq. E.7 (in blue and red, respectively) and the analytic approximation for large energies (in cyan and green, respectively).

Fig. E. 2 -

 2 Fig. E.2 -Comparison between the values of T 1,2 /T 0 predicted from the Toda oscillator model and the data measured during the full system evolution. The simulation was done for the following set of parameters: s = 3, b = 10 -4 , r = 1, µ = µ h + 0.5 ≈ 3 (a), µ = µ h -0.5 ≈ 2 (b), and initial conditions (I 0 , n 0 , n0 ) = (0.001, 0, 0). The system has been numerically integrated on 10 6 temporal units. We indicated with a cross the data related to the solution at regime in (b). (c) Same comparison for the same parameters for 30 simulations with 30 different values of a, always below the Hopf bifurcation.

F

  η , τ)e ik η η dk η = ik η F +∞ -∞ Ê(k η , τ)e ik η η dk η = ik η Ê(k η , τ) η , τ)e ik η η dk η = -k 2 η Ê(k η , τ) then we have ∂ Ê(k η , τ) ∂τ = -k η i + k η d Ê(k η , τ) , (F.2)The differential Eq. (F.2) can be solved through separation of variables and, unless constants, we obtain∂ Ê(k η , τ) Ê(k η , τ) = -k η i + k η d ∂τ log Ê(k η , τ) = -k η i + k η d τ Ê(k η , t) = e -k η (i+kηd)τ .

F

  ∂D(x, y, t) ∂t = ∂ D(k x , k y , t) ∂t = -dk 2 D(k x , k y , t) , (G.5)The two differential Eqs. (G.4) and (G.5) can be solved through separation of variables and, unless constants, we obtain, for Eq. (G.4),∂ Ê(k x , k y , t) Ê(k x , k y , t) = -iσk 2 ∂t log Ê(k x , k y , t) = -iσk 2 t Ê(k x , k y , t) = e -iσk 2 t ,and, for Eq. (G.5),∂ D(k x , k y , t) D(k x , k y , t) = -dk 2 ∂t log D(k x , k y , t) = -dk 2 t D(k x , k y , t) = e -dk 2 t .

  ] . . . . . . . . . . 4 Rogue wave formation in [Liu 2015] . . . . . . . . . . . . . . . . . 5 Extreme event dynamics in [Zamora-Munt 2013] . . . . . . . . . . 6

  2.33 Roll pattern in the spatiotemporal diagram and phase space . . . 2.34 HSS and roll pattern solution . . . . . . . . . . . . . . . . . . . . . 2.35 Roll pattern breaking for decreasing values of injection . . . . . . 2.36 Experimental spatiotemporal diagram in the case of abnormally high events emerging from an unstable roll background . . . . . . 2.37 HSS for the intensity as a function of the detuning . . . . . . . . . 2.38 Fourier spectrum for a stable and an unstable roll pattern . . . . . 2.39 Spatiotemporal diagram and temporal profiles around an abnormal event . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.40 PDFs of the peak height and total intensity . . . . . . . . . . . . . 2.41 Experimental histogram of the peak heights and PDF of the total intensity . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.42 Spatiotemporal diagram and PDF of peak height for different values of detuning . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.43 Time trace of the carrier density and its adiabatic value . . . . . . 2.44 Experimental and numerical histograms for the electric field in the complex plane . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.45 Experimental phase time trace . . . . . . . . . . . . . . . . . . . . . 2.46 Spatiotemporal diagram and temporal profiles around an abnormal event as seen in the experiment . . . . . . . . . . . . . . . . . 2.47 Experimental phase portraits when approaching an abnormal

1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 1.2 Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 1.3 HSSs and linear stability analysis . . . . . . . . . . . . . . . . . . . 27

	1.3.1 Bistability . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
	1.3.2 Plane-wave instability . . . . . . . . . . . . . . . . . . . . . . 30
	1.3.3 Pattern-forming instabilities . . . . . . . . . . . . . . . . . . . 32
	1

.4 Zoology of the stationary and nonstationary solutions . . . . . . . 33

	1.4.1 Cavity solitons . . . . . . . . . . . . . . . . . . . . . . . . . . 34
	1.4.2 Extended spatiotemporal chaos . . . . . . . . . . . . . . . . . 35
	1.4.3 Chaotic cavity solitons . . . . . . . . . . . . . . . . . . . . . . 38
	1.4.4 Oscillatory cavity solitons . . . . . . . . . . . . . . . . . . . . 40
	1.4.5 Role of the diffusion coefficient d . . . . . . . . . . . . . . . . 42
	1

.5 Extreme event detection . . . . . . . . . . . . . . . . . . . . . . . . . 43

	1.5.1 Maxima individuation method . . . . . . . . . . . . . . . . . 43
	1.

6 Statistical analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45

	1.7
	1.6.1 PDF of the spatiotemporal maxima . . . . . . . . . . . . . . . 45
	1.6.2 Rayleigh and negative exponential PDF . . . . . . . . . . . . 48
	1.6.3 Weibull PDF . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
	1.6.4 Gumbel PDF . . . . . . . . . . . . . . . . . . . . . . . . . . . 51

Dependence on the laser parameters . . . . . . . . . . . . . . . . . 53 1.8 Comparison with other works in the field . . . . . . . . . . . . . . 57 1.9 Waiting time statistics . . . . . . . . . . . . . . . . . . . . . . . . . . 60 1.10 Profiles and cavity soliton comparison . . . . . . . . . . . . . . . . 61 1.11 Future perspectives . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64

	1.11.1 Conservative limit of a LSA . . . . . . . . . . . . . . . . . . . 64
	1

.12 Conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

  2with transparency currents I 0,1;2 = eV 1;2 N 0,1;2 /τ nr,1;2 (1 + ζ 1;2 ). Note that in the case studied in this Thesis the absorber is not pumped hence I 2 = 0. Finally b, r and s in Eqs. (1.1) are defined as decay rates redefined as τ 1;2 = τ nr,1;2 /(1 + 2ζ 1;2 ).F(k x , k y , t)e i(k x x+k y y) dk x dk y F(k x , k y , t)e i(k x x+k y y) dk x dk y = -k 2 x F(k x , k y , t)and, in the same manner,F ∂ 2 F(x, y, t) ∂y 2 = -k 2 y F(k x , k y , t) ,then we have∂ F(k x , k y , t) ∂t = -k 2 (d + i) F(k x , k y , t) , (B.4) with k 2 = k 2 x + k 2 y .The differential Eq. (B.4) can be solved through separation of variables and, unless constants, we obtain

	Since
	F	∂ 2 F(x, y, t) ∂x 2		= F	∂ 2 ∂x 2	+∞ -∞
				= -k 2 x F	+∞ -∞
		b =	2τ p τ 1	,	r =	τ 1 τ 2	,	s =	a 2 b 1 a 1 b 2
	F	∂F(x, y, t) ∂t	= (d + i) F	∂ 2 F(x, y, t) ∂x 2	+ F	∂ 2 F(x, y, t) ∂y 2
								∂x 2	+ F	∂ 2 F(x, y, t) ∂y 2	.

with ∂ F(k x , k y , t) ∂t = (d + i) F ∂ 2 F(x, y, t)

  t w , t w + t) =Sibani 2003], given α sequences of recordings, the probability of having k events in the time span between t w and t w + t is the following P α,k (t w , t w + t) =

	1 k!	t w + t t w	-1	log	t w + t t w	k	,
	furthermore [1 k!	t w + t t w	-α	α log	t w + t t w		k	,	(D.1)

  α,k (t w , t w + t) = 1 k! (e) -αx(t w ,t) [αx(t w , t)] k ∼ e -αx(t w ,t) with null integration constants, since2 (D s n + Ds n) = 2b (D s + s Ds r n) I s

	∼	(αe) k k k √ 2πk LSA x(t ω 2	= ẋ .

w , t) k e -αx(t w ,t) (D.2)

  t 0 < t < t 0 x max + log sech 2 ω LSA e x max 2 (tt max ) t 0 < t < -t 0 + 2t maxwhere we defined t 0 = t| x=0 and t max = t| x=x max as follows e -x max neglecting terms of order e -exp (x max ) . These analytic expression can also be generalized for all periods,

	    	-e x max +	ω 2 LSA t 2 2
	   
		t 0 = t max =	√ ω LSA 2e x max √ 2e x max ω LSA	+	1 ω LSA	2 e x max	arctanh	√	1 -e -x max
		= t 0 +      -e x max + ω LSA 1 ω 2 LSA (t -2mt max ) 2 2 e x max arctanh 2 1 x(t) = √     x max + log sech 2 ω LSA e x max 2 [t -(2m + 1)t max ]

x(t) =

  Ê(k x , k y , t)e i(k x x+k y y) dk x dk yÊ(k x , k y , t)e i(k x x+k y y) dk x dk y = -k 2 x Ê(k x , k y , t)and, in the same manner,

	F	∂E(x, y, t) ∂t	= iσ F	∂ 2 E(x, y, t) ∂x 2	+ F	∂ 2 E(x, y, t) ∂y 2
		∂ Ê(k x , k y , t) ∂t	= iσ F	∂ 2 E(x, y, t) ∂x 2	+ F	∂ 2 E(x, y, t) ∂y 2	.
	Since
	F	∂ 2 E(x, y, t) ∂x 2	= F	∂ 2 ∂x 2	+∞ -∞
			= -k 2 x F		+∞ -∞
		F	∂ 2 E(x, y, t) ∂y 2		= -k 2

y Ê(k x , k y , t) , then we have ∂ Ê(k x , k y , t) ∂t = -iσk 2 Ê(k x , k y , t) , (G.4) with k 2 = k 2 x + k 2 y .

Similarly for Eq. (G.3b) we obtain

Cristina Rimoldi Valbonne, FranceDecember 2017

and Ê(k x , k y , t + ∆t) = Ê(k x , k y , t)e -iσk 2 ∆t , D(k x , k y , t + ∆t) = D(k x , k y , t)e -dk 2 ∆t .

In our programs we implement the split-step method executing first a half step ∆t/2 in the Fourier space, then a full step ∆t of the purely temporal part of the equations through a Runge-Kutta method at second order and then completing with another half step in the Fourier space.

Acknowledgments

No achievement in someone's life is ever reached alone. This is no different. There are many people to thank that helped me reach this point where a chapter of my life ends, giving way to new experiences and adventures.

First of all I would like to thank the people that made it possible for me to follow my dreams and supported me every step of the way: my parents and my family. Thanks for always encouraging me to follow my interests in my study choices as well as in life experiences and for always being there for me, in a bond of trust and love. Secondly I would like to thank my supervisor for following me along this Thesis path, for always stimulating curiosity and interest for new topics in the people around her through many remarks as well as new questions and explanations, for her availability, her patience and very friendly advice, well beyond the obligations of a supervisor, aimed to enrich my formation. Many thanks also go to my co-supervisor, who followed me since the beginning of my university experience, especially during the master years, and has always being available to discuss and advice me on my formation and the topics of this Thesis. Then I would like to thanks my fellow PhD students and friends at the lab, whether they are still here or they left for new places, you know who you are. Thanks for making this experience so much funnier, interesting and overall amazing with your jokes, your ideas and your attitude toward work and life in general. Thanks for all the chats and discussions whether it was over a barbecue, a trip or a conference together. Many thanks also to the administrative personnel and the other researchers of the lab for making it a nice place to work in, which is not something that goes for granted. Finally I would like to thank all the friends in my hometown for always being available to meet up during the weekends, for a chat or a drink or both, and for the traditional amazing holiday hikes together.

Many thanks basically to all the people who shared with me a piece of the road with their support and their friendship.

Appendix A

LSA equations in physical variables

In terms of the physical electric field E and the carrier densities in the active (N 1 ) and passive (N 2 ) medium, Eqs. (1.1) can be written [Wilmsen 2001, Carr 2001] as:

where [Vahed 2012] τ p is the photon lifetime, a 1;2 is the differential gain, N 0,1;2 is the transparency carrier density, L 1;2 is the effective length of each region, α and β are the linewidth enhancement factors and the cavity losses, such as the mirror losses and the linear absorption, are included through T. I 1;2 is the injected current, e is the electron charge and V 1;2 is the volume of, respectively, the gain and the absorption region, τ nr,1;2 is the non-radiative recombination time, B sp,1;2 is the coefficient of bimolecular radiative recombination, v g is the group velocity and N p = 0 n 2 |E| 2 / (2hω) is the photon density with n the background refractive index and hω the photon energy).

Let us now define η 1 as the ratio of gain to losses and η 2 as the ratio of absorption to losses. Furthermore we introduce ζ 1;2 as the ratio of non-radiative to radiative decay times respectively in the active and passive media.

In this way we can observe that the radiative recombination parameter B that appears in Eqs. (1.1) corresponds to

Split-step method for Eqs. (1.1)

Let us consider the following model, same as in Eqs. (1.1),

The split-step method consists here in separating the purely temporal part of the evolution of the dynamical equations from the spatial part. Furthermore the spatial part of the evolution is computed in the Fourier space. In this case the temporal part is the following

which is integrated through a Runge-Kutta method at the second order, and the spatial part is contained just in the differential equation for

Now, the Fourier transform of F is defined as follows

F(x, y, t)e -i(k x x+k y y) dxdy with F Fourier transform operator, and its inverse transform is

Applying the operator F to Eq. (B.3), we obtain

Appendix C

The Routh-Hurwitz stability criterion

We will give here proof of the criterion for third and fourth degree polynomial.

Let us first consider the following characteristic equation:

where J is a 3 × 3 Jacobian matrix for a set of three differential equations describing some kind of system. Then we need to find the roots of

We are interested in studying where the system passes from a stable to an unstable domain, which occurs when the real part of the eigenvalues λ from negative becomes positive. On the boundary of the instability domain λ results purely imaginary, we can then substitute λ = iν in Eq. C.2 and obtain

If λ is real then ν = 0 and the stationary instability boundary is given by

in particular if c 0 < 0 the system is stationary unstable. On the other hand if λ is complex then the Hopf instability boundary is given by

where we have separated real and imaginary part of Eq. (C.3), and we therefore obtain

in particular if c 2 c 1 < c 0 then the system is Hopf unstable.

The Routh-Hurwitz stability criterion

Let us now consider Eq. (C.1) where J is now a 4 × 4 Jacobian matrix for a set of four differential equations describing some kind of system. Then we need to find the roots of:

Proceeding in the same manner as in the previous case, we substitute λ = iν in Eq. (C.4), obtaining

Once again the stationary instability boundary is given by the condition c 0 = 0 in particular if c 0 < 0 the system is stationary unstable.

Regarding instead the Hopf instability boundary, we separate real and imaginary part of Eq. (C.5) and obtain

Appendix D

Log-Poissonian distribution

Let us consider [Sibani 1993, Sibani 2003] a discrete time sequence of events (for instance rogue waves) described by identically distributed positive stochastic variables N i , with i = 1, ..., t. A record occurs at a time j if and only if the variable N j is greater than any N i previous variable. We immediately point out that with this definition the first record will always happen at i = 1. Let us derive now the probability P n (t) that n records occur in a time between 0 and t.

It is clear that the probability of one and only one record at time 0 < m < t is

therefore the probability of having exactly two records, the first at i = 1 and the second at i = m is

.

We can then observe that the probability of having the first record at i = 1 and the second record at any subsequent time m is given by the sum over all possible m values and in the limit t 1 we obtain [Sibani 1993]

Now the probability of having k records, the first at i = 1 and the remaining k -1 exactly at times 1 < m

therefore the probability of having k records at any time between 0 and t is

LSA conservative limit

Large energy analytic approximation in the conservative limit

In Chapter 1 we observed that it is possible to describe the system of a purely temporal (semiconductor) laser with an intracavity saturable absorber as a particle of unitary mass subject to a Toda potential. In particular the conservative limit of this system is given by

the Toda potential, and

Eq. (E.2) can be rewritten as ẍ = -ω 2 LSA (e x -1) which implies, from Eqs. (E.1)

Further notes

It is possible to reduce the number of the equations of the conservative LSA system from three to two: in fact starting from

we can notice that s ṅ -ṅ/r = 0 hence snn/r is a constant of the system.

Therefore if we substitute z 1 = sn + n/r, z 2 = snn/r we obtain the following system

If we then define as initial conditions z 1,0 = 0, z 2,0 = 0, which is n 0 = 0 and n0 = 0, we can apply the following reduction

where we have normalized the variable z 1 to 2sI s .

If we insert optical injection into the model and we apply the same initial conditions to the system in its reversible limit it is possible once again to reduce the number of real variables, here from four to three. Interestingly enough, through some substitutions and rescaling we can formally reach the same set of equations studied in [Politi 1986] for a class B laser with optical injection

What was found in this system was a coexistence of conservative and dissipative behavior for certain values of the parameters à and f (in [Politi 1986] for à = 4/3 and f = 0.42). We can then conclude that also in the case of a laser with saturable absorber, the optical injection destroys the conservative limit of the laser but preserves its reversible limit, giving rise to the same coexistence of conservative and dissipative structures for specific sets of parameters and different initial conditions of the model under study.

Appendix G

Split-step method for Eqs. (3.2)

Let us consider the following model, same as in Eqs. (3.2),

As in Appendix B we separate the purely temporal part of the dynamical equation evolution from the spatial part, to be computed in the Fourier space.

In this case the temporal part is the following

which is integrated through a Runge-Kutta method at the second order, and the spatial part is given by

The Fourier transforms of E and D are defined as follows

List of acronyms