Attractor hopping between polarization dynamical states in a vertical-cavity surface-emitting laser subject to parallel optical injection

Florian Denis-le Coarer,1,2 Ana Quirce,3 Angel Valle,1,* Luis Pesquera,1 Miguel A. Rodríguez,1 Krassimir Panajotov,3,4 and Marc Sciamanna2 1Instituto de Física de Cantabria, Consejo Superior de Investigaciones Científicas (CSIC), Universidad de Cantabria, 39005 Santander, Spain 2Chair in Photonics, LMOPS Laboratory, CentraleSupélec, Université de Paris-Saclay and Université de Lorraine, 57070 Metz, France 3Vrije Universiteit Brussel, Faculty of Engineering Sciences, Brussels Photonics Team (B-PHOT), Pleinlaan 2, 1050 Brussels, Belgium 4Institute of Solid State Physics, 72 Tzarigradsko, Chaussee Blvd., 1784 Sofia, Bulgaria


I. INTRODUCTION
Multistable systems are characterized by the coexistence of different attractors for a fixed set of parameters.These systems are extremely sensitive to noise because even small strength fluctuations may switch the system from one attractor to another [1].Multistable systems can have complex dynamics in which there is coexistence of competing behaviors on which the trajectory alternates [2][3][4].An example of these dynamics is the attractor hopping [5][6][7][8] in which, as the system tries to follow a regular motion in the neighborhood of one attractor, then noise causes jumping among the different states.An example of attractor hopping is the two-state on-off intermittency that can be found in bistable systems [9].In the on-off intermittency there is a laminar phase in which a variable of the system takes small values for times that can be very long [10].This regular state is interrupted by irregular bursts (turbulent phase) [10].On-off intermittency is characterized by a −3/2 power law for the probability distribution of the laminar phases versus the laminar length [11].This scaling relation has been found in experiments with electronic circuits [12], nematic liquid crystals [13], distributed-feedback semiconductor lasers (DFB) [14], diode lasers with external cavity [15], optically injected dual-mode semiconductor lasers [16], and mutually coupled diode lasers [17].Noise-induced attractor hopping has been also experimentally observed in multistable systems like fiber lasers subject to externally applied noise in the diode * valle@ifca.unican.espump current [8].A multistate intermittent regime in which the trajectory visits two or more periodic states has been observed [8].For the case of hopping between two states the probability distribution of periodic windows was characterized by the −3/2 power law [8] as in the case of two-state on-off intermittency.Attractor hopping phenomena have been also observed in CO 2 lasers [18,19].Quasi-isotropic CO 2 lasers can present polarization bistability with spontaneous flips between the two orthogonal linear polarization eigenmodes [18].Transverse multimode CO 2 lasers with optical feedback also have bistability and alternation of different polarized patterns [19].
Many of the previous experimental studies have been performed using semiconductor lasers.Dynamics of diode lasers is quite susceptible to external influences, such as optical feedback or injection [20].Nonlinear dynamics of semiconductor lasers subject to optical injection have been investigated for different types of devices such as edge emitters (DFB, discrete mode lasers, multisection) [21] and verticalcavity surface-emitting lasers (VCSELs) [22][23][24][25].VCSELs are intrinsically single-longitudinal mode devices, and when emitting in the fundamental transverse mode, two orthogonally polarized modes can be observed [26].Polarization switching (PS) between these two orthogonal modes can be achieved with temperature or bias current changes [27][28][29][30][31] but also through optical injection [32][33][34].PS in VCSELs (see [35], and references therein) has been widely characterized with regard to its bistable properties when subject to optical injection [36][37][38], and possible applications in all-optical signal processing and data transmission have been demonstrated [38][39][40].The PS process is intrinsically connected to noise, which plays an important role in laser dynamics.Indeed, spontaneous emission noise can strongly modify laser dynamics [41], induces random spikes and intermittent switches in bistable regions [42], or causes switchings between bistable linear polarizations in VCSELs following Kramers hopping behavior [43].Noise can also play a positive role inducing coherence resonance in edge-emitter semiconductor lasers subject to optical injection [42,44], in the polarization mode hopping regime of a VCSEL subject to optical feedback [45], or stochastic resonance in the polarized emission of a pump modulated VCSEL [46].
Recently, the polarization nonlinear dynamics of a singlemode VCSEL subject to parallel optical injection has been studied in experimental and theoretical ways [47][48][49][50].In the parallel optical injection the directions of linear polarization of the injected light and the free-running VCSEL are parallel.Rich nonlinear dynamics have been observed involving the polarization of the solitary VCSEL (parallel) and the orthogonal polarization (fixed point, limit cycle, period doubling, chaos) [50].A novel state in which injection locking of the parallel polarization and excitation of the free-running orthogonal polarization occur simultaneously (IL+PS) was found [47].Bistability between the IL+PS and the single polarization mode injection locked (IL) solutions was described [49].Also random switchings between some of the previous dynamical states were experimentally observed [50].In particular, we found switchings between (i) a period-2 in the parallel polarization (P2) and a period-1 dynamic in both polarizations (P1both) and (ii) a period-1 in the parallel polarization (P1) and IL+PS.These switchings are examples of attractor hopping between polarization dynamical states induced by noise.
In this work we characterize these switchings and shed some light on the origins of the noise or noises that induce this attractor hopping.We first study the hopping dynamics between P1 and IL+PS by analyzing long experimental time series.Close to the boundary of the region in which this behavior appears we find that the residence time in one of those states can take values spanning several orders of magnitude in such a way that its probability distribution is characterized by a −1.83 ± 0.17 power law.We also identify coherence enhancement of the attractor hopping in which transitions between P1 and IL+PS occur regularly.We report another example of attractor hopping between two polarized dynamical states: chaos in the parallel polarization (CH) and chaos in both polarizations (CH-both).We also experimentally observe two examples of attractor hopping between more than two states: at negative frequency detuning we observe switching between CH-both, IL, and CH while at positive detuning we find switching between P1, IL+PS, and P1-both.Good agreement between our experiments and the numerical results of a rate equations model allows us to identify the origin of the noise that triggers these switchings.
There are several characteristics that differentiate our results from previous works on on-off intermittency.In the typical on-off intermittency a Hopf bifurcation is crossed in such a way that there is an intermittent regime in the laminar phase associated to the fixed point (off) while no intermittent regime appears in the phase corresponding to the limit cycle [10,11].This is an asymmetric behavior due to the asymmetry of the Hopf bifurcation.In contrast, our experiment reports a symmetric behavior in which intermittent regimes can be observed in the laminar state associated to our fixed point (IL+PS, off) or in the limit cycle (P1, on).We also obtain heavytailed residence time distributions for both previous cases with an intermediate region with coherence enhancement.These heavy-tailed residence time distributions are characterized by power laws that span a much larger time range (up to three decades) than those observed in previous experiments with lasers [8,16,17].In our experiment the laser can emit in two different polarizations.The dimensionality of our system is larger than that observed in previous works with on-off intermittency in lasers in which only a single linear polarization was observed.Our experimental residence time distributions are characterized by a −1.83 ± 0.17 exponent, instead of the −3/2 observed in typical on-off intermittency.This difference in exponents and the previous symmetric behavior are possibly due to the fact that the Hopf bifurcation occurs in a space with a higher dimensionality.
The article is organized as follows.First we describe in Sec.II the experimental setup used to attain parallel optical injection, then we present in Sec.III the fundamental states that are involved in the different multistate hopping dynamics that appear in this work.In Sec.IV we focus on the attractor hopping between P1, IL+PS, and P1-both.Section V is devoted to switching between CH-both, IL, and CH.In Sec.VI we present our theoretical results and finally in Sec.VII we present some discussion and conclusions.

