Laser chimeras as a paradigm for multistable patterns in complex systems

A chimera state is a rich and fascinating class of self-organized solutions developed in high-dimensional networks. Necessary features of the network for the emergence of such complex but structured motions are non-local and symmetry breaking coupling. An accurate understanding of chimera states is expected to bring important insights on deterministic mechanism occurring in many structurally similar high-dimensional dynamics such as living systems, brain operation principles and even turbulence in hydrodynamics. Here we report on a powerful and highly controllable experiment based on an optoelectronic delayed feedback applied to a wavelength tuneable semiconductor laser, with which a wide variety of chimera patterns can be accurately investigated and interpreted. We uncover a cascade of higher-order chimeras as a pattern transition from N to N+1 clusters of chaoticity. Finally, we follow visually, as the gain increases, how chimera state is gradually destroyed on the way to apparent turbulence-like system behaviour.


Supplementary Note 1: Experimental setup
It consists mainly of 5 functional physical elements: 1. The tunable laser diode is expected to provide a laser light beam which wavelength can be continuousy and linearly tuned through an electrical drive current. It is a linear electrical-to-optical converter, more precisely a current-to-laser wavelength converter.
2. The Fabry-Pérot interferometer is ruling the nonlinear transformation occurring from the wavelength deviation, into optical intensity fluctuations. It is a "nonlinear" opticalto-optical converter. Though it is commonly considered as a linear optical filter for an incoming optical wave amplitude, it has to be considered here as nonlinear in our delay dynamics problem, since the input variable is the instantaneous wavelength of laser, and not the electromagnetic field amplitude (the output being the optical intensity).
3. The photodiode is required since we expect to have consistent physical quantities in such a closed feedback architecture: we need to recover back an electrical signal before closing the oscillation loop. This is fulfilled by the photodiode, converting linearly, and instantaneously the incoming optical intensity into a detected photo-current (photodiode bandwidth, as well as laser tuning bandwidth, are far above the characteristic times of the next elements of the electrical feedback, typically beyond 10 MHz).
4. An electrical digital delay line is providing the delay in the feedback, moreover in a very flexible way. The delay can be indeed easily adjusted whether through the depth of the FIFO (First In First Out) memory used to generate the delay, or through the digital clock frequency (0.6 MHz to 40 MHz) ruling the speed at which the signal is travelling through the memory. This bandwidth (≃ 120 kHz), though not that important, is also assumed to be negligeable (instantaneous delay of any relevant Fourier frequency component of the travelling signal) compared to the last element.
5. An electronic filter is finally used to force the dominant time scales ruling the dynamics of the oscillator loop. The filter is of bandpass type, with the simplest second order profile. A (variable) gain function can also be attached to this last element, allowing for the tuning of the overall closed loop gain of the nonlinear delay oscillator. This variable gain is simply provided by a multiplier circuit, one input being the signal to amplify, and the other being a tunable constant voltage controlling thus linearly the normalized feedback β around the values of interest (typ. from 0 to 5).

Supplementary Note 2: Modeling
The wavelength tunable laser This element is the central one in a wavelength dynamics implemented through an optoelectronic feedback loop. It allows for a physical variable conversion from an electrical signal (a current) into a laser wavelength deviation. The device is a DBR (distributed Bragg reflector) multi-electrode laser diode as depicted in Fig. 1. It consists of a Fabry-Pérot cavity having natural reflector on one side, and a tunable DBR reflector on the other side. The latter reflector has a maximum Bragg reflection wavelength determined by the optical periodicity nΛ B of a Bragg reflector, which can be tuned through the change of the refractive index of the Bragg material. This variation is caused by the modulation of the free carrier occupying the diode junction where the Bragg grating is designed, through the modulation of a forward injection current I B . The lasing wavelength is thus directly tuned according to the maximum of reflection of the DBR grating. The laser gain is provided by a gain section next to (but electrically indepedent from) the Bragg section, through an electrical pumping current I A (lasing threshold typ. around 8 mA, operating value at 20 mA, leading to a useful lightbeam of ca. 1 mW optical power).
For a DBR current modulated around 11 mA (λ 0 ≃ 1549.8 nm), the device exhibits a continuous tuning range with an efficiency of about S λ ≃ 0.19 nm/mA. This continuous tuning range spans over 1.45 nm, thus corresponding to a peak to peak current modulation amplitude of ca. 7 mA (see Fig. 1 Right).
If one assumes the current I B is delivered by a voltage source v N = v B + α v in series with a "strong" (Norton model) resistor R N (much stronger than the static resistance of the forwardly biaised DBR junction), so that the current I B can be confidently determined as: where v B /R N is the DC contribution to the junction current bias < I B > (around 11 mA), such that the other v−contribution to I B verifies < v >= 0. α is a voltage-linearly adjustable factor used to tune the overall gain in the oscillator loop.
Notice finally that the laser delivers a variable wavelength light beam, with an optical power P 0 . This laser beam is passed then through an optical filter in order to perform a nonlinear transformation of the input wavelength into the output light intensity.

