- Email: [email protected]
Complex intermittent dynamics in large-aspect-ratio homogeneously broadened single-mode lasers
Physica D 203 (2005) 185–197
Complex intermittent dynamics in large-aspect-ratio homogeneously broadened single-mode lasers D. Amrouna , M. Brunela , C. Letellierb,∗ , H. Leblondc , F. Sanchezc a
b
Groupe d’Optique et d’Optronique, CORIA UMR 6614, Universit´e de Rouen, Av. de l’Universit´e, BP 12, F-76801 Saint-Etienne du Rouvray Cedex, France Groupe d’Analyse Topologique et de Mod´elisation de Syst`emes Dynamiques, CORIA UMR 6614, Universit´e de Rouen, Av. de l’Universit´e, BP 12, F-76801 Saint-Etienne du Rouvray Cedex, France c Laboratoire POMA, UMR 6136, Universit´ e d’Angers, 2 Bd Lavoisier, F-49045 Angers Cedex 01, France Received 10 October 2003; received in revised form 21 December 2004; accepted 24 March 2005 Communicated by C.K.R.T. Jones
Abstract Spatio-temporal dynamics of a homogeneously broadened single-mode laser with large Fresnel number is investigated above the second laser threshold. The system is described by the two-level Maxwell–Bloch equations. A simple quasi-periodic regime is observed when the cavity is tuned below resonance, and very uncommon dynamics is obtained when the cavity is tuned above resonance. In the latter case, the laser intensity presents “plateaux” of nearly constant values which are interrupted by bursts of large amplitude oscillations. The underlying dynamics is described in terms of heteroclinic connections between unstable periodic orbits associated with constant intensities. A very surprising characteristic of this dynamics is that the periodic orbits are always visited in a given order which is related to the sequence of wave vectors selected by the laser. This new type of behavior presents many characteristics of intermittency. © 2005 Elsevier B.V. All rights reserved. PACS: 42.55.−f ; 42.65.Sf ; 05.45.−a Keywords: Laser dynamics; Intermittency
1. Introduction The understanding of complex behaviors is directly connected to the development of the nonlinear dynam∗
Corresponding author. Tel.: +33 2 32 95 37 15; fax: +33 2 32 91 04 85. E-mail address: [email protected] (C. Letellier). 0167-2789/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.physd.2005.03.015
ical systems theory, i.e. the study of ordinary differential equations when the system is purely temporal or the study of partial derivative equations when it has a spatial dependence. Lasers are widely known for their propensity to produce a large variety of temporal chaos since the work of Haken [1], who showed that instability can occur in single-mode class-C lasers in the bad cavity configuration. Lasers are further well known for
186
D. Amroun et al. / Physica D 203 (2005) 185–197
generating spatio-temporal chaos. A Maxwell–Bloch model was first formulated by Lugiato et al. [2] to describe the phenomenon of spontaneous transverse spatial pattern formation in lasers. Since this work, many investigations have been devoted to pattern formation and related transverse effects in lasers and other systems [3–13]. The control parameter that determines the mechanisms which give rise to pattern formation is the Fresnel number [7]. Optical resonator with low Fresnel number imposes the geometry of the laser field which can be expressed in terms of empty cavity modes. The behavior of the system is less constrained by boundary conditions when the Fresnel number increases. Bulk parameters and nonlinearities of the active medium then dominate the structure formation. Homogeneously broadened two-level lasers with large Fresnel number are known to exhibit dynamical behaviors depending on the sign of the detuning of the field frequency from the atomic resonance transition frequency [5,9,10,14]. For positive detuning, the laser selects a transverse spatially homogeneous solution above the lasing threshold, whereas a traveling wave is selected for negative detuning. One of the first studies of a single-mode laser above the second laser threshold has been discussed for a good cavity configuration by Jakobsen et al. [5]. The analysis was performed using the full Maxwell–Bloch model with one transverse spatial dimension. Moreover, the stability diagram of the solutions in the plane of detuning and pumping parameter was computed. Time evolutions of the power spectrum of the electric field were also provided to describe the spatial dynamics above the second laser threshold. In the present paper, the underlying dynamics of a similar laser is investigated by using a topological analysis of the phase portraits reconstructed from the time evolution of the intensity or from the electric field. In addition to that, a first link with the spatial character of the dynamics is discussed. Above the second laser threshold, a quasiperiodic regime is observed when the cavity is tuned below resonance and a very uncommon dynamics is observed when the cavity is tuned above resonance. In particular, the time evolution of the laser intensity presents phases of nearly constant intensity which are interrupted by bursts characterized by large amplitude oscillations. This behavior resemble “cycling chaos” as described by Ashwin and co-workers [15,16] in the
sense that heteroclinic cycles between saddles are observed but the underlying dynamics is here very different. Other heteroclinic connections were also investigated using a six-dimensional approximation of the inertial manifold for the Kuramoto–Sivashinsky equation with periodic boundary conditions [17]. The main difference between the behaviors observed in the laser system investigated here and those by Dawson and Mancho is that the former behaviors are structured around tori and the latter behaviors are not. Such a deep departure comes from the continuous rotation induced by the detuning involved in the laser equations but not present in the Kuramoto–Sivashinsky equation. The subsequent part of this paper is organized as follows. In Section 2, the full Maxwell–Bloch equations with one transverse spatial dimension are described. The stability analysis of lasing solutions is also briefly discussed. In Section 3, two typical behaviors are investigated. Section 3.2 is devoted to the dynamics observed for negative detuning. Section 3.3, which is the main part of this paper, presents the very particular behavior observed for positive detunings. Such a behavior has many characteristics of intermittency and its topological structure is discussed in terms of phase portraits and first-return maps. Section 4 gives a conclusion.
2. Model equations and stationary solutions We consider here a set of two-level atoms interacting with a quasi-resonant electromagnetic wave. This system is described theoretically by the normalized Maxwell–Bloch equations, which take into account diffraction in the transverse direction x [2,3,9,5]: ∂d 1 ∗ ∗ = −γ d − r + (ep + e p) , ∂τ 2 ∂p = −(1 − iδ)p + ed, ∂τ ∂e ∂2 e = −σ(e − p) + iA 2 . ∂τ ∂x
(1)
The quantities d, p and e are obtained through normalization of the population inversion D, the macroscopic polarization P and the electric field E, respectively, ac-
D. Amroun et al. / Physica D 203 (2005) 185–197
cording to: ωN0 µ2 d= D, 0 γl γ⊥ 2ωN0 µ2 p=i P, √ 0 γl γ γ⊥ 2µ E, e = √ γ γ⊥
(2)
where µ is the off-diagonal dipole matrix element and N0 is the density of resonant atoms. γl is the inverse of the photon lifetime in the cavity, and is related to the cavity losses per round-trip; r is the pumping parameter; γ and γ⊥ are the relaxation rates of inversion and polarization, respectively; δ = (ω − ω0 )/γ⊥ is the normalized detuning of the field frequency ω from the atomic resonant frequency ω0 . We write σ = γl /(2γ⊥ ) and γ = γ /γ⊥ . The time τ is normalized to the coherence lifetime through τ = γ⊥ t. Parameter A is the diffraction parameter. In the next sections, we will consider the following set of parameters: σ = 0.01 and γ = 0.2. These parameters correspond to a good cavity configuration as in [5]. Parameter A is equal to 0.05 which corresponds to a large Fresnel number cavity. The first stationary solution of system (1) is the “off” state e = p = 0 and d = r, i.e. when the laser is below threshold. When pumping passes over the first laser ¯ p= threshold, lasing solution is of the form d = d, p¯ exp i(kx + Ωτ) and e = e¯ exp i(kx + Ωτ), where e¯ , p¯ and d¯ are given by: ¯ e¯ = r − d, i(Ω + Ak2 ) p¯ = 1 + e¯ , (3) σ 2 2 Ak + δ d¯ = 1 + . 1+σ The plane-wave pulsation Ω and the wavenumber k are linked through the dispersion relation: Ω=
δσ − Ak2 . 1+σ
187
(4)
The wavenumber k selected by the laser above the first laser threshold is obtained by seeking for the value of k for which the lasing threshold is minimum. This value √ has been demonstrated to be k = ± −δ/A when δ <
Fig. 1. Stability diagram in the plane (δ, r). The other parameters are: σ = 0.01, γ = 0.2 and A = 0.05.
0 leading to a traveling wave, and k = 0 when δ > 0, corresponding to a spatially homogeneous solution [5]. For negative detunings, the laser selects one of the two possible values of k. Fig. 1 summarizes the stability diagram of the lasing solution (as reported in [5]). The threshold of the lasing solution (first laser threshold) has been reported for sake of clarity. The instability threshold of this lasing solution has been obtained using the standard method described in [5]. With our set of parameters, when the detuning is greater than δ = 2.46, the lasing solution remains unstable over a finite range of the pumping parameter r, since this solution becomes stable again through an inverse bifurcation for higher pumping parameters, as previously mentioned by Jakobsen et al. [5]. 3. Dynamical behaviors above the second laser threshold 3.1. Introduction We have performed numerical integration of system (1) in order to examine precisely the dynamical regimes in the different instability domains of Fig. 1. The numerical method is described in Appendix A. Periodic boundary conditions are used in our computations. Numerical integrations have been first tested in the stability domains by checking that the asymptotic behavior of the Maxwell–Bloch equation (1) corresponds to a homogeneous state for δ > 0, whereas it is a traveling
188
D. Amroun et al. / Physica D 203 (2005) 185–197
wave for δ < 0, as predicted by the analytical study. In the latter case, the laser selects one of the two wave vectors k predicted by such an analytical study, which corresponds physically to an off-axis far-field emission. A plane-wave solution with k = k0 is selected since close to the threshold this mode grows the fastest [5]. We have then performed numerical integration of the system in order to study the unstable regimes beyond the second laser threshold. Initial conditions are chosen by slightly perturbing the lasing solutions of system (1) which are unstable above the second laser threshold. Two kinds of unstable regimes are presented in the next sections for two sets of values of the detuning δ and the pumping parameter r. 3.2. Negative detuning case and high pumping power We start by investigating the dynamics corresponding to the √ case where the traveling-wave solution with k0 = ± −δ/A loses its stability for negative detunings δ (Fig. 1). This occurs for high values of r (above the second laser threshold). We computed the time evolution of the electric field Re[e(x = 0, τ)] at the point x = 0 with parameter values δ = −2 and r = 100. Similar time series are obtained for any other value of x. The laser dynamics is a limit cycle in the reconstructed phase space spanned by the derivative coordinates from the laser intensity [18]. In the reconstructed phase space from Re[e(x = 0, τ)], the attractor of the system is a torus (Fig. 2), which is characteristic of a quasi-periodic dynamics. Typically, a quasi-periodic behavior results from the presence of two (or more) incommensurate frequencies, i.e. the ratio between them is irrational. Consequently, the trajectory evolves on a torus in the phase space. The Poincar´e section of a torus is obviously isotopic to a simple circle. A first-return map to a Poincar´e section of the attractor shown in Fig. 2 is used to unfold the dynamics. Such a first-return map reveals a closed curve (Fig. 3), thus confirming the quasi-periodic behavior. To check the presence of incommensurate frequencies, a temporal Fourier spectrum |ˆe(x = 0, ν)|, where eˆ (x, ν) is the Fourier transform of e(x, τ) with respect to the time variable τ, is computed from the time series of the real part of the electric field (Fig. 4). This spectrum is characterized by two main frequencies, ν1 = 0.314 and ν2 = 1.085 (in units of γ⊥ ), which
Fig. 2. Phase portrait in the reconstructed phase space using the derivative coordinates from the real part of the electric field X = ˙ Parameters are now δ = −2 and r = 100. Re[e(x = 0, τ)] (Y = X).
Fig. 3. First-return map to a Poincar´e section of the attractor reconstructed in Fig. 2.
Fig. 4. Temporal Fourier spectrum of Re[e(x = 0, τ)]. Two linear combinations of the main frequencies ν1 and ν2 are identified in this spectrum.
D. Amroun et al. / Physica D 203 (2005) 185–197
189
Fig. 6. Evolution of the laser intensity I(τ) vs. time at the point x = 0. Parameters are δ = 3 and r = 25.
Fig. 5. Spatial Fourier spectrum |˜e(k, τ)| vs. the wave vector k. It does not vary with time.
are incommensurate (as far we can conclude from a numerical estimation). The spectrum contains only two additional components which are linear combinations of the two main frequencies. They are 2ν2 − ν1 and ν2 − 2ν1 , respectively. We thus have the characteristics of a quasi-periodic regime. Note that 2πν1 corresponds to the detuning |δ|. The frequencies 2ν2 − ν1 and ν2 − 2ν1 correspond to those produced by a purely cubic nonlinearity. Let us consider now the transverse aspect of the behavior investigated here. The dynamics exhibited by the system corresponds to a solution with two wave vectors, as it is demonstrated in Fig. 5 which represents the spatial Fourier spectrum |˜e(k)| versus the wave vector k, where e˜ (k, τ) is the Fourier transform of e(x, τ) with respect to the space variable x. This spectrum is not timedependent for this regime. The main component √in this spectrum is one of the two wave vectors k = ± −δ/A corresponding to the lasing solution that loses its stability above the second threshold. The second component occurs at the second laser threshold. Each one is associated with one of the two frequencies ν1 and ν2 , in such a way that both satisfy the dispersion relation (4). In this regime, the laser emits simultaneously in two directions in the far-field region.
k0 = 0) and for relatively low pumping power. Fig. 6 shows the laser intensity I(x = 0, τ) = |e(x = 0, τ)|2 versus τ on the particular point x = 0 with δ = 3 and r = 25. Note that the regime considered here is set to 2.55 times the first laser threshold roughly equal to 10.0 (Fig. 1). Time series of the intensity observed for other spatial positions x are similar to the one shown in Fig. 6. Time intervals during which the intensity is nearly constant are interrupted by bursts of large amplitude oscillations. During each of these intervals, the dynamics is very close to a sinusoidal evolution with a unique pulsation Ω of the fields confirmed by the temporal Fourier spectrum |ˆe(x = 0, ν)|. From the time evolution of the laser intensity (Fig. 6), five different nearly sinusoidal evolutions are identified and denoted by A–E, respectively. The time evolution of the spatial Fourier spectrum |˜e(k, τ)| (Fig. 7) clearly reveals that each nearly sinu-
3.3. Positive detuning and low pumping power Let us now consider the instability occurring for positive detunings (unstable homogeneous solution with
Fig. 7. Time evolution of the components of the spatial Fourier spectrum |˜e(k, τ)| for the first part of the time series shown in Fig. 6 (0 < τ < 8 × 104 ).
190
D. Amroun et al. / Physica D 203 (2005) 185–197
Fig. 8. Pulsations |Ω| vs. the wave vectors |k| for the different traveling waves (triangles) and analytic dispersion curve (solid curve).
soidal evolution is associated with a unique wave vector k. Thus, each of them corresponds to a specific transverse traveling wave. Physically, during transitions between one traveling wave to another one, the laser spontaneously modifies its direction of emission according to the value of the wave vector k. Values of pulsation Ω (Eq. (4)) associated with traveling waves are plotted versus the corresponding wave vectors k (obtained from the corresponding spatial Fourier spectrum) in Fig. 8. The dispersion curve Ω(k) computed analytically (relation (4)) is also reported for comparison. The pulsation corresponding to each traveling wave satisfies the dispersion relation. This confirms that each of them is associated with a unique Ω and with a specific norm of wave vector |k|. At the end of each time interval during which the intensity remains constant, a competition between different wave vectors k (and thus different traveling waves) appears, but only one vector k wins and emerges during a new finite time interval corresponding to the neighborhood of a new periodic solution. It is possible to predict partially the sequence of appearance of the spatial wave vectors using the linear stability analysis described in [5]. The eigenvalues (λi ) of each traveling wave are computed versus a perturbing wave vector "k. We thus determine the value of "k corresponding to the highest growth rate, i.e. the highest Re(λ), leading to the prediction of the next traveling wave (k, Ω). In this way, the whole evolution may be predicted with a relatively good precision.
Fig. 9. Real part of the eigenvalues max[Re(λ)] vs. the perturbation wavenumber "k.
An example is presented in Fig. 9 where max[Re(λ)] versus "k is plotted. In this case, the system has been linearized around the transverse traveling-wave solution noted A corresponding to (k = 0, Ω(k = 0)). We can observe that a maximum growth rate is obtained for a perturbation wave vector "k = 2.775, which is very close to the wave vector associated with traveling wave B (k = 2.69) which is hereafter selected. The difference is smaller than the numerical discretization step δk = 0.14 × 2π. Note however that an uncertainty exists about the sign of "k. Moreover, it may appear that two different perturbation wave vectors present high growth rates. The system then selects one of these perturbation wave vectors, but it cannot be predicted which one. Such a feature is rather similar for each type of transition observed. The nature of this complex behavior is now investigated from a temporal point of view. One of the most useful results from the nonlinear dynamical systems theory is that the dynamics is better investigated in the phase space than using a single time series. Thus, the five traveling waves associated with the five plateau described by the laser intensity (Fig. 6) are clearly identified as periodic orbits when a plane projection of the phase space is used (Fig. 10). These five periodic orbits are organized around the fixed point of system (1) located at the origin of the phase space (the “off” state). All these orbits are period-1. In the phase space, the trajectory switches from the neighborhood of one orbit to the neighborhood of another one through a so-called “burst”.
D. Amroun et al. / Physica D 203 (2005) 185–197
Fig. 10. The five periodic orbits shown in the plane projection Re(e)– Im(e) of the original phase space. Bursts are omitted.
These orbits can never be reached by the trajectory since they are all unstable. According to a stability analysis, each of the five orbits is characterized by: • two complex conjugate eigenvalues with positive real parts; • two complex conjugate eigenvalues with negative real parts; • one negative real eigenvalue. These orbits are therefore saddles. A trajectory close to one of these unstable periodic orbits is ejected in spiralling within its 2D unstable manifold and approaches another periodic orbit in spiralling toward it within its 2D stable manifold (Fig. 11). In order to show these heteroclinic connections, the plane Re(e)–Im(e) should be used as for displaying the unstable periodic orbits (Fig. 10). Unfortunately, the bursts are too long with respect to the time period of one orbit (a burst between C and D is typically thousands of times longer than the time period of orbit D) and the path of the burst is blurred. Consequently, we used the intensity I = |e|2 which modds out the phase effect and plotted the trajectory in the plane spanned by u = I cos(wt) and v = I sin(wt), where pulsation w is chosen arbitrarily to ensure a representation of a subsequence CDE over 2π radians. Typically, w = T2π , where TCDE is CDE the time duration of a subsequence CDE. Bursts from C to D and from D to E are shown in Fig. 11, where periodic orbits C–E appear as simple ellipses. Each burst starts with an oscillation with an increasing amplitude.
191
Fig. 11. The (u, v)-plane projection of one subsequence CDE where the pulsations w is set to T2π to enlarge the typical structure of CDE bursts. Periodic orbits are explicitly represented as closed loops.
Once the oscillation amplitude grows in such a way that the trajectory is close to the “off” state (u = 0 and v = 0), the amplitude of oscillations decreases around the next orbit. We have further carried a study of the different bursts and intensity-constant phases over many sequences ABCDE. A global signature of the dynamics in this case can be obtained by computing a first-return map to a Poincar´e section. In order to do that, it is convenient to use the phase portrait reconstructed from the intensity using the derivative coordinates, where rotation symmetry is modded out [19,20]. Topological analysis is thus simplified since each periodic orbit becomes a fixed point. Bursts associated with the transitions are consequently evidenced. A natural choice for the Poincar´e section is to use the plane I˙ = 0, since it contains all fixed points corresponding to the periodic orbits of system (1). The first-return map thus obtained is confined along the first bisecting line. Since a Nthreturn map to a Poincar´e section may be viewed as a reconstruction of this section using a single coordinate, we then use a quite large value of N to sufficiently unfold the underlying structure. A 10th-return map is adequate to spread out the structure from the bisecting line (Fig. 12). From this map, it appears clearly that bursts between orbits B and C, C and D, D and E, E and A describe always the same path in the differential embedding from the intensity I. For theses transitions, the heteroclinic connections are unique in this space and,
192
D. Amroun et al. / Physica D 203 (2005) 185–197
Fig. 12. The 10th-return map to a Poincar´e section of the attractor reconstructed using the derivative coordinates from I(x = 0, τ). Fixed points corresponding to unstable periodic orbits A–E are drawn, respectively. They are located in the bisecting line (drawn by hand) as fixed points must be.
manifold of orbit B. Usually, the structure of such heteroclinic tangles is very complicated and the dynamics is very sensitive to initial conditions [21]. Consequently, bursts between orbits A and B may be considered as chaotic. We must however note that bursts between orbits B and C, C and D, D and E, E and A are not identical from one sequence to the other in the plane projection Re(e)–Im(e), in contrast to the intensity representation. Two different bursts between orbits C and D are shown in Fig. 14 in the Re(e)–Im(e) plane. In this figure, one burst is the symmetric of the other under a rotation around the d-axis, although they are identical in the reconstructed intensity phase space. Bursts associated with these transitions are not sensitive to initial
consequently, there is a single path between the corresponding fixed points. Only bursts between A and B evolve with different trajectories in the differential embedding from the intensity. This is exhibited by a blow up of the 10th-return map around fixed point B (Fig. 13) where the reinjections of the trajectory toward fixed point B differ at each sequence. This means that there are many heteroclinic connections between orbits A and B, i.e. there are many intersections between the 2D unstable manifold of orbit A and the 2D stable
Fig. 13. Blow up of the 10th-return map around the fixed point corresponding to orbit B. The trajectory tends toward the fixed point by spiralling in and moves away along the straight line (indicated by the arrow). All reinjections toward orbit B follow different paths.
Fig. 14. Two different transitions from orbit C to orbit D. These two transitions are invariant under a rotation symmetry. (b) is obtained from (a) by rotating through about π3 radians in the clockwise direction.
D. Amroun et al. / Physica D 203 (2005) 185–197
conditions, since they are exactly reproduced from one sequence to another. In summary, quite long phases during which the trajectory is nearly periodic are interrupted by bursts which may be sensitive to initial conditions or not. This is therefore a kind of intermittent behavior which is observed. Usually, intermittencies are classified according to the reinjection process from one laminar phase to the other and using the distribution of the laminar lengths (time duration of laminar phases). A type-II intermittency is characterized by a reinjection path which spirals around the periodic solution [22]. This is the case for the reinjection process at the end of the burst between periodic points A and B (Fig. 13) and, consequently, a type-II intermittency could be associated with these bursts. The intermittency here observed is structured around five periodic orbits. Between two successive bursts, the trajectory is nearly periodic during a time interval which fluctuates from one sequence to another. This is quite surprising since such a “laminar phase” is bounded by bursts which are identical from one se-
193
quence to the other in the differential embedding from the intensity. Such a point is not yet clearly understood. Fluctuations of laminar lengths are exhibited by computing their distributions (Fig. 15). Laminar lengths are delimited using a small threshold value of oscillation amplitude under which the intensity is assumed constant. The lengths of these laminar phases are averaged over more than 500 cycles and reported below with their standard deviation: "τA = 14696, "τB = 11939, "τC = 8302, "τD = 16891, "τE = 2523,
σA = 785, σB = 509, σC = 413, σD = 355, σE = 99,
where the times are expressed in units of time τ⊥ = 1/γ⊥ . The distributions of laminar lengths associated with orbits C–E have similar shapes and present a long queue toward the long times (Fig. 15). This is a rather characteristic feature of type-III intermittency. Note
Fig. 15. Histograms of laminar lengths associated with orbits A–E, respectively.
194
D. Amroun et al. / Physica D 203 (2005) 185–197
that the burst from C to D follows exactly the same trajectory each time. The same for D to E and E to A. Distributions of laminar phases associated with orbits A and B are not typical since they have roughly a symmetric bell shape. These phases are bounded by one side by chaotic bursts which may explain the differences observed with the distributions of laminar lengths for phases C–E. These departures in the distributions suggest that successive laminar lengths are decorrelated. This was checked by computing lengths for phases B versus the lengths for phases A, and so on for other successive phases. A cloud of points was obtained in each case thus revealing the lack of correlation between successive laminar phases in a sequence ABCDE. This would suggest that this is a multichannel intermittency [23] although the periodic orbits are always visited according to the sequence ABCDE. Thus, the dynamical behavior associated with the very particular time series shown in Fig. 6 is an intermittency structured around five periodic orbits which are visited always in the same order. The phase portrait made of the laminar phases and the bursts is quite complicated and may be viewed within the surface of a nontrivial torus. Such a torus must be embedded within a space with a dimension at least equal to 4 since heteroclinic connections involve one 2D unstable and one 2D stable manifolds. This is in agreement with the embed-
Fig. 16. Intensities associated with the unstable periodic orbits vs. the pumping parameter r. Other parameter values: δ = 3, σ = 0.01, γ = 0.2 and A = 0.05. Periodic orbits are detected by a threshold related to the amplitude of oscillations under which the intensity is assumed to be constant.
ding dimension we estimated to 4 using the algorithm written by Cao [24]. Such intermittent behavior can be observed over a quite large interval of the pumping parameter but the number of unstable periodic orbits involved varies. In order to have a global view of how this number depends on the pumping parameter, we plotted the intensities associated with the unstable periodic or-
Fig. 17. Two windows of a time series of the laser intensity with r = 37. Other parameter values as for Fig. 16.
D. Amroun et al. / Physica D 203 (2005) 185–197
bits versus r (Fig. 16). Roughly, the number of periodic orbits increases versus parameter r. Thus, there are three orbits for r ∈ [15; 18], four for r ∈ [18; 24], five for r ∈ [24; 35]. For greater values, up to six unstable periodic orbits are observed (r ∈ [35; 38]) and then decreases to five. For r > 41, a single stable periodic orbit remains because the lasing solution is stable again, according to the stability diagram shown in Fig. 1. The number of unstable periodic orbits can increase or decrease. It increases when a burst is splitted into two bursts when the pumping parameter is varied. Once these two bursts are sufficiently separated, a new plateau is created. Plateaux disappear by the reverse process. A time series where six plateau are clearly identified is shown in Fig. 17a for r = 37. We designate the sixth plateau by the letter B to indicate that it comes from the burst at the end of the plateau B. The chaotic character of the bursts between B and C implies very different lengths for plateau B . Thus, in a different window of the time series computed for r = 37, the plateau B is skipped since the two bursts are too close (Fig. 17b). This time series well exemplifies the way in which plateaux—or equivalently, the unstable periodic orbits—appear or disappear.
195
finite number of directions, and in a well-defined sequence. In the phase space reconstructed from the intensity where the symmetry associated with the phase of the electric field is modded out, it is shown that a single heteroclinic connection has a chaotic character, the others being associated with a unique path. This very uncommon behavior presents many characteristics of intermittency. In particular, the reinjection process near the periodic orbit visited after chaotic bursts is typical of type-II intermittency. It is expected that such a new type of dynamics could be observed in other spatio-temporal systems. Obviously, the number of periodic orbits around which the dynamics is structured is parameter-dependent but the global mechanism should be quite general. It seems that such spatio-temporal behaviors could be obtained for systems having traveling waves for solutions. Thus, perhaps the most appropriate experimental device for observing the complex intermittent behavior here reported could correspond to the semiconductor laser device where experimental traveling waves have been observed [27]. In the unstable domain of its parameter space, similar intermittent behaviors could be reproduced.
Acknowledgement 4. Conclusion Using the two-level Maxwell–Bloch equations with diffraction, complex dynamical behaviors of a homogeneously broadened single-mode laser with large Fresnel number have been investigated above the second laser threshold. When the cavity detuning δ is negative, there is a bi-angular laser emission in the farfield region. The underlying dynamics is quasi-periodic with two main frequencies. However, for positive cavity detuning δ, the dynamics is more complex. With the parameter values investigated here, a very particular evolution of the laser intensity has been observed. It has been shown that the dynamics is structured around five unstable periodic orbits. Heteroclinic connections between them are associated with bursts which interrupt phases during which the laser intensity is nearly constant. We showed that the order in which periodic orbits are visited is always the same. This is related to a sequence of five wave vectors. Thus, the laser spontaneously modifies its direction of emission in a
We wish to thank Robert Gilmore for helpful comments and encouraging this work.
Appendix A We describe in this section the numerical method used to solve the Maxwell–Bloch system (1). It is an adaptation of the Picard iteration method [25]. The latter is frequently used to prove the existence of solutions of the Cauchy problem for hyperbolic partial derivative equations (PDEs). In addition, the equation is split into a linear and a nonlinear part. In fact, only the nonlinear part is solved by means of Picard iterations, while the linear part is solved analytically by means of a Fourier transform with regard to the space variable. In the present situation, this modified scheme runs reasonably fast, and is very stable. A previous implementation of a scheme of this type is presented in [26]. We need to solve the Cauchy problem for the following system
196
D. Amroun et al. / Physica D 203 (2005) 185–197
of PDEs: 1 ∗ ∂d
∗ = −γ d + + e p) , (ep ∂t 2 ∂p = −(1 − iδ)p + e(d + r), ∂t 2 ∂e = −σ(e − p) + iA ∂ e , ∂t ∂x2
According to Picard’s method, we put Eq. (A.5) into the integral form: t
ˆ ˆ exp(−Lt )B(u(t )) dt , u(t) ˆ = exp(Lt) uˆ i + (A.8)
where we put d = d − r for sake of simplicity. This system can be split into the sum of a linear and a nonlinear part, as: ∂u = L(u) + B(u), (A.2) ∂t with −γd d −(1 − iδ)p + re u = p , L(u) = , ∂2 e e −σ(e − p) + iA 2 ∂x (A.3) γ − (ep∗ + e∗ p) 2 . B(u) = ed
(A.4)
0 Thus, L is a linear operator, while B is quadratic. We denote by u(k, ˆ t) the Fourier transform of u(x, t) with regard to the space variable x. We chose "x = 7 and δx = 0.1, where "x is the total length of the transverse dimension x and δx the sampling step. This fixes the maximum value of the wave vector to kmax = 10π. Then, the Fourier transform of Eq. (A.2) reads:
σ
where u(t ˆ = 0) = uˆ i is the initial data. An important point of the method is the analytical ˆ We computation of the exponential of the matrix Lt. find:
e−γ∗t
0 ˆ exp(Lt)= 0
0
0
e−Xt/2 Rt 2r sinh R 2 , −Xt/2 −Xt/2 Rt e e 2s sinh Q R R 2 (A.9) e−Xt/2 Q R
with X = 1 − iδ + ik2 A + σ, R= 2 2 A + σ)] X + 4[σr − (1 − iδ)(ik and Q = Rt R cosh Rt + (X − 2(1 − iδ)) sinh . Now we 2 2 have to solve Eq. (A.5) with the initial data uˆ i . We build, following the Picard iteration scheme, the recurrent sequence: • uˆ 0 is given (in practice, uˆ 0 (t) ≡ uˆ i ), • uˆ n known, uˆ n+1 is the solution of: ∂uˆ n+1 = A ˆ uˆ n+1 + B(u n ), ∂t uˆ n+1 (t = 0) = uˆ i .
(A.10)
It follows from Eq. (A.8) that: t
ˆ ˆ exp(−Lt )B(un (t )) dt . uˆ n+1 = exp(Lt) uˆ i + 0
∂uˆ ˆ uˆ + B(u), =L ∂t with −γ 0 ˆ = 0 −1 + iδ L 0
0
(A.1)
(A.5)
0 r −σ
,
(A.6)
− iAk2
γ ∗ p) − (e p∗ + e 2 . = B(u) ed 0
(A.7)
(A.11) Since at this stage B(u n ) is a known function, the computation of uˆ n+1 reduces to that of an integral. As soon as the sequence (u) ˆ n converges to a fixed point, its limit uˆ is the solution of Eq. (A.5). By computation of the inverse Fourier transform, we get u, the solution of system (A.1). In fact, the sequence uˆ n converges on some time interval [0, T ], where T is bounded from above. Thus, to calculate the solution for a longer range (a time series), we proceed as follows: we choose t0 < T and apply the
D. Amroun et al. / Physica D 203 (2005) 185–197
previous algorithm which gives u(t) in the range [0, t0 ]. Then, we apply it again with ui = u(t0 ) to obtain u in the interval [t0 , 2t0 ], and so on (as long as T is greater than t0 ). The main interest of this procedure is that the fixed point is reached with some chosen precision, on the whole interval [0, t0 ]. This way, we are able to control directly and quantitatively the difference between the numerical and the exact solution in the L∞ -norm.
References [1] H. Haken, Analogy between higher instabilities in fluids and lasers, Phys. Lett. A 53 (1) (1975) 77–78. [2] L.A. Lugiato, C. Oldano, L.M. Narducci, Cooperative frequency locking and stationary spatial structures in lasers, J. Opt. Soc. Am. B 5 (5) (1988) 879–888. [3] P. Coullet, L. Gil, F. Rocca, Optical vortices, Opt. Commun. 73 (5) (1989) 403–408. [4] N.B. Abraham, W.J. Firth, Special issue on transverse effects in optical systems, J. Opt. Soc. Am. 7 (1990) 951. [5] P.K. Jakobsen, J.V. Moloney, A.C. Newell, R. Indik, Space– time dynamics of wide-gain-section lasers, Phys. Rev. A 45 (1992) 8129–8137. [6] D. Dangoisse, D. Hennequin, C. Lepers, E. Louvergneaux, P. Glorieux, Two-dimensional optical lattices in a CO2 laser, Phys. Rev. A 46 (1992) 5955. [7] F.T. Arecchi, S. Boccaletti, P.L. Ramazza, S. Residori, Transition from boundary- to bulk-controlled regimes in optical pattern formation, Phys. Rev. Lett. 70 (1993) 2277. [8] M.C. Cross, P.H. Hohenberg, Pattern formation outside of equilibrium, Rev. Mod. Phys. 65 (1993) 851. [9] P.K. Jakobsen, J. Lega, Q. Feng, M. Staley, J.V. Moloney, A.C. Newell, Nonlinear transverse modes of large-aspect-ratio homogeneously broadened lasers. I. Analysis and numerical simulation, Phys. Rev. A 49 (1994) 4189. [10] J. Lega, P.K. Jakobsen, J.V. Moloney, A.C. Newell, Nonlinear transverse modes of large-aspect-ratio homogeneously broadened lasers. II. Pattern analysis near and beyond threshold, Phys. Rev. A 49 (1994) 4201.
197
[11] G. Huyet, M.C. Martinoni, J.R. Tredicce, S. Rica, Spatiotemporal dynamics of lasers with a large Fresnel number, Phys. Rev. Lett. 75 (1995) 4027–4030. [12] F. Encinas-Sanz, I. Leyva, J.M. Guerra, Time resolved pattern evolution in a large aperture laser, Phys. Rev. Lett. 84 (2000) 883–886. [13] Y.F. Chen, Y.P. Lan, Spontaneous transverse pattern formation in a microchip laser excited by a doughnut pump profile, Appl. Phys. B 75 (2002) 453–456. [14] D. Amroun, M. Brunel, H. Leblond, F. Sanchez, Space–time laser instabilities in homogeneously broadened dense media, J. Mod. Opt. 50 (2003) 1487–1495. [15] P. Ashwin, A.M. Rucklidge, Cycling chaos: its creation, persistence and loss of stability in a model of nonlinear magnetoconvection, Physica D 122 (1998) 134–154. [16] P. Ashwin, A.M. Rucklidge, R. Sturman, Infinities of stable periodic orbits in systems of coupled oscillators, Phys. Rev. E 66 (2002) 035201(R). [17] S.P. Dawson, A.M. Mancho, Collections of heteroclinic cycles in the Kuramoto–Sivashinsky equation, Physica D 100 (1997) 231–256. [18] N.H. Packard, J.P. Crutchfield, J.D. Farmer, R.S. Shaw, Geometry from a time series, Phys. Rev. Lett. 45 (9) (1980) 712–716. [19] C. Letellier, R. Gilmore, Covering dynamical systems: two-fold covers, Phys. Rev. E 63 (2001) 016206. [20] C. Letellier, Modding out a continuous rotation symmetry for disentangling a laser dynamics, Int. J. Bifurcation Chaos 13 (6) (2003) 1573–1577. [21] S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos, Springer-Verlag, 1990. [22] J.-Y. Huang, J.-J. Kim, Type-II intermittency in a coupled nonlinear oscillator: experimental observation, Phys. Rev. A 36 (1987) 1495–1497. [23] C. Letellier, P. Werny, J.-M. Malasoma, R. Gilmore, Multichannel intermittencies induced by symmetries, Phys. Rev. E 66 (2002) 036220. [24] L. Cao, Practical method for determining the minimum embedding dimension of a scalar time series, Physica D 110 (1997) 43–52. [25] F.B. Hildebrand, Introduction to Numerical Analysis, second ed., McGraw-Hill, New Delhi, 1974. [26] H. Leblond, Transverse stability of solitons and moving domain walls, J. Phys. A 33 (2000) 8105. [27] J. Houlihan, J.R. O’Callaghan, V. Voignier, G. Huyet, J.G. McInerney, B. Corbett, Experimental observation of traveling waves in a semiconductor laser device, Opt. Lett. 26 (2001) 1556–1558.
Complex intermittent dynamics in large-aspect-ratio homogeneously broadened single-mode lasers D. Amrouna , M. Brunela , C. Letellierb,∗ , H. Leblondc , F. Sanchezc a
b
Groupe d’Optique et d’Optronique, CORIA UMR 6614, Universit´e de Rouen, Av. de l’Universit´e, BP 12, F-76801 Saint-Etienne du Rouvray Cedex, France Groupe d’Analyse Topologique et de Mod´elisation de Syst`emes Dynamiques, CORIA UMR 6614, Universit´e de Rouen, Av. de l’Universit´e, BP 12, F-76801 Saint-Etienne du Rouvray Cedex, France c Laboratoire POMA, UMR 6136, Universit´ e d’Angers, 2 Bd Lavoisier, F-49045 Angers Cedex 01, France Received 10 October 2003; received in revised form 21 December 2004; accepted 24 March 2005 Communicated by C.K.R.T. Jones
Abstract Spatio-temporal dynamics of a homogeneously broadened single-mode laser with large Fresnel number is investigated above the second laser threshold. The system is described by the two-level Maxwell–Bloch equations. A simple quasi-periodic regime is observed when the cavity is tuned below resonance, and very uncommon dynamics is obtained when the cavity is tuned above resonance. In the latter case, the laser intensity presents “plateaux” of nearly constant values which are interrupted by bursts of large amplitude oscillations. The underlying dynamics is described in terms of heteroclinic connections between unstable periodic orbits associated with constant intensities. A very surprising characteristic of this dynamics is that the periodic orbits are always visited in a given order which is related to the sequence of wave vectors selected by the laser. This new type of behavior presents many characteristics of intermittency. © 2005 Elsevier B.V. All rights reserved. PACS: 42.55.−f ; 42.65.Sf ; 05.45.−a Keywords: Laser dynamics; Intermittency
1. Introduction The understanding of complex behaviors is directly connected to the development of the nonlinear dynam∗
Corresponding author. Tel.: +33 2 32 95 37 15; fax: +33 2 32 91 04 85. E-mail address: [email protected] (C. Letellier). 0167-2789/$ – see front matter © 2005 Elsevier B.V. All rights reserved. doi:10.1016/j.physd.2005.03.015
ical systems theory, i.e. the study of ordinary differential equations when the system is purely temporal or the study of partial derivative equations when it has a spatial dependence. Lasers are widely known for their propensity to produce a large variety of temporal chaos since the work of Haken [1], who showed that instability can occur in single-mode class-C lasers in the bad cavity configuration. Lasers are further well known for
186
D. Amroun et al. / Physica D 203 (2005) 185–197
generating spatio-temporal chaos. A Maxwell–Bloch model was first formulated by Lugiato et al. [2] to describe the phenomenon of spontaneous transverse spatial pattern formation in lasers. Since this work, many investigations have been devoted to pattern formation and related transverse effects in lasers and other systems [3–13]. The control parameter that determines the mechanisms which give rise to pattern formation is the Fresnel number [7]. Optical resonator with low Fresnel number imposes the geometry of the laser field which can be expressed in terms of empty cavity modes. The behavior of the system is less constrained by boundary conditions when the Fresnel number increases. Bulk parameters and nonlinearities of the active medium then dominate the structure formation. Homogeneously broadened two-level lasers with large Fresnel number are known to exhibit dynamical behaviors depending on the sign of the detuning of the field frequency from the atomic resonance transition frequency [5,9,10,14]. For positive detuning, the laser selects a transverse spatially homogeneous solution above the lasing threshold, whereas a traveling wave is selected for negative detuning. One of the first studies of a single-mode laser above the second laser threshold has been discussed for a good cavity configuration by Jakobsen et al. [5]. The analysis was performed using the full Maxwell–Bloch model with one transverse spatial dimension. Moreover, the stability diagram of the solutions in the plane of detuning and pumping parameter was computed. Time evolutions of the power spectrum of the electric field were also provided to describe the spatial dynamics above the second laser threshold. In the present paper, the underlying dynamics of a similar laser is investigated by using a topological analysis of the phase portraits reconstructed from the time evolution of the intensity or from the electric field. In addition to that, a first link with the spatial character of the dynamics is discussed. Above the second laser threshold, a quasiperiodic regime is observed when the cavity is tuned below resonance and a very uncommon dynamics is observed when the cavity is tuned above resonance. In particular, the time evolution of the laser intensity presents phases of nearly constant intensity which are interrupted by bursts characterized by large amplitude oscillations. This behavior resemble “cycling chaos” as described by Ashwin and co-workers [15,16] in the
sense that heteroclinic cycles between saddles are observed but the underlying dynamics is here very different. Other heteroclinic connections were also investigated using a six-dimensional approximation of the inertial manifold for the Kuramoto–Sivashinsky equation with periodic boundary conditions [17]. The main difference between the behaviors observed in the laser system investigated here and those by Dawson and Mancho is that the former behaviors are structured around tori and the latter behaviors are not. Such a deep departure comes from the continuous rotation induced by the detuning involved in the laser equations but not present in the Kuramoto–Sivashinsky equation. The subsequent part of this paper is organized as follows. In Section 2, the full Maxwell–Bloch equations with one transverse spatial dimension are described. The stability analysis of lasing solutions is also briefly discussed. In Section 3, two typical behaviors are investigated. Section 3.2 is devoted to the dynamics observed for negative detuning. Section 3.3, which is the main part of this paper, presents the very particular behavior observed for positive detunings. Such a behavior has many characteristics of intermittency and its topological structure is discussed in terms of phase portraits and first-return maps. Section 4 gives a conclusion.
2. Model equations and stationary solutions We consider here a set of two-level atoms interacting with a quasi-resonant electromagnetic wave. This system is described theoretically by the normalized Maxwell–Bloch equations, which take into account diffraction in the transverse direction x [2,3,9,5]: ∂d 1 ∗ ∗ = −γ d − r + (ep + e p) , ∂τ 2 ∂p = −(1 − iδ)p + ed, ∂τ ∂e ∂2 e = −σ(e − p) + iA 2 . ∂τ ∂x
(1)
The quantities d, p and e are obtained through normalization of the population inversion D, the macroscopic polarization P and the electric field E, respectively, ac-
D. Amroun et al. / Physica D 203 (2005) 185–197
cording to: ωN0 µ2 d= D, 0 γl γ⊥ 2ωN0 µ2 p=i P, √ 0 γl γ γ⊥ 2µ E, e = √ γ γ⊥
(2)
where µ is the off-diagonal dipole matrix element and N0 is the density of resonant atoms. γl is the inverse of the photon lifetime in the cavity, and is related to the cavity losses per round-trip; r is the pumping parameter; γ and γ⊥ are the relaxation rates of inversion and polarization, respectively; δ = (ω − ω0 )/γ⊥ is the normalized detuning of the field frequency ω from the atomic resonant frequency ω0 . We write σ = γl /(2γ⊥ ) and γ = γ /γ⊥ . The time τ is normalized to the coherence lifetime through τ = γ⊥ t. Parameter A is the diffraction parameter. In the next sections, we will consider the following set of parameters: σ = 0.01 and γ = 0.2. These parameters correspond to a good cavity configuration as in [5]. Parameter A is equal to 0.05 which corresponds to a large Fresnel number cavity. The first stationary solution of system (1) is the “off” state e = p = 0 and d = r, i.e. when the laser is below threshold. When pumping passes over the first laser ¯ p= threshold, lasing solution is of the form d = d, p¯ exp i(kx + Ωτ) and e = e¯ exp i(kx + Ωτ), where e¯ , p¯ and d¯ are given by: ¯ e¯ = r − d, i(Ω + Ak2 ) p¯ = 1 + e¯ , (3) σ 2 2 Ak + δ d¯ = 1 + . 1+σ The plane-wave pulsation Ω and the wavenumber k are linked through the dispersion relation: Ω=
δσ − Ak2 . 1+σ
187
(4)
The wavenumber k selected by the laser above the first laser threshold is obtained by seeking for the value of k for which the lasing threshold is minimum. This value √ has been demonstrated to be k = ± −δ/A when δ <
Fig. 1. Stability diagram in the plane (δ, r). The other parameters are: σ = 0.01, γ = 0.2 and A = 0.05.
0 leading to a traveling wave, and k = 0 when δ > 0, corresponding to a spatially homogeneous solution [5]. For negative detunings, the laser selects one of the two possible values of k. Fig. 1 summarizes the stability diagram of the lasing solution (as reported in [5]). The threshold of the lasing solution (first laser threshold) has been reported for sake of clarity. The instability threshold of this lasing solution has been obtained using the standard method described in [5]. With our set of parameters, when the detuning is greater than δ = 2.46, the lasing solution remains unstable over a finite range of the pumping parameter r, since this solution becomes stable again through an inverse bifurcation for higher pumping parameters, as previously mentioned by Jakobsen et al. [5]. 3. Dynamical behaviors above the second laser threshold 3.1. Introduction We have performed numerical integration of system (1) in order to examine precisely the dynamical regimes in the different instability domains of Fig. 1. The numerical method is described in Appendix A. Periodic boundary conditions are used in our computations. Numerical integrations have been first tested in the stability domains by checking that the asymptotic behavior of the Maxwell–Bloch equation (1) corresponds to a homogeneous state for δ > 0, whereas it is a traveling
188
D. Amroun et al. / Physica D 203 (2005) 185–197
wave for δ < 0, as predicted by the analytical study. In the latter case, the laser selects one of the two wave vectors k predicted by such an analytical study, which corresponds physically to an off-axis far-field emission. A plane-wave solution with k = k0 is selected since close to the threshold this mode grows the fastest [5]. We have then performed numerical integration of the system in order to study the unstable regimes beyond the second laser threshold. Initial conditions are chosen by slightly perturbing the lasing solutions of system (1) which are unstable above the second laser threshold. Two kinds of unstable regimes are presented in the next sections for two sets of values of the detuning δ and the pumping parameter r. 3.2. Negative detuning case and high pumping power We start by investigating the dynamics corresponding to the √ case where the traveling-wave solution with k0 = ± −δ/A loses its stability for negative detunings δ (Fig. 1). This occurs for high values of r (above the second laser threshold). We computed the time evolution of the electric field Re[e(x = 0, τ)] at the point x = 0 with parameter values δ = −2 and r = 100. Similar time series are obtained for any other value of x. The laser dynamics is a limit cycle in the reconstructed phase space spanned by the derivative coordinates from the laser intensity [18]. In the reconstructed phase space from Re[e(x = 0, τ)], the attractor of the system is a torus (Fig. 2), which is characteristic of a quasi-periodic dynamics. Typically, a quasi-periodic behavior results from the presence of two (or more) incommensurate frequencies, i.e. the ratio between them is irrational. Consequently, the trajectory evolves on a torus in the phase space. The Poincar´e section of a torus is obviously isotopic to a simple circle. A first-return map to a Poincar´e section of the attractor shown in Fig. 2 is used to unfold the dynamics. Such a first-return map reveals a closed curve (Fig. 3), thus confirming the quasi-periodic behavior. To check the presence of incommensurate frequencies, a temporal Fourier spectrum |ˆe(x = 0, ν)|, where eˆ (x, ν) is the Fourier transform of e(x, τ) with respect to the time variable τ, is computed from the time series of the real part of the electric field (Fig. 4). This spectrum is characterized by two main frequencies, ν1 = 0.314 and ν2 = 1.085 (in units of γ⊥ ), which
Fig. 2. Phase portrait in the reconstructed phase space using the derivative coordinates from the real part of the electric field X = ˙ Parameters are now δ = −2 and r = 100. Re[e(x = 0, τ)] (Y = X).
Fig. 3. First-return map to a Poincar´e section of the attractor reconstructed in Fig. 2.
Fig. 4. Temporal Fourier spectrum of Re[e(x = 0, τ)]. Two linear combinations of the main frequencies ν1 and ν2 are identified in this spectrum.
D. Amroun et al. / Physica D 203 (2005) 185–197
189
Fig. 6. Evolution of the laser intensity I(τ) vs. time at the point x = 0. Parameters are δ = 3 and r = 25.
Fig. 5. Spatial Fourier spectrum |˜e(k, τ)| vs. the wave vector k. It does not vary with time.
are incommensurate (as far we can conclude from a numerical estimation). The spectrum contains only two additional components which are linear combinations of the two main frequencies. They are 2ν2 − ν1 and ν2 − 2ν1 , respectively. We thus have the characteristics of a quasi-periodic regime. Note that 2πν1 corresponds to the detuning |δ|. The frequencies 2ν2 − ν1 and ν2 − 2ν1 correspond to those produced by a purely cubic nonlinearity. Let us consider now the transverse aspect of the behavior investigated here. The dynamics exhibited by the system corresponds to a solution with two wave vectors, as it is demonstrated in Fig. 5 which represents the spatial Fourier spectrum |˜e(k)| versus the wave vector k, where e˜ (k, τ) is the Fourier transform of e(x, τ) with respect to the space variable x. This spectrum is not timedependent for this regime. The main component √in this spectrum is one of the two wave vectors k = ± −δ/A corresponding to the lasing solution that loses its stability above the second threshold. The second component occurs at the second laser threshold. Each one is associated with one of the two frequencies ν1 and ν2 , in such a way that both satisfy the dispersion relation (4). In this regime, the laser emits simultaneously in two directions in the far-field region.
k0 = 0) and for relatively low pumping power. Fig. 6 shows the laser intensity I(x = 0, τ) = |e(x = 0, τ)|2 versus τ on the particular point x = 0 with δ = 3 and r = 25. Note that the regime considered here is set to 2.55 times the first laser threshold roughly equal to 10.0 (Fig. 1). Time series of the intensity observed for other spatial positions x are similar to the one shown in Fig. 6. Time intervals during which the intensity is nearly constant are interrupted by bursts of large amplitude oscillations. During each of these intervals, the dynamics is very close to a sinusoidal evolution with a unique pulsation Ω of the fields confirmed by the temporal Fourier spectrum |ˆe(x = 0, ν)|. From the time evolution of the laser intensity (Fig. 6), five different nearly sinusoidal evolutions are identified and denoted by A–E, respectively. The time evolution of the spatial Fourier spectrum |˜e(k, τ)| (Fig. 7) clearly reveals that each nearly sinu-
3.3. Positive detuning and low pumping power Let us now consider the instability occurring for positive detunings (unstable homogeneous solution with
Fig. 7. Time evolution of the components of the spatial Fourier spectrum |˜e(k, τ)| for the first part of the time series shown in Fig. 6 (0 < τ < 8 × 104 ).
190
D. Amroun et al. / Physica D 203 (2005) 185–197
Fig. 8. Pulsations |Ω| vs. the wave vectors |k| for the different traveling waves (triangles) and analytic dispersion curve (solid curve).
soidal evolution is associated with a unique wave vector k. Thus, each of them corresponds to a specific transverse traveling wave. Physically, during transitions between one traveling wave to another one, the laser spontaneously modifies its direction of emission according to the value of the wave vector k. Values of pulsation Ω (Eq. (4)) associated with traveling waves are plotted versus the corresponding wave vectors k (obtained from the corresponding spatial Fourier spectrum) in Fig. 8. The dispersion curve Ω(k) computed analytically (relation (4)) is also reported for comparison. The pulsation corresponding to each traveling wave satisfies the dispersion relation. This confirms that each of them is associated with a unique Ω and with a specific norm of wave vector |k|. At the end of each time interval during which the intensity remains constant, a competition between different wave vectors k (and thus different traveling waves) appears, but only one vector k wins and emerges during a new finite time interval corresponding to the neighborhood of a new periodic solution. It is possible to predict partially the sequence of appearance of the spatial wave vectors using the linear stability analysis described in [5]. The eigenvalues (λi ) of each traveling wave are computed versus a perturbing wave vector "k. We thus determine the value of "k corresponding to the highest growth rate, i.e. the highest Re(λ), leading to the prediction of the next traveling wave (k, Ω). In this way, the whole evolution may be predicted with a relatively good precision.
Fig. 9. Real part of the eigenvalues max[Re(λ)] vs. the perturbation wavenumber "k.
An example is presented in Fig. 9 where max[Re(λ)] versus "k is plotted. In this case, the system has been linearized around the transverse traveling-wave solution noted A corresponding to (k = 0, Ω(k = 0)). We can observe that a maximum growth rate is obtained for a perturbation wave vector "k = 2.775, which is very close to the wave vector associated with traveling wave B (k = 2.69) which is hereafter selected. The difference is smaller than the numerical discretization step δk = 0.14 × 2π. Note however that an uncertainty exists about the sign of "k. Moreover, it may appear that two different perturbation wave vectors present high growth rates. The system then selects one of these perturbation wave vectors, but it cannot be predicted which one. Such a feature is rather similar for each type of transition observed. The nature of this complex behavior is now investigated from a temporal point of view. One of the most useful results from the nonlinear dynamical systems theory is that the dynamics is better investigated in the phase space than using a single time series. Thus, the five traveling waves associated with the five plateau described by the laser intensity (Fig. 6) are clearly identified as periodic orbits when a plane projection of the phase space is used (Fig. 10). These five periodic orbits are organized around the fixed point of system (1) located at the origin of the phase space (the “off” state). All these orbits are period-1. In the phase space, the trajectory switches from the neighborhood of one orbit to the neighborhood of another one through a so-called “burst”.
D. Amroun et al. / Physica D 203 (2005) 185–197
Fig. 10. The five periodic orbits shown in the plane projection Re(e)– Im(e) of the original phase space. Bursts are omitted.
These orbits can never be reached by the trajectory since they are all unstable. According to a stability analysis, each of the five orbits is characterized by: • two complex conjugate eigenvalues with positive real parts; • two complex conjugate eigenvalues with negative real parts; • one negative real eigenvalue. These orbits are therefore saddles. A trajectory close to one of these unstable periodic orbits is ejected in spiralling within its 2D unstable manifold and approaches another periodic orbit in spiralling toward it within its 2D stable manifold (Fig. 11). In order to show these heteroclinic connections, the plane Re(e)–Im(e) should be used as for displaying the unstable periodic orbits (Fig. 10). Unfortunately, the bursts are too long with respect to the time period of one orbit (a burst between C and D is typically thousands of times longer than the time period of orbit D) and the path of the burst is blurred. Consequently, we used the intensity I = |e|2 which modds out the phase effect and plotted the trajectory in the plane spanned by u = I cos(wt) and v = I sin(wt), where pulsation w is chosen arbitrarily to ensure a representation of a subsequence CDE over 2π radians. Typically, w = T2π , where TCDE is CDE the time duration of a subsequence CDE. Bursts from C to D and from D to E are shown in Fig. 11, where periodic orbits C–E appear as simple ellipses. Each burst starts with an oscillation with an increasing amplitude.
191
Fig. 11. The (u, v)-plane projection of one subsequence CDE where the pulsations w is set to T2π to enlarge the typical structure of CDE bursts. Periodic orbits are explicitly represented as closed loops.
Once the oscillation amplitude grows in such a way that the trajectory is close to the “off” state (u = 0 and v = 0), the amplitude of oscillations decreases around the next orbit. We have further carried a study of the different bursts and intensity-constant phases over many sequences ABCDE. A global signature of the dynamics in this case can be obtained by computing a first-return map to a Poincar´e section. In order to do that, it is convenient to use the phase portrait reconstructed from the intensity using the derivative coordinates, where rotation symmetry is modded out [19,20]. Topological analysis is thus simplified since each periodic orbit becomes a fixed point. Bursts associated with the transitions are consequently evidenced. A natural choice for the Poincar´e section is to use the plane I˙ = 0, since it contains all fixed points corresponding to the periodic orbits of system (1). The first-return map thus obtained is confined along the first bisecting line. Since a Nthreturn map to a Poincar´e section may be viewed as a reconstruction of this section using a single coordinate, we then use a quite large value of N to sufficiently unfold the underlying structure. A 10th-return map is adequate to spread out the structure from the bisecting line (Fig. 12). From this map, it appears clearly that bursts between orbits B and C, C and D, D and E, E and A describe always the same path in the differential embedding from the intensity I. For theses transitions, the heteroclinic connections are unique in this space and,
192
D. Amroun et al. / Physica D 203 (2005) 185–197
Fig. 12. The 10th-return map to a Poincar´e section of the attractor reconstructed using the derivative coordinates from I(x = 0, τ). Fixed points corresponding to unstable periodic orbits A–E are drawn, respectively. They are located in the bisecting line (drawn by hand) as fixed points must be.
manifold of orbit B. Usually, the structure of such heteroclinic tangles is very complicated and the dynamics is very sensitive to initial conditions [21]. Consequently, bursts between orbits A and B may be considered as chaotic. We must however note that bursts between orbits B and C, C and D, D and E, E and A are not identical from one sequence to the other in the plane projection Re(e)–Im(e), in contrast to the intensity representation. Two different bursts between orbits C and D are shown in Fig. 14 in the Re(e)–Im(e) plane. In this figure, one burst is the symmetric of the other under a rotation around the d-axis, although they are identical in the reconstructed intensity phase space. Bursts associated with these transitions are not sensitive to initial
consequently, there is a single path between the corresponding fixed points. Only bursts between A and B evolve with different trajectories in the differential embedding from the intensity. This is exhibited by a blow up of the 10th-return map around fixed point B (Fig. 13) where the reinjections of the trajectory toward fixed point B differ at each sequence. This means that there are many heteroclinic connections between orbits A and B, i.e. there are many intersections between the 2D unstable manifold of orbit A and the 2D stable
Fig. 13. Blow up of the 10th-return map around the fixed point corresponding to orbit B. The trajectory tends toward the fixed point by spiralling in and moves away along the straight line (indicated by the arrow). All reinjections toward orbit B follow different paths.
Fig. 14. Two different transitions from orbit C to orbit D. These two transitions are invariant under a rotation symmetry. (b) is obtained from (a) by rotating through about π3 radians in the clockwise direction.
D. Amroun et al. / Physica D 203 (2005) 185–197
conditions, since they are exactly reproduced from one sequence to another. In summary, quite long phases during which the trajectory is nearly periodic are interrupted by bursts which may be sensitive to initial conditions or not. This is therefore a kind of intermittent behavior which is observed. Usually, intermittencies are classified according to the reinjection process from one laminar phase to the other and using the distribution of the laminar lengths (time duration of laminar phases). A type-II intermittency is characterized by a reinjection path which spirals around the periodic solution [22]. This is the case for the reinjection process at the end of the burst between periodic points A and B (Fig. 13) and, consequently, a type-II intermittency could be associated with these bursts. The intermittency here observed is structured around five periodic orbits. Between two successive bursts, the trajectory is nearly periodic during a time interval which fluctuates from one sequence to another. This is quite surprising since such a “laminar phase” is bounded by bursts which are identical from one se-
193
quence to the other in the differential embedding from the intensity. Such a point is not yet clearly understood. Fluctuations of laminar lengths are exhibited by computing their distributions (Fig. 15). Laminar lengths are delimited using a small threshold value of oscillation amplitude under which the intensity is assumed constant. The lengths of these laminar phases are averaged over more than 500 cycles and reported below with their standard deviation: "τA = 14696, "τB = 11939, "τC = 8302, "τD = 16891, "τE = 2523,
σA = 785, σB = 509, σC = 413, σD = 355, σE = 99,
where the times are expressed in units of time τ⊥ = 1/γ⊥ . The distributions of laminar lengths associated with orbits C–E have similar shapes and present a long queue toward the long times (Fig. 15). This is a rather characteristic feature of type-III intermittency. Note
Fig. 15. Histograms of laminar lengths associated with orbits A–E, respectively.
194
D. Amroun et al. / Physica D 203 (2005) 185–197
that the burst from C to D follows exactly the same trajectory each time. The same for D to E and E to A. Distributions of laminar phases associated with orbits A and B are not typical since they have roughly a symmetric bell shape. These phases are bounded by one side by chaotic bursts which may explain the differences observed with the distributions of laminar lengths for phases C–E. These departures in the distributions suggest that successive laminar lengths are decorrelated. This was checked by computing lengths for phases B versus the lengths for phases A, and so on for other successive phases. A cloud of points was obtained in each case thus revealing the lack of correlation between successive laminar phases in a sequence ABCDE. This would suggest that this is a multichannel intermittency [23] although the periodic orbits are always visited according to the sequence ABCDE. Thus, the dynamical behavior associated with the very particular time series shown in Fig. 6 is an intermittency structured around five periodic orbits which are visited always in the same order. The phase portrait made of the laminar phases and the bursts is quite complicated and may be viewed within the surface of a nontrivial torus. Such a torus must be embedded within a space with a dimension at least equal to 4 since heteroclinic connections involve one 2D unstable and one 2D stable manifolds. This is in agreement with the embed-
Fig. 16. Intensities associated with the unstable periodic orbits vs. the pumping parameter r. Other parameter values: δ = 3, σ = 0.01, γ = 0.2 and A = 0.05. Periodic orbits are detected by a threshold related to the amplitude of oscillations under which the intensity is assumed to be constant.
ding dimension we estimated to 4 using the algorithm written by Cao [24]. Such intermittent behavior can be observed over a quite large interval of the pumping parameter but the number of unstable periodic orbits involved varies. In order to have a global view of how this number depends on the pumping parameter, we plotted the intensities associated with the unstable periodic or-
Fig. 17. Two windows of a time series of the laser intensity with r = 37. Other parameter values as for Fig. 16.
D. Amroun et al. / Physica D 203 (2005) 185–197
bits versus r (Fig. 16). Roughly, the number of periodic orbits increases versus parameter r. Thus, there are three orbits for r ∈ [15; 18], four for r ∈ [18; 24], five for r ∈ [24; 35]. For greater values, up to six unstable periodic orbits are observed (r ∈ [35; 38]) and then decreases to five. For r > 41, a single stable periodic orbit remains because the lasing solution is stable again, according to the stability diagram shown in Fig. 1. The number of unstable periodic orbits can increase or decrease. It increases when a burst is splitted into two bursts when the pumping parameter is varied. Once these two bursts are sufficiently separated, a new plateau is created. Plateaux disappear by the reverse process. A time series where six plateau are clearly identified is shown in Fig. 17a for r = 37. We designate the sixth plateau by the letter B to indicate that it comes from the burst at the end of the plateau B. The chaotic character of the bursts between B and C implies very different lengths for plateau B . Thus, in a different window of the time series computed for r = 37, the plateau B is skipped since the two bursts are too close (Fig. 17b). This time series well exemplifies the way in which plateaux—or equivalently, the unstable periodic orbits—appear or disappear.
195
finite number of directions, and in a well-defined sequence. In the phase space reconstructed from the intensity where the symmetry associated with the phase of the electric field is modded out, it is shown that a single heteroclinic connection has a chaotic character, the others being associated with a unique path. This very uncommon behavior presents many characteristics of intermittency. In particular, the reinjection process near the periodic orbit visited after chaotic bursts is typical of type-II intermittency. It is expected that such a new type of dynamics could be observed in other spatio-temporal systems. Obviously, the number of periodic orbits around which the dynamics is structured is parameter-dependent but the global mechanism should be quite general. It seems that such spatio-temporal behaviors could be obtained for systems having traveling waves for solutions. Thus, perhaps the most appropriate experimental device for observing the complex intermittent behavior here reported could correspond to the semiconductor laser device where experimental traveling waves have been observed [27]. In the unstable domain of its parameter space, similar intermittent behaviors could be reproduced.
Acknowledgement 4. Conclusion Using the two-level Maxwell–Bloch equations with diffraction, complex dynamical behaviors of a homogeneously broadened single-mode laser with large Fresnel number have been investigated above the second laser threshold. When the cavity detuning δ is negative, there is a bi-angular laser emission in the farfield region. The underlying dynamics is quasi-periodic with two main frequencies. However, for positive cavity detuning δ, the dynamics is more complex. With the parameter values investigated here, a very particular evolution of the laser intensity has been observed. It has been shown that the dynamics is structured around five unstable periodic orbits. Heteroclinic connections between them are associated with bursts which interrupt phases during which the laser intensity is nearly constant. We showed that the order in which periodic orbits are visited is always the same. This is related to a sequence of five wave vectors. Thus, the laser spontaneously modifies its direction of emission in a
We wish to thank Robert Gilmore for helpful comments and encouraging this work.
Appendix A We describe in this section the numerical method used to solve the Maxwell–Bloch system (1). It is an adaptation of the Picard iteration method [25]. The latter is frequently used to prove the existence of solutions of the Cauchy problem for hyperbolic partial derivative equations (PDEs). In addition, the equation is split into a linear and a nonlinear part. In fact, only the nonlinear part is solved by means of Picard iterations, while the linear part is solved analytically by means of a Fourier transform with regard to the space variable. In the present situation, this modified scheme runs reasonably fast, and is very stable. A previous implementation of a scheme of this type is presented in [26]. We need to solve the Cauchy problem for the following system
196
D. Amroun et al. / Physica D 203 (2005) 185–197
of PDEs: 1 ∗ ∂d
∗ = −γ d + + e p) , (ep ∂t 2 ∂p = −(1 − iδ)p + e(d + r), ∂t 2 ∂e = −σ(e − p) + iA ∂ e , ∂t ∂x2
According to Picard’s method, we put Eq. (A.5) into the integral form: t
ˆ ˆ exp(−Lt )B(u(t )) dt , u(t) ˆ = exp(Lt) uˆ i + (A.8)
where we put d = d − r for sake of simplicity. This system can be split into the sum of a linear and a nonlinear part, as: ∂u = L(u) + B(u), (A.2) ∂t with −γd d −(1 − iδ)p + re u = p , L(u) = , ∂2 e e −σ(e − p) + iA 2 ∂x (A.3) γ − (ep∗ + e∗ p) 2 . B(u) = ed
(A.4)
0 Thus, L is a linear operator, while B is quadratic. We denote by u(k, ˆ t) the Fourier transform of u(x, t) with regard to the space variable x. We chose "x = 7 and δx = 0.1, where "x is the total length of the transverse dimension x and δx the sampling step. This fixes the maximum value of the wave vector to kmax = 10π. Then, the Fourier transform of Eq. (A.2) reads:
σ
where u(t ˆ = 0) = uˆ i is the initial data. An important point of the method is the analytical ˆ We computation of the exponential of the matrix Lt. find:
e−γ∗t
0 ˆ exp(Lt)= 0
0
0
e−Xt/2 Rt 2r sinh R 2 , −Xt/2 −Xt/2 Rt e e 2s sinh Q R R 2 (A.9) e−Xt/2 Q R
with X = 1 − iδ + ik2 A + σ, R= 2 2 A + σ)] X + 4[σr − (1 − iδ)(ik and Q = Rt R cosh Rt + (X − 2(1 − iδ)) sinh . Now we 2 2 have to solve Eq. (A.5) with the initial data uˆ i . We build, following the Picard iteration scheme, the recurrent sequence: • uˆ 0 is given (in practice, uˆ 0 (t) ≡ uˆ i ), • uˆ n known, uˆ n+1 is the solution of: ∂uˆ n+1 = A ˆ uˆ n+1 + B(u n ), ∂t uˆ n+1 (t = 0) = uˆ i .
(A.10)
It follows from Eq. (A.8) that: t
ˆ ˆ exp(−Lt )B(un (t )) dt . uˆ n+1 = exp(Lt) uˆ i + 0
∂uˆ ˆ uˆ + B(u), =L ∂t with −γ 0 ˆ = 0 −1 + iδ L 0
0
(A.1)
(A.5)
0 r −σ
,
(A.6)
− iAk2
γ ∗ p) − (e p∗ + e 2 . = B(u) ed 0
(A.7)
(A.11) Since at this stage B(u n ) is a known function, the computation of uˆ n+1 reduces to that of an integral. As soon as the sequence (u) ˆ n converges to a fixed point, its limit uˆ is the solution of Eq. (A.5). By computation of the inverse Fourier transform, we get u, the solution of system (A.1). In fact, the sequence uˆ n converges on some time interval [0, T ], where T is bounded from above. Thus, to calculate the solution for a longer range (a time series), we proceed as follows: we choose t0 < T and apply the
D. Amroun et al. / Physica D 203 (2005) 185–197
previous algorithm which gives u(t) in the range [0, t0 ]. Then, we apply it again with ui = u(t0 ) to obtain u in the interval [t0 , 2t0 ], and so on (as long as T is greater than t0 ). The main interest of this procedure is that the fixed point is reached with some chosen precision, on the whole interval [0, t0 ]. This way, we are able to control directly and quantitatively the difference between the numerical and the exact solution in the L∞ -norm.
References [1] H. Haken, Analogy between higher instabilities in fluids and lasers, Phys. Lett. A 53 (1) (1975) 77–78. [2] L.A. Lugiato, C. Oldano, L.M. Narducci, Cooperative frequency locking and stationary spatial structures in lasers, J. Opt. Soc. Am. B 5 (5) (1988) 879–888. [3] P. Coullet, L. Gil, F. Rocca, Optical vortices, Opt. Commun. 73 (5) (1989) 403–408. [4] N.B. Abraham, W.J. Firth, Special issue on transverse effects in optical systems, J. Opt. Soc. Am. 7 (1990) 951. [5] P.K. Jakobsen, J.V. Moloney, A.C. Newell, R. Indik, Space– time dynamics of wide-gain-section lasers, Phys. Rev. A 45 (1992) 8129–8137. [6] D. Dangoisse, D. Hennequin, C. Lepers, E. Louvergneaux, P. Glorieux, Two-dimensional optical lattices in a CO2 laser, Phys. Rev. A 46 (1992) 5955. [7] F.T. Arecchi, S. Boccaletti, P.L. Ramazza, S. Residori, Transition from boundary- to bulk-controlled regimes in optical pattern formation, Phys. Rev. Lett. 70 (1993) 2277. [8] M.C. Cross, P.H. Hohenberg, Pattern formation outside of equilibrium, Rev. Mod. Phys. 65 (1993) 851. [9] P.K. Jakobsen, J. Lega, Q. Feng, M. Staley, J.V. Moloney, A.C. Newell, Nonlinear transverse modes of large-aspect-ratio homogeneously broadened lasers. I. Analysis and numerical simulation, Phys. Rev. A 49 (1994) 4189. [10] J. Lega, P.K. Jakobsen, J.V. Moloney, A.C. Newell, Nonlinear transverse modes of large-aspect-ratio homogeneously broadened lasers. II. Pattern analysis near and beyond threshold, Phys. Rev. A 49 (1994) 4201.
197
[11] G. Huyet, M.C. Martinoni, J.R. Tredicce, S. Rica, Spatiotemporal dynamics of lasers with a large Fresnel number, Phys. Rev. Lett. 75 (1995) 4027–4030. [12] F. Encinas-Sanz, I. Leyva, J.M. Guerra, Time resolved pattern evolution in a large aperture laser, Phys. Rev. Lett. 84 (2000) 883–886. [13] Y.F. Chen, Y.P. Lan, Spontaneous transverse pattern formation in a microchip laser excited by a doughnut pump profile, Appl. Phys. B 75 (2002) 453–456. [14] D. Amroun, M. Brunel, H. Leblond, F. Sanchez, Space–time laser instabilities in homogeneously broadened dense media, J. Mod. Opt. 50 (2003) 1487–1495. [15] P. Ashwin, A.M. Rucklidge, Cycling chaos: its creation, persistence and loss of stability in a model of nonlinear magnetoconvection, Physica D 122 (1998) 134–154. [16] P. Ashwin, A.M. Rucklidge, R. Sturman, Infinities of stable periodic orbits in systems of coupled oscillators, Phys. Rev. E 66 (2002) 035201(R). [17] S.P. Dawson, A.M. Mancho, Collections of heteroclinic cycles in the Kuramoto–Sivashinsky equation, Physica D 100 (1997) 231–256. [18] N.H. Packard, J.P. Crutchfield, J.D. Farmer, R.S. Shaw, Geometry from a time series, Phys. Rev. Lett. 45 (9) (1980) 712–716. [19] C. Letellier, R. Gilmore, Covering dynamical systems: two-fold covers, Phys. Rev. E 63 (2001) 016206. [20] C. Letellier, Modding out a continuous rotation symmetry for disentangling a laser dynamics, Int. J. Bifurcation Chaos 13 (6) (2003) 1573–1577. [21] S. Wiggins, Introduction to Applied Nonlinear Dynamical Systems and Chaos, Springer-Verlag, 1990. [22] J.-Y. Huang, J.-J. Kim, Type-II intermittency in a coupled nonlinear oscillator: experimental observation, Phys. Rev. A 36 (1987) 1495–1497. [23] C. Letellier, P. Werny, J.-M. Malasoma, R. Gilmore, Multichannel intermittencies induced by symmetries, Phys. Rev. E 66 (2002) 036220. [24] L. Cao, Practical method for determining the minimum embedding dimension of a scalar time series, Physica D 110 (1997) 43–52. [25] F.B. Hildebrand, Introduction to Numerical Analysis, second ed., McGraw-Hill, New Delhi, 1974. [26] H. Leblond, Transverse stability of solitons and moving domain walls, J. Phys. A 33 (2000) 8105. [27] J. Houlihan, J.R. O’Callaghan, V. Voignier, G. Huyet, J.G. McInerney, B. Corbett, Experimental observation of traveling waves in a semiconductor laser device, Opt. Lett. 26 (2001) 1556–1558.