II. EXPERIMENTAL SETUP
Experimentally we attain parallel optical injection in the VCSEL by using the all-fiber setup illustrated in Fig. 1(a) and described in detail in [48][49][50].More specifically, we inject the parallely polarized light from our tunable master laser into the slave laser, which is a long-wavelength 1550-nm VCSEL, through an optical circulator.This slave laser is a commercially available quantum-well VCSEL (RayCan TM ), based on the InAlGaAs-active region.
The threshold current of the VCSEL at 298 K is I th = 1.66 mA, and the temperature is held constant to this value.Measurements presented in Sec.IV (Sec.V) are performed at a value of the bias current I bias = 5.0 mA (I bias = 3.0 mA). Figure 1 shows the experimental optical spectra, obtained with a high-resolution (10 MHz) optical spectrum analyzer (BOSA), for (b) I bias = 5.0 mA and (c) I bias = 3.0 mA.At these bias currents, the device emits in a single transverse mode, and the two modes appearing in the spectrum are the two orthogonal polarizations of this single transverse mode.The lasing mode, called parallel (x), emits for I bias = 5.0 mA (I bias = 3.0 mA) at λ = 1541.82nm (λ = 1540.91nm), chosen to correspond to the zero frequency.The subsidiary mode, that is, the orthogonal polarization (y), is shifted by 33.36 GHz toward the high frequencies, corresponding to λ ⊥ = 1541.56nm (λ ⊥ = 1540.65nm) with a side mode suppression ratio of 36 dB (34 dB).We note that these values of suppression ratios do not correspond to the spectra shown in Fig. 1 but to spectra measured with another optical spectrum analyzer with a better power resolution although with a lower frequency resolution (8 GHz).Several PS can be found in the free-running VCSEL for corresponding bias currents: 2.25, The optical injection is characterized by its strength P i measured by a power meter after a 50:50 optical coupler which other branch is connected to the slave VCSEL, and by its frequency detuning ν i , defined as the difference between the frequency of the injected light and the frequency of the free-running lasing mode.The value of ν i is determined at the minimum value of P i at which the spectral peak corresponding to the optical injection becomes visible at BOSA.In this way, optical injection affects in a minimum way the value of ν i because the frequency of the parallel polarization is barely shifted.
Time series are obtained with a real-time high-speed (12 GHz) oscilloscope with two high-speed photodetectors (9 GHz bandwidth).We have measured a 0.4 ns electrical delay between the two time traces in the oscilloscope.Time series corresponding to the orthogonal polarization are then shifted 0.4 ns in order to have simultaneous time traces.They have also been scaled in order to compensate the polarization-dependent losses induced by the polarization beam splitter.The parallel polarization signal is shifted upward (0.057 a.u) in order to avoid overlapping of the two signals, for clarity purposes.
We present for I bias = 5.0 mA a period-1 behavior in the parallel polarization (P1) in Figs.2(a) and 2(b), period-1 appearing simultaneously in both polarizations (P1-both) in Figs.2(c) and 2(d), and the state where both injection locking of the parallel polarization mode and excitation of the orthogonal polarization occur simultaneously (IL+PS) in Figs.2(e) and 2(f).These three states have been obtained consecutively for a fixed value of the detuning (ν i = 0.9 GHz) and by decreasing the injection level at, respectively, P i = 303.8μW, P i = 144.8μW, and P i = 31.62μW.P1 is characterized by a limit cycle in the parallel polarization at 3.3 GHz, which is close to the relaxation oscillation frequency at this bias current, and centered on the injected light frequency ν i = 0.9 GHz.P1-both is characterized by an anticorrelated limit cycle appearing simultaneously in both polarizations with frequencies of 3.3 and 3.4 GHz, and centered on the injected light frequency and the free-running orthogonal mode, respectively.Finally, IL+PS, which has been widely discussed in our previous work [47,48], shows the locking of the parallel polarization to the injected light, and a simultaneous excitation of the orthogonal polarization mode.These three types of nonlinear dynamics will appear in the attractor hopping discussed in Sec.IV.
We present for I bias = 3.0 mA a chaotic behavior in the parallel polarization (CH) in Figs.2(g) and 2(h), and a chaotic behavior appearing simultaneously in both the parallel and the orthogonal polarization modes in Figs.2(i) and 2(j).These two chaotic behaviors have been obtained consecutively for a fixed value of the detuning ν i = −2.2GHz, by decreasing the level of injection at, respectively, P i = 19.87μW (CH) and P i = 14.21 μW (CH-both) [50].Chaotic behaviors are characterized by broadened spectra and aperiodic time series, respectively, in the parallel polarization mode for CH, and simultaneously in both polarization modes for CH-both.These two dynamics, along with injection locking of the parallel polarization, will take part in the hopping dynamics that will be discussed in Sec.V.