The nonlinear transformation
Contrarily to the original wavelength chaos generator reported in [1], the observation of chimera states requires two important differences. One of these relates to an asymmetric shape of the nonlinear transformation involved in the delay dynamics. With the usual two-wave interference provided e.g. by a birefringent interferometer, one indeed observes symmetric maxima and minima described analytically by the well known cos 2 or sin 2 nonlinear profile.
Since the wavelength delay dynamics involves a nonlinear transformation from a wavelength variation into the intensity output of a spectral optical filter, a rather straightforward solution to obtain asymmetry between minima and maxima is to make use of a multiple wave interferometer, such as a Fabry-Pérot one: where 2ne is the optical path corresponding to one round trip in the cavity, m = 4R/(1−R) 2 is a function of the reflection coefficient in intensity R of each mirror end of the cavity influencing the "contrast" ([(1 + R)/(1 − R)] 2 ) of the Fabry-Pérot transfer function, A = 1/(1 − R) is a vertical scaling factor, and Φ(t) = 2πne/λ(t) is a phase modulated due to the variation of the laser wavelength λ(t) = λ 0 + δλ(t). Considering the very small deviation of the wavelength around a central value The fine tuning of the central laser wavelength (through a mean DBR and/or active laser current, and/or through the controlled laser temperature) allows to adjust the actual rest point of the Fabry-Pérot modulation transfer function around which the modulation operates (e.g. positive or negative slope, around a broad minimum or a sharp maximum of the nonlinear function formed by this wavelength-to-intensity conversion through the Fabry-Pérot effect).
One could notice here the conceptual analogy between the FM delay dynamics [2,3], and the wavelength delay dynamics [1]. The voltage controlled oscillator in the first case corresponds to the tunable laser in the second. The resonant RLC filters which filtering profile determines the nonlinear function in the first case, corresponds to the optical filtering function provided by the two-wave or multiple-wave interference in the second case.
Physical values of the Fabry-Pérot parameters have been properly chosen to exhibit a very few broad resonances and anti-resonances within the wavelength tuning range of the laser. The Fabry-Pérot was designed with a glass plate of e = 5 mm thickness, and intensity reflection coefficients R = 50 % (20 GHz free spectral range, or 0.16 nm around 1.55 µm).
The Fabry-Pérot output intensity fluctuations are pratically converted into an electrical signal by a simple photodiode, thus providing an electrical signal which can be used for the electronic feedback onto the DBR electrode of the tunable laser: where S = 0.9 A/W is the photodiode sensitivity, and R is the equivalent resistor of a typical transimpedance amplifier typically used to convert the photocurrent into a voltage. The product R S (in V/W) is then the conversion efficiency of an integrated amplified photodiode.