IV. ATTRACTOR HOPPING: P1-P1-BOTH-IL+PS A. Experimental observations
We now focus on the hopping dynamics found between the P1 state and the IL+PS state.Figure 3 shows the experimental optical spectrum and the corresponding time traces of this dynamics, named (P1 & IL+PS).The experimental optical spectrum of Fig. 3(a) shows a period-1 behavior in the parallel polarization mode at 3.0 GHz, centered on the injected light frequency, and a constant excitation of the orthogonal polarization.This spectrum cannot be explained without looking at the time trace of Fig. 3(b).From Fig. 3(b) we conclude that the spectral behavior results from a random switching between the P1 state and the IL+PS state.The optical spectrum is obtained by using a temporal window in which many of these switching events occur.So this spectrum is a superposition of the spectra corresponding to the P1 and IL+PS states.These switchings do not show periodicity, hence a complete characterization of this hopping statistics with regard to residence time on each state will be performed in the next section.P1 & IL+PS is a good example of a dynamics that needs the time series of the power of both linear polarizations in order to understand the structure of its optical spectrum.Characterization of the dynamics by optical spectra was used to obtain the experimental stability maps presented in [50].We present in Fig. 4(a) a map-in the injected power-frequency detuning plane-in which we report a variety of two-state or multistate hopping dynamics, for which spectra are not different from the spectrum of the fundamental state P1-both presented in Fig. 2(c).This map was obtained by decreasing the injection level from a high value for a fixed value of the detuning, and repeating the process for several values of the detuning.We only focus on the part of the P i -ν i plane where hopping dynamics has been observed.At these levels of injection for lower detuning, P1 or P1-both dynamics can be found, and for higher values of the detuning, the IL+PS state appears [50].Maps included in [50] cover a much wider region of the P i -ν i plane.However, they give less information in the hopping region because different dynamics identified with the oscilloscope have similar optical spectrum.In fact, only P1-both, P1 & IL+PS, and IL dynamics are identified in [50] for the region considered in Fig. 4(a).
Figure 4(a) shows three regions where hopping dynamics appear.First, P1 & IL+PS dynamics results from a switching between P1 and IL+PS, for which optical spectra and time series have been studied in Fig. 3. Then P1 & P1-both dynamics results from a switching between P1 and P1-both states, for which time traces are given in Fig. 4(b).These time traces show an aperiodicity in the hopping between these two states.In this example, the residence time in the P1-both state can vary quite a lot, while the residence time on the P1 state is almost constant.This aperiodicity is similar to that found for the case of the P1 & IL+PS dynamics.And finally we find P1 & P1-both & IL+PS dynamics, in which the system jumps between three states: P1, P1-both, and IL+PS.The time trace of this state is given in Fig. 4(c).Figure 4(d) shows the orthogonal polarization mode power as a function of the parallel polarization mode power, giving a representation of the three attractors.In this representation, the P1 state is when the power of the orthogonal polarization is close to zero and the parallel signal oscillates, the IL+PS state is when both polarizations are almost constant, and the P1-both state is the closed curve, showing that both powers oscillate and are anticorrelated.We do not show experimental optical spectra for either P1 & P1-both or P1 & P1-both & IL+PS since the spectra would be very similar to the P1-both spectrum.
In the P1 & P1-both & IL+PS dynamics where the system jumps between three attractors, the analysis of time traces shows that the hopping always follows the same sequence, which is P1, then IL+PS, and finally P1-both.This sequence is seen in the time traces of Fig. 4(c) and is illustrated by arrows in Fig. 4(d).The origin of this fixed sequence of hopping dynamics will be discussed in Sec.VII.

B. Residence time distributions
We now characterize the hopping statistics focusing on the switching between the two attractors P1 and IL+PS, which is the P1 & IL+PS dynamics described in Fig. 3. Figure 5  .The time series only show the orthogonal polarization signal for clarity purposes.Thus the system is in the P1 state when the signal is almost at 0 (a.u), and the system is in the IL+PS state when the signal is at its higher values, around 0.012 (a.u) as shown in Fig. 3(b).We define τ 1 and τ 2 as the residence time in the P1 and IL+PS state, respectively.These residence times are obtained by detecting when the time trace of the orthogonal power crosses a fixed level placed at the middle of the variation range, as illustrated in Fig. 5(e).In order to obtain enough statistics, we have saved 2.0-s-long time traces with a 0.1 μs sampling time.This sampling time is small enough to detect with good precision the values of both residence times.
We focus at first on the left panels, Figs.5(a), 5(d) and 5(g).The time traces, and the corresponding residence time distribution were obtained at ν i = 0.8 GHz and P i = 48.01μW, which is close to the region where P1 is observed [50].The time trace, Fig. 5(a), and its zoomed part, Fig. 5(d), show that the system mainly stays on the P1 state, and randomly jumps for a few tens of microseconds to the IL+PS state.We plot in a log-log scale the residence time distribution on each state in Fig. 5(g), obtained from one 2.0-s-long time trace.We see that the residence time on the IL+PS state, τ 2 in red, can vary from 10 to 100 μs, and that τ 1 in black can vary from 8 μs to 10 ms, which is about three orders of magnitude.
The centered panel, which is obtained for injection parameters ν i = 0.9 GHz and P i = 53.83μW, shows time series [Figs.5(b) and 5(e)] where the switching between the P1 and the IL+PS states is more predictable and seems almost periodic.The residence time distribution, Fig. 5(h), for a 2.0 s time trace confirms that the stay on each state is almost constant around 8.0 μs.This behavior appears in the major part of the P1 & IL+PS region of the injected power-frequency plane, while the situations described by left and right panels are obtained close to the regions where the fundamental major state can be found, which is the P1 or the IL+PS attractor, respectively.
A very long tail distribution can also be found in the opposite case in the right panel where the system mainly stays on the IL+PS state and jumps to the P1 state, as shown in Figs.5(c) and 5(f).The residence time distribution for one 2.0-s-long time trace is given in Fig. 5(i) and shows that the system stays for some microseconds on the P1 state (τ 1 ), whereas the stay on the IL+PS state, given by τ 2 , can vary from 5 μs to 10 ms.This panel was obtained for the injection parameters ν i = 1.1 GHz and P i = 64.71μW, which is close to the region where IL+PS is observed [47,48,50].
In order to characterize the decay of one of these heavytailed distributions, we have chosen the case illustrated in Fig. 5(i).We have obtained 24 consecutive 2.0-s-long time traces for the same conditions of Fig. 5(i).The corresponding 24 residence time distributions are plotted with dashed colored lines in Fig. 6.Most of them follow the same trend shown in Fig. 5(i).Figure 6 shows with open dots the residence time distribution obtained with all the data (around half a million).A linear fit of the distribution over three orders of magnitude (from τ 2 = 10 −5 s to 10 −2 s) is also shown with a thick dashed line.This fit shows that the probability distribution is characterized by a −1.85 power law.In order to quantify the uncertainty of this number, we calculate the average and standard deviation of the 24 power laws obtained with our 24 time traces.We obtain that the probability distribution is characterized by a −1.83 ± 0.17 power law.Straight lines corresponding to −3/2 and −2 power-law decays are also included in Fig. 6 for the purpose of comparison.
We now use interspike-interval fluctuations to quantitatively characterize the coherent behavior shown in Fig. 5(h).We can define interspike intervals (ISIs) between consecutive optical switches as T i = t i+1 − t i , where t i is the time when a spike occurs [42].In terms of our consecutive residence times τ 1,i and τ 2,i we can also write T i = τ 1,i + τ 2,i .ISI fluctuations are defined as the difference between consecutive natural logarithms, i = ln T i+1 − ln T i .The standard deviation of i for the cases illustrated in Figs.5(a), 5(b) and 5(c) are 1.4,0.05, and 1.16, respectively.In this way the system dynamics exhibits maximum coherence for the case illustrated in Fig. 5(b).In coherence enhancement the regularity is optimized for a certain value of a parameter of the system, while the noise intensity is constant [42].That situation has also been reported in the literature as deterministic coherence resonance [51,52].This is sometimes distinguished from coherence resonance in which maximum regularity is obtained for a certain value of the noise strength [42].In our case we show coherence enhancement of the attractor hopping as we do not have control over the strength of noises present in our experimental system.