The delay
The natural physical delay line with an optoelectronic setup is obviously a fiber spool. The magnitude of the delay is however to be compared with the slowest characteristic response time involved in the feedback loop. In the present setup, we are planning to observe dynamics in the µs range, thus implying much greater time delays (up to ms Special attention needs however to be drawn on the possible response times τ which can be set, due to basic sampling theory limitations (the sampling rate being defined by f CLK , thus imposing a maximum non-zero Fourier component up to f CLK /2).
The Fourier frequency filter ruling the integro-differential dynamics The motion of the feedback loop oscillator, differently speaking the differential equation ruling the dynamics of the designed nonlinear delayed oscillator, is originating from the limiting Fourier filtering performed by the oscillation loop in the linear regime. More precisely, it has been shown in [2] that observation of virtual chimera motions requires a bandpass filtering, which results in an integro-differential process in terms of equation of motion.
Assuming a bandpass filtering of the simplest form, the corresponding Fourier filter is trivially expressed as: where τ = (2πf h ) −1 is the high cut-off frequency of the corresponding low pass filter, and θ = (2πf l ) −1 is the low cut-off frequency of the corresponding high-pass filter, and G 0 is the gain of the filter. Equation (6)-Left is a Fourier domain second order bandpass filter, which illustrates the global frequency response of the filter (relative phase shift and attenuation between output and input, when a sinusoidal signal is used at the input). Equation (6)-Right is a local (in time) description through a differential law ruling the evolution of the filter output v at time t, knowing the input signal u D . The global description for v in the time domain can be also expressed, as a convolution (integral) operation involving the so-called impulse reponse of the filter h(t), practically defined as the inverse Fourier transform of the Fourier domain filter transfer function in Eq. (6)-Left: where h H (t) is the stepwise Heavyside function (null for t < 0 and unity for t ≥ 0), reflecting the causal feature of any impulse response. The impulse response of the second order bandpass filter can be analyzed as follows: • At the origine, it has a stepwise variation from zero to τ −1 , identically to the first order low pass filter with the same characteristic time τ (low-pass case: H(ω) = 1/(1 + iωτ ) and FT −1 [H(ω)] = h(t) = τ −1 e −t/τ h H (t)); • From the origine, the function is decreasing linearly with a slope −τ −2 (1 + τ /θ) (which is steeper than the low pass one, −τ −2 ); • It crosses zero at time t 0 = [τ · ln(θ/τ )]/(1 − τ /θ); • It reaches a minimum at 2t 0 , with an amplitude close to −θ −1 (1 − α) where α is a small quantity scaling as 2(τ /θ)[ln(θ/τ )]/(1 − τ /θ); • It finally decreases slowly exponentially to zero from negative values as (−θ −1 )e −t/θ /(1 − τ /θ).
Since τ ≪ θ in the definition of the bandpass filter, it is worth noticing that (τ /θ) ≪ 1 is a very small quantity (of the order of 10 −4 ).
It is worth noticing here the main difference in terms of impulse response, between a low pass filter delay dynamics (the usual in most of the popular delay dynamics, e.g. Mackeygalss or Ikeda ones), and the bandpass of highest importance for the observation of sustained chimeras patterns. The previous description of the bandpass impulse response reveals nonmonotonus feature of its "Mexican hat"-like profile through the negative values it takes, on the contrary to the usual impulse response of the low pass case which is indeed monotonuously decreasing to zero.

Supplementary Note 3: Delay dynamics normalization
A very common way for the time normalization in delay equations is to introduce a new time variable s = t/τ D (dt = τ D ds) normalized to the delay. The differential equation (6), once re-written in s and combined with the time domain delay process (5), reads as: where δ = τ D /θ and ε = τ /τ D are two (generally small in the large delay and broadband situation of concern) temporal parameters, which are crucial for the stability of the virtual chimera states. This is precisely one of the purpose in the present article. The experimental setup is indeed intended to explore physically the influence of these temporal parameters on the chimera states stability. Numerical simulations as well as analytical investigations, have estimated the practical range of interest for the scanning of different chimera patterns in the experiment.
When combining Eqs. (1) to (4) together with Eq. (8), one is led to the amplitude normalization of the dynamics, considering the various conversion efficiencies and gains in the optoelectronic feedback loop. The unit free variable x(s) can be defined as the varying argument δΦ inside the sin 2 −function appearing in Eq. (3). In order to convert v into δΦ in the left hand side of Eq. (8), one needs to multiply both sides of the differential equation with the constant factor (α S λ /R N ) · (2πne/λ 2 0 ) (where α(U G ) = U G /U G 0 is the linearly tuneable gain already described in Eq.(2)) leading to: As a direct consequence, two amplitude parameters are revealed through this normalization, the feedback loop gain β and the offset phase Φ 0 ruling the mean operating point for the Fabry-Pérot nonlinear modulation transfer function, One should notice in these last two normalized expressions with respect to physical parameters, that both appear as independently adjustable experimentally: • β through the voltage controlled factor α(U G ), • and Φ 0 through the offset voltage v B .
The normalized dynamics written as a vectorial (2 variables, x and y = x ds) delay differential equation reads finally as follows: dy ds (s) = x(s).
One could slightly simplify the previous equation when neglecting the quantity δε = τ /θ ≪ 1.

Supplementary Note 4: Numerical simulations
The previous normalized model was numerically investigated using a fourth order Runge-Kutta algorithm with constant integration time step, as described in the "Methods" section at the end of the main article.
Despite the fact that experiments obviously showed sustained patterns and thus not transient chimera states (10 6 time delays are indeed corresponding to easily testable duration in practice, with time duration of the order of hours free running experiment with fixed parameter), it is important for the numerical viewpoint to also indicate how long, in terms of number of time delays integration duration, one has to consider the numerics in order to get the asymptotic domains in the (ε, δ)−plane. The actually required time for chimeras to emerge spontaneously from random initial conditions, is practically less than 10 2 time delays. Thus, one could in principle choose a calculation slightly longer than this typical duration. For a safe calcultation margin ensuring the removal of any transient, our numerical simulations have systematically considered a duration of 10 5 . This calculation duration was used to obtain the borders of the domains reported in Fig. 2 in the main article. Calculation duration values of about 10 3 to 10 4 time delays as in Ref. [1] of the main article are however already enough to check that the obtained solution indeed numerically persists at any point of the tested parameter plane.
One could finally also remark that small quantitative disagreement can be observed in the upper-left part of Fig.2 in the main article. This has to be discussed however with respect to the experimental difficulty to accurately calibrate the parameters in this particular region of the hyperbola curves. Qualitative disagreement however only holds when one consider two neighboring experimental hyperbola curves. When staying on the same hyperbola (as it is actually the case in the experimental scanning conditions), the agreement between numerics and experiments is indeed very good (the sequence of successive number of possible heads). Changing the hyperbola constant (jumping from one hyperbola to another) is unfortunately rather critical experimentally in this upper-left part of the (ε, δ)-plane. The consequence is that the relative positions of the hyperbola might not be controlled with high enough accuracy to get "perfect" agreement. Good qualitative is nevertheless globally obtained with a very good accuracy, confirming that the model is well capturing the experimental observations.

Supplementary Note 5: Analytics highlighting the space-time analogy
The dynamics, instead of a local differential equation as in (6) or (8), can be expressed as a global convolution product as already proposed in Eq. (7) when describing the impulse response h(t) of the loop filtering as the inverse Fourier transform of filter transfer function H(ω). Considering the dimensionless time unit, one has: The principle of the space time analogy is based on the mutiple time scale character of the delay dynamics, leading to the definition of the normalized time s as being decomposed into a short time scale contribution of the order of ε = τ /τ D , and into a long time scale contribution of the order of the unit delay: where n ∈ N counts roughly the number of time delay interval since the origine of time (since γ = o(ε) ≪ 1), and σ ∈ [0; η] is a small time variable identifying the fine temporal position within one approximate time delay interval. σ is then represented as a virtual space coordinate spanning over a finite virtual space interval, namely [0; η].
With a simple change of variable ζ = ξ + η − 1 − nη = ξ + γ − nη, one can rewrite the spatial coupling integral as: With the previous expression for I n (σ), this contribution to the dynamics appears as a nonlinear non local "spatial" coupling, according to a nonlinear function defined as f [x n−1 (ζ)], and a "spatially" extended coupling strength with a spatial extension profile determined by the time reversed (due to the convolution integral) impulse response h(s), as discussed for the modeling issues (Suppl. Note 2). The spatial non locallity depends on the actual spread of the impulse response. The local coupling at σ thus amounts to 1 + h(γ), and the spatial extension around the central position σ is spread from (σ − ∆) to (σ + γ), with coupling weights of h(0) and h(∆), respectively. The extension ∆ (≪ 1, thus fitting within the integral domain from σ − 1 to σ + γ) is an arbitrary distance (or short time) defined by the features of the chosen h(s), above which the impulse response has a negligeable coupling amplitude. This has been analytically derived in the modeling section (Suppl. Note 2), Eq. (7), for a second order bandpass filter. It was moreover shown there, that such a bandpass filter is responsible for the integral term in the integro-differential delay equation of motion, and thus is responsible for some specific stability features of the positive feedback delay dynamics, compared to the usual "differentialonly" delay equation (e.g. Mackey-Glass or Ikeda models). In terms of nonlocal spatial coupling, one realizes that such bandpass feature is equivalent to an increase by a factor of ln(θ/τ ) = − ln(εδ) of the original extension ∆, compared to the differential-only situation.