V. ATTRACTOR HOPPING: CH-CH-BOTH-IL
In this section we focus on some situations found when decreasing the VCSEL current to I bias = 3.0 mA.Complete experimental maps corresponding to this bias current can be found in [50].We present in Fig. 7(a) experimental time traces obtained for ν i = −2.2GHz, and when decreasing the injected power to P i = 13.56 μW.This figure shows attractor hopping between a chaotic behavior in the parallel polarization and a chaotic behavior appearing simultaneously in both polarization modes.Again the residence times in each of these dynamical states are random in the microsecond time scale.Similarly to the attractor hopping involving the P1-both state [Figs.4(b) and 4(c)] for which the optical spectra were the same as the optical spectrum of the P1-both state, the optical spectra of attractor hopping dynamics involving the CH-both and CH states are similar to the optical spectrum of the CH-both behavior of Fig. 2(i), hence we do not give these spectra.
An example of hopping between three different states is given in Fig. 7(b) that has been found in the same conditions as Fig. 7(a) but for a slightly larger P i = 15.08 μW.This figure illustrates switching between CH, CH-both, and injection locking of the parallel polarization state.The system jumps from a CH-both state to IL, and then to CH.Once in CH we observe switching between this state and IL, until the system finally jumps back from CH to CH-both.Repetition of this cycle is also observed in Fig. 7(b).
Length of time series in Fig. 7 is limited to 500 μs because the sampling time is much smaller than in Fig. 5 in order to check that the signal has the fast chaotic oscillations that characterize CH and CH-both.This limits the amount of residence time data that we can get and so experimental residence time distributions have not been obtained, in contrast to the P1 & IL+PS case.

VI. THEORETICAL ANALYSIS
The attractor hopping described in previous sections is due to noise.The optical injection system is inherently noisy due to spontaneous emission in both lasers, noise in the current source, and the temperature controller.These effects produce random variations of power, wavelength, and phase for both lasers.One of the aims of this section is to investigate the type of noise responsible for the observed switchings.Another objective is to distinguish if multistability is important for the explanation of our results.Switchings can occur because noise induces hoppings between different multistable states.Also, hoppings can be due to changes of the parameters of the system (deterministic or random) that cause crossings between the frontiers of the different dynamical regions, as in [24].
The tool that we will use in order to shed some light on this question is the theoretical analysis using a physical model of our system.We use the spin-flip model [28,49], which is a widely used rate equation model describing the polarization modes of a single-mode VCSEL, in which we have added an optical injection term.The model equations are given in (1)- (5).In these equations, E x and E y are the two linearly polarized slowly varying components of the field in the x and y direction, corresponding to parallel and orthogonal polarizations, respectively.D and n are two carrier variables.D = (N − N t )/(N th − N t ) where N, N th , and N t are the carrier number, and the carrier number at threshold and transparency, respectively [53].n corresponds to the difference between the population inversions for the spin-up and spindown radiation channels.μ is the normalized bias current, given in terms of I and I th , which are the bias current and the threshold current, respectively.
The meaning of the remaining parameters is the following: κ is the field decay rate, α is the linewidth enhancement factor, γ is the decay rate of D, τ n is the differential carrier lifetime at threshold, τ e is the carrier lifetime at threshold, γ a is the linear dichroism, γ p is the birefringence parameter, and γ s is the spin-flip rate.Spontaneous emission noise is included by using Gaussian white noises, ξ ± .The statistical properties of these noises and the expressions of spontaneous emission noise rates, R ± , are given in [49].E inj is the amplitude of the field of the injected light and ν inj is the detuning between the injected light and the intermediate frequency between those of the x and y polarization, ν x and ν y , where 2πν x = αγ a − γ p and 2πν y = γ p − αγ a .Thus the frequency detuning used to characterize the injected light in our experiment is ν i = ν inj − ν x .We numerically integrate the model equations using the numerical values of the parameters given in [49].
The model parameters that we use were extracted for a similar VCSEL device [53,54].Using these parameters, the agreement between experimental and theoretical results describing some of the dynamics in our system was very good (see Fig. 7 of [49]).We note that our device is the same as that used in [49], although only similar to the other device characterized in [53,54].

A. Attractor hopping between CH and CH-both
We first analyze the hopping between CH and CH-both by using numerical simulations of the model.We integrate Eqs. ( 1)-( 4) for the values of I bias = 3.0 mA and ν i = −2.2GHz corresponding to Fig. 7(a).In our simulations we take the value of the injected power P i = E 2 inj = 0.024 (a.u.) that corresponds to the injected power of Fig. 7(a) after applying the conversion factor between the theoretical and experimental injected power in our setup [1 (a.u.) ↔ 565.76 μW] [49].We also take the value of the linear dichroism, γ a = −0.21ns −1 , that corresponds to I bias = 3.0 mA [48], and a fraction of spontaneously emitted photons that are coupled into the laser mode, β SF = 2.7 × 10 −4 .
Figure 8(a) shows the theoretical time traces of the power of both linear polarizations.This figure has been obtained after a decrease in the injected power, similarly to the experiment.The numerical simulation process is as follows.Equations ( 1)-( 4) are integrated by using a continuous decrease of the injected power from a large value to the desired value P i = 0.024 (a.u.), at which the time series is recorded for the desired length, with a 0.01 ps integration time, and a 5.0 ps sampling time.Black (red) traces correspond to the parallel (orthogonal) polarization signal, and the parallel signal has been shifted upward in order to avoid overlapping, for clarity purposes.proportional to the fraction of spontaneously emitted photons that are coupled to the lasing mode, β SF [49].We note that the numerical value of β SF in Fig. 8(a), 2.7 × 10 −4 , for which transitions between CH and CH-both occur in similar time scales as in the experiment, is in the same order of magnitude as the value measured in [53] for a similar device.We also note that the simulation of the almost deterministic version of our model (β SF = 10 −30 ) gives no transitions between CH FIG. 9. Simulated residence time distributions in log-log scale obtained for I bias = 3.0 mA, ν i = −2.2GHz, P i = 0.024 (a.u), and a noise strength of β SF = 6.5 × 10 −4 .Black (red, dotted) trace is the parallel (orthogonal) polarization mode signal.and CH-both.This is shown in Figs.8(e) and 8(f) in which results corresponding to two different simulations of the almost deterministic model for two different sets of initial conditions have been plotted.Therefore the spontaneous emission noise induces transitions between CH and CH-both.These transitions are observed for parameter values in which both solutions, CH and CH-both, are stable, as no model parameters are changed during the simulation Fig. 8(a).this way, this is bistable system in which transitions are induced by emission noise in VCSEL.The amount of residence time data that we can numerically get for β SF = 2.7 × 10 −4 is small due to long computational times.A comparison between Figs. 8(a) and 8(b) shows that with β SF = 6.5 × 10 −4 we obtain a larger number of residence time data that we can use to calculate the residence time distributions.These distributions are shown in Fig. 9, in which τ 1 and τ 2 are the residence times in CH-both and CH, respectively.Both distributions are plotted in a log-log scale showing a power-law decay, although over a smaller range of time than in Fig. 5 due to the small number of numerical data.Both distributions are characterized by a −1.4 power law.

B. Attractor hopping: P1-P1-both-IL+PS
Attractor hopping between P1, P1-both, and IL+PS states happens at a different I bias = 5.0 mA, so in our simulations we take into account this value and its corresponding linear dichroism γ a = −0.36ns −1 [48].In order to see if multistability is important for the explanation of our results, we first calculate bifurcation diagrams in which the control parameter is ν i .In this way we will try to identify possible multistability regions.The bifurcation diagram is calculated by fixing a value of the injected power and increasing linearly ν i , step by step.The duration of each step, in which ν i is constant, is 201 ns.In the last nanosecond of each step all the maxima and minima of the power of both linear polarizations are plotted in Fig. 10.The results corresponding to a repetition of this process but now decreasing ν i back to its initial value are also shown in Fig. 10.Note that this situation does not correspond to that illustrated in FIG.10.Simulated bifurcation diagram showing the maximum and minimum of the power of both linear polarizations when increasing (squares), and decreasing (open dots) ν i .Results are obtained for I bias = 5.0 mA, P i = 0.022 (a.u), and β SF = 10 −30 .Black (red, lower) symbols correspond to the parallel (orthogonal) polarization mode.
Fig. 4(a) because frequency sweeping is performed for a fixed injected power.In order to have almost deterministic results, we have made these calculations with a value of β SF = 10 −30 .Three different solutions appear in Fig. 10: P1, P1-both, and IL+PS.P1 is characterized by oscillations of the parallel polarization power at 3.25 GHz, very close to the theoretical value of the relaxation oscillation frequency (3.29 GHz) [48].The power of both linear polarizations oscillates at 3.25 GHz when the state is P1-both.In Fig. 10, for increasing ν i , P1, P1both, and IL+PS appear for 0.75 ν i 0.774, 0.774 < ν i 0.784, and ν i > 0.784 GHz, respectively.Figure 10 also shows the hysteresis region in which we find bistability between P1 and P1-both for 0.768 ν i 0.774 GHz.
We now fix the value of ν i in the bistable region in order to see if spontaneous emission noise is able to induce switchings between stable states.We take ν i = 0.774 GHz and we choose the value of β SF = 2.7 × 10 −4 with which we have found results similar to the experimental ones in the previous section.Figure 11 shows the time traces corresponding to the power of both linear polarizations and the total population inversion, D. This figure illustrates that there are switchings between the different states.Figure 11(a) shows that the system can be in the P1 state; for instance, between 1146 and 1150 ns the power of the parallel polarization oscillates from 2.2 to 4.2, while the power in the orthogonal polarization is negligible (see also Fig. 10).Also, the system can oscillate in the P1-both state; for instance, between 1238 and 1243 ns the amplitude of the oscillation of the power of the parallel polarization (between 2.4 and 3.7) is smaller than for P1, while the orthogonal polarization oscillates with small amplitude (see also Fig. 10).We also observe some switchings to the IL+PS state (see, for instance, from 1217 to 1226 ns): oscillations in the parallel polarization power shrink close to the analytical value, 2.78, calculated in [48].Also, the power of the orthogonal polarization and the value of D are close to those calculated with the theory of Ref. [48] (0.29 and 1.011, respectively).However, although excitation of P1, FIG.11.Simulated time series of (a) of both linear polarizations, (b) the total population inversion.Results are obtained for I bias = 5.0 mA, ν i = 0.774 GHz, P i = 0.022 (a.u), and β SF = 2.7 × 10 −4 .Black (red, lower) lines correspond to the parallel (orthogonal) polarization mode.P1-both, and IL+PS is possible in this region, we note that the temporal dependence is quite different from that reported in our experimental results of Fig. 3(b), Fig. 4(b), or Fig. 4(c).First, experimental switchings occur in the microsecond time scale, while the scale of theoretical switchings is in the nanosecond range.Second, in the case of excitation of the three states there is a definite experimental order in which states appear, while in theory the order is random and so experimentally observed order, P1, IL+PS, P1-both, is not observed.These results suggest that only the effect of spontaneous emission noise is not enough to explain the observed dynamics.
We now explore the possibility in which hopping is mainly due to changes of the parameters of the system (deterministic or random) that cause crossings between the frontiers of the different dynamical regions, as in [24].This is suggested by the experimental time dependence observed in the power of both linear polarizations.For instance, Fig. 3(b) shows that just after switching from P1, the power of each of the polarizations is not constant, in contrast to the theoretical results of [48].The same happens just after switching from IL+PS: the amplitude of the power oscillations increases in time.This increase saturates in a microsecond time scale [see Fig. 4(b)].Also, the power of each linear polarization in the IL+PS state reaches a constant value in a similar time scale.The fact that changes in the amplitude of solutions occur in the microsecond time scale suggests that thermal effects can play an important role in this behavior.
A simple way to include temperature changes in our model is by considering a time dependence of ν i .The VC-SEL resonance frequency, ν || , is inversely proportional to the refractive index inside the device.An increase of the temperature causes an increase of the refractive index and so ν i = ν inj − ν || increases.Noise in the current and temperature controllers cause changes of the temperature of the device in the microsecond thermal time scale.Changes in ν i are smaller than 0.05 GHz (0.15 GHz) for changes in current (temperature) smaller than the current (temperature) controller resolution.As a first approach to include this effect we consider a periodic change of ν i between two values for which only one of the states is stable.We have chosen ν i = 0.75 and 0.79 GHz for which only P1 and IL+PS are stable, respectively (see Fig. 10).We choose a periodic time dependence of ν i between two values, ν a and ν b , characterized by an exponential approach with a thermal time constant, τ th , and a period P , given by This dependence is illustrated in Fig. 12(a) in which ν a = 0.75 GHz, ν b = 0.79 GHz, τ th = 1 μs, and P = 10 μs.Initially, at 120 μs, ν i is close to 0.75 GHz and P1 state is observed (see also Fig. 10).In Fig. 12 the power of the parallel linear polarization has been shifted upward in four units for a better comparison with our experiments.When ν i begins to increase at 120 μs there is a transient in which P1 with an increasing amplitude appears.After this state P1-both is briefly excited between 120.74 and 120.76 μs.These values correspond to ν i ∼ 0.771 GHz, which is slightly below the value at which transition to P1-both appears in Fig. 10.After P1-both there is a transition to IL+PS.The state appearing after 120.76 μs is an IL+PS state because its total population inversion is very close to the constant value, 1.011, that characterizes this state [48].This is a transient IL+PS state because both linearly polarized powers change in time in such a way that the orthogonal polarization is dominant.Experimental observation of this IL+PS transient state is clear from Fig. 3(b).The transient IL+PS state is opposite to the IL+PS steady state in which, for this range of ν i the parallel polarization is dominant because P x = 2.60 and P y = 0.47 for ν i = 0.79 GHz [48].The observed decrease (increase) of the orthogonal (parallel) polarization power during IL+PS is the result of the transition to the steady state corresponding to ν i = 0.79 GHz, which is almost reached at 125 μs.When ν i begins to decrease at 125 μs the system goes to a transient P1 state after which IL appears at 125.9 μs.This state is IL because D goes well below 1, a characteristic situation of single-mode injection locking.This solution is not stable in this ν i range and the system goes to the stable solution, which is the P1 state.In this way the ν i cycle closes and dynamics is repeated.We note that two P1 oscillations with two different amplitudes when leaving the IL+PS state are also observed in our experimental [see Fig. 3(b) at t ∼ 15 μs and more in Fig. 4(f) of Ref. [50]].The sequence of dynamical states discussed above also appears when changing the frequency and amplitude of the ν i cycle and the value of τ th .Also, the same sequence appears at smaller values of time at which the permanent regime has not been reached yet.The reason for the delayed bifurcations observed in Fig. 12 is the critical slowing down that appears in dynamical systems with time-dependent parameters [55].
We have obtained numerical residence time distributions obeying a power law for the case of CH and CH-both hopping (see Fig. 9).These numerical results have been obtained with a rate equation model in which the only noise that is included is the spontaneous emission noise.In this case, the spontaneous emission noise is enough to cause hopping between CH and CH-both states.Taking into account these results, it seems plausible that the inclusion of random changes of parameter values (frequency detuning) in time would also result in the observation of time distributions obeying a power law, similar to that shown in Fig. 6.However, this simulation is beyond the scope of this work because intensive calculations would be required in order to get numerical residence time data similar to the experimental ones.We note that the experimental data can be of the order of several milliseconds, while the numerical integration time step, 0.01 ps, is almost 12 orders of magnitude smaller.In this way obtaining enough data to get a residence time distribution directly comparable to the experimental one would require computational resources beyond our reach.

VII. DISCUSSION AND CONCLUSIONS
We have analyzed in more depth two cases of attractor hopping.First, we have shown that CH and CH-both hopping is induced by spontaneous emission noise that causes switchings between these two stable states.In the second case, hopping between P1, P1-both, and IL+PS, spontaneous emission fluctuations in a bistable region do not seem enough to explain the observed behavior.Changes in some system parameter, such as the frequency detuning, seem necessary to explain the behavior.The fixed sequence of switchings between IL+PS, P1-both, and P1 shown in Fig. 4(c) can be explained by looking at Figs. 2(c) and 2(d) of [50].In these figures the same sequence is observed at ν i = 0.9 GHz increasing and decreasing P i around 100 μW.In this way, bistability between these solutions can explain the fixed sequence observed in Fig. 4(c), which is increasing and decreasing ν i for a fixed value of P i .Also, changes in the frequency detuning permits us to explain some characteristics of the nonlinear dynamics of our system.However, we have considered a periodic dependence of ν i for which most of the transitions between states happen in times determined by that period.This occurs because these transitions used to occur when there are crossings between regions in which only one of these solutions is stable.A good description of the statistics of the time the system spends in each of these states requires not only the consideration of spontaneous emission noise but also of realistic fluctuations in some system parameters, such as the frequency detuning, possibly induced by fluctuations in the current and temperature controllers.
Both fluctuations and excursions of the frequency detuning can help to interpret the results shown in Fig. 5. Excursions of the frequency detuning around its mean value can have different amplitudes.For small amplitudes the system is in a bistable region and switchings are due to spontaneous emission noise.For large amplitudes switchings are due to the loss of stability of a solution.Both switching mechanisms can play a relevant role in our system.The detailed temporal dependence of the frequency detuning determines which is the dominant mechanism.Heavy-tailed residence time distributions of Figs.5(g) and 5(i) appear when the frontier of the region in which only one of the solutions is stable is approached.Coherence enhancement of hopping between IL+PS and P1 states has also been described in Fig. 5(e).The power of the orthogonal polarization has a period close to 16 μs, which corresponds to a frequency around 60 kHz.The origin of this coherence enhancement could be linked to coherence or stochastic resonance phenomena since the bistable behavior, observed in our system, is essential for both.We think that coherence resonance effects seem less plausible since we are not aware of any internal process in the VCSEL characterized by a frequency close to 60 kHz.In contrast, stochastic resonance could be the reason for the observed coherence enhancement because current laser controllers are characterized by current noise density spectra with narrow and strong peaks in the 10-100 kHz range [56].One of these oscillations could be the weak periodic signal 60 kHz necessary to stochastic resonance.The origin of oscillations at 60 kHz could also be attributed to resonances in the current noise density spectra of the temperature controller or to the effect of electromagnetic radiation noise on the connection cable to the laser.
Summarizing, we have reported noise-induced attractor hopping between polarization dynamical states in a single transverse mode VCSEL subject to parallel optical injection.We have experimentally characterized several types of attractor hopping.We have obtained an experimental map identifying regions where attractor hopping between two or more states occurs.We have found multistability regions that are characterized by heavy-tailed residence time distributions in which residence times can take values spanning more than three orders of magnitude.These distributions are characterized by a −1.83 ± 0.17 power law.Between these regions we have obtained coherence enhancement of attractor hopping in which transitions between states occur regularly.Our theoretical analysis using a rate equation model describing the VCSEL polarization has shed some light on the origin of the noises that trigger these hoppings.Simulations of this model have shown that frequency detuning variations and spontaneous emission noise play a role in causing switching between attractors.This is shown for hopping between P1, P1-both, and IL+PS states.For the case of P1 and IL+PS hopping, it has not been possible to obtain numerically the residence time distribution obeying a power law because of difficulty in numerical simulations.However, in the case of CH and CH-both hopping the effect of spontaneous emission noise alone is enough to explain the observed dynamics.In this case, we have obtained numerical residence time distributions that obey a power law.

FIG. 3 .
FIG. 3. Experimental optical spectrum (a) and corresponding time trace (b) of the (P1 & IL+PS) state, resulting from an aperiodic switching between P1 and IL+PS.Black (red, lower) trace corresponds to the parallel (orthogonal) polarization signal.The injection parameters are I bias = 5.0 mA, ν i = 0.9 GHz, and P i = 53.83μW.

FIG. 4 .
FIG. 4. (a) Experimental mapping of the nonlinear dynamics in which hopping is observed in the injected power-frequency detuning plane.Time series of the (P1 & P1-both) state (b) and the (P1 & P1-both & IL+PS) state (c) are obtained, respectively, for ν i = 0.7 GHz, P i = 120.6 μW and ν i = 0.9 GHz, P i = 106.2μW.Black (red, lower) traces correspond to the parallel (orthogonal) polarization signal.(d) Representation in the orthogonal-parallel polarization plane of the P1 & P1-both & IL+PS time series.The three attractors are P1-both, P1, and IL+PS, and the arrows show the order of the sequence.Cases illustrated in Figs.5(a), 5(b) and 5(c) are indicated in the diagram of (a) by using circle, square, and star, respectively.
shows a panel of time series corresponding to different behaviors in the P1 & IL+PS region [Figs.5(a)-5(c)], zooms of the previous series [Figs.5(d)-5(f)], and their corresponding residence time distributions on each attractor in a log-log scale [Figs.5(g)-5(i)]

5 FIG. 6 .
FIG. 6. Residence time distribution of τ 2 for I bias = 5.0 mA, ν i = 1.1 GHz, and P i = 64.71μW (open dots) and data fit from τ 2 = 10 −5 to 10 −2 s giving a decay characterized by a −1.85 power law (thick dashed line).Individual distributions obtained with fewer data are plotted with colored dashed lines.Decays characterized by −3/2 and −2 power laws (thin dashed lines) are also plotted.

Figure 8 (
FIG. 8. Simulated time series of attractor hopping between CH and CH-both, obtained for I bias = 3.0 mA, ν i = −2.2GHz, P i = 0.024 (a.u), and a noise strength of (a) β SF = 2.7 × 10 −4 , and (b) β SF = 6.5 × 10 −4 .(c) and (d) are zooms of (a) on a CH and a CH-both time window, respectively, such that (c) has been shifted in 50 μs.(e) and (f) are results obtained with β SF = 1 × 10 −30 for two different sets of initial conditions.Black (red, lower) trace is the parallel (orthogonal) polarization mode signal.