Dissipative structures in the resonant interaction of laser radiation with nonlinear dispersive medium

The article presents the results of studies on the stability of dissipative structures (DS) arising in the resonant interaction of laser radiation with a nonlinear medium. Resonant interaction is modeled by the one dimensional complex Ginzburg-Landau equation with a nonconservative cubic–quintic nonlinearity. The areas of existence of stable DS solutions have been determined analytically using a variational approach and confirmed numerically by extensive numerical simulations.


Introduction
The class of DS solutions draws a lot of attention due to promising applications in many scientific fields such as nonlinear optics, mode-locked lasers (Grelu and Akhmediev 2012), fibre lasers (Kalashnikov et al. 2006;Turitsyn et al. 2016). The CGLE represents a well known model for different types of dinamycal problems. It is used in a variety of physics fields like superconductivity (Ginzburg and Landau 1950), nonequilibrium phenomena (Aranson 2002), phase transitions (Riste 1975), binary fluid convection (Kolodner 1992), liquid crystals, Bose-Einstein condensates, nonlinear optics, but also in other sciences like biology, chemistry, etc (Aranson 2002;Nicolis and Prigogine 1977). The CGLE is similar to the nonlinear Schrodinger equation (NLSE), but contrary to NLSE which describes conservative systems, the CGLE can model dissipative systems and hence better describe "real-world" problems.
In order to obtain solitons in the NLSE, it is necessary to find the balance between dispersion and nonlinearity, while in addition to that solitons in the CGLE must balance between (linear or nonlinear) gain and energy loss. Because of that soliton solutions of the CGLE do not form a continuous family like NLSE, but discrete solutions for given set of parameters (Uvarova and Burenok 2016;Aleksić et al. 2012b;Skarka et al. 2014;Aleksić et al. 2015). Studies so far have revealed areas of stability for DS in CGLE with a nonzero conservative part of nonlinearity-the real part of nonlinear susceptibility (Kalashnikov et al. 2006). Under conditions of resonant interaction, the real part of the nonlinear susceptibility can be neglected (Fedorov et al. 2000). Model without conservative nonlinearity is presented in (Malomed et al. 2000) and the solitons are obtained using numerical methods.
In this work, we analyze analytically, using the variational method, and numerically, using finite difference and spectral methods, evolution and stability of one dimensional DS formed in the resonant interaction of laser radiation with nonlinear dispersive medium.

Model
From the Maxwell-Bloch equations, when the electromagnetic pulse duration is much greater than the relaxation time of the upper working level, in the corresponding dimensionless variables, the one dimensional CGLE for the electrical field has the form: where E is slowly varying amplitude of the envelope field, z is normalized length of propagation and t is normalized time. = ± 1 for anomalous and normal dispersion, respectively In two level atom model (Boyd 2008) the nonlinear functions N and G have the following forms : where ̄= , ̄= . The constants , are dimensionless detunings between the gain (amplification) and loss (absorption), with the respect to the spectral line center and the frequency of empty cavity mode. The coefficients and represent the linear amplification of active atoms and the absorption of passive atoms, (normalized to the absorption coefficient of the background medium), respectively (Fedorov et al. 2000). The linear part of the polarization determines the linear losses, temporal dispersion of the media, and spectral filtering ( ). Under conditions of resonant interactions, when → 0 and → 0 , term N(|E| 2 ) in Eq. (1) can be neglected.
In order to exploit variational methods which have been studied and well developed for CGLE with a polynomial nonlinearity, we approximate the saturable nonlinearity given by Eq.
When coefficients and satisfy the following conditions function G(I) has two roots: for As shown in paper (Aleksić et al. 2020) under these conditions for function G, it is possible to obtain a good approximation using a second order polynomial where , , , and represent the linear loss, nonlinear loss, nonlinear gain, and coefficient of spectral filtering, respectively.

Variational approach
Approximate solution of Eq. (11) is obtained using the variational method extended to dissipative systems (Skarka and Aleksić 2006), which provides a systematic way of finding approximate analytical solutions. The trial function for the variational Ritz-Kantorovich method can be written as: phase, which all depend on the coordinate z.
After optimizing the functions obtained from the standard procedure (Skarka and Aleksić 2006), we get corresponding differential Euler-Lagrange equations: The stationary point of the system of Eqs. (20)-(22) corresponds to the stationary solution of (11), and is determined by the system of equations where e = ∕ √ 1 + 2 . Eq. (22), expresses the energy balance.
Here A − and A +1,2 correspond to C < 0 and C > 0 , respectively. The bifurcation point is defined by the following expressions: = (4 − e ) 2 2 ∕(8(3 − e ) * ) and . Physically, the essential characteristic of these solutions is that the dissipative nonlinearity creates a nonzero parameter of chirp C , which is in balance with dispersion. The balance is achieved regardless of whether the dispersion is normal or anomalous. There are no solutions when C = 0 . When C → 0 then R → 0 and the solution approaches homogeneous distribution.

Stability analysis
The stationary point (20)-(22) is stable according to Lyapunov criteria if all the eigenvalues of the Jacobian matrix F i ∕ j , i = A, T, C have negative real parts. In general the equation det F i ∕ j − I = 0 is a polynomial of degree 3: where and F i j = F i ∕ j . According to the Hurwitz criterion for the stability of a stationary point, the necessary and sufficient conditions to have Re( ) < 0 are: The stability region, where all the coefficients are positive, is bounded by the following relations: and corresponds to the solution A +2 where C > 0 . Thus, analytically stable fundamental solutions in the case of normal dispersion have a negative chirp parameter C < 0 and positive C > 0 in the case of anomalous dispersion.
Note that Eq. (11) has two homogeneous nontrivial solutions where |E| + is stable and |E| − unstable. Figure 1 shows the analytical values for the amplitudes in unstable cases A − (dashed black) and A +1 (dashed orange), and in the stable case A +2 (red) in terms of . Markers ▪ − e = 0.0 , ▴ − e = 0.2 , ⧫ − e = 0.5 represent the final value of the amplitude of stable localized DS obtained in the numerical simulations. One can see good agreement between predicted and confirmed amplitude values as the markers are almost perfectly aligned with the corresponding red curves. However, not all points on the red line remain stable and for bigger the initial structures are attracted to the homogeneous solution Eq. (34) which is represented by the green line A at = |E| + . Marker • corresponds to the amplitude of the homogeneous solutions obtained numerically, regardless of the value of e . (33)

Numerical Simulations
The stability regions obtained from the variational approach are approximate due to the fact that the solution is approximate. The conclusion about stability was made on the basis of the analysis of the finite number of ordinary differential equations (ODE), however the CGLE is a partial differential equation, which corresponds to a system consisted of infinitely large number of ODEs. Exhaustive numerical simulations are necessary to confirm the stability of the obtained analytical solutions of the CGLE. The simulations are performed using modified finite difference time domain (FDTD) method (Aleksić et al. 2012), as well as spectral split-step method (Aleksić and LA 2018). To make sure that the predicted DS solutions remain stable, every numerical simulation has been propagated for more than 100 dispersion lengths using the values for A, T and C given by the stationary point. Without loss of generality, we use fixed parameters, = 1 and * = 1 , in the all numerical simulations. Figure 2 shows results of the numerical analysis of the analytically predicted stability region. The analytically obtained region is the area between the two black lines and corresponds to the relation given in (33). Red dots represent the stable solutions where the parameters of the initial pulse change slightly and quickly reach their stable value. Figure 3 shows the evolution of the properties: (a) amplitude A , (b) pulse width T and power P = ∫ |E| 2 dt , and (c) initial and final field distribution, of the solution obtained as results of the numerical simulation with parameters e = 0.26 and = 0.6 (located in the red area of Fig. 2). Starting from initial values in the parameter space, which are obtained analytically, each property exhibits small variations before it converges to a steady state value, which remains constant until the end of the simulation. Figure 3(d), which represents dependence of the field distribution |E| on z, shows almost imperceptible changes in the shape and properties of the solution confirming the accuracy of the used analytical method.
Green dots in Fig. 2 represent unstable solutions where the initial pulse turns into a stable homogeneous distribution given by Eq. (34). In the green area the homogeneous distribution acts like an attractor. Figure 4 shows the evolution of the same properties like in fig. 3 just in case of an unstable solution. The result are obtained in the numerical simulation with parameters e = 0.74 and = 0.55 (located in the green area of fig. 2). The amplitude of the Yellow dots in Fig. 2 represent unstable region of the solutions that rapidly lose power and disappear.
Blue dots in Fig. 2 lie between stable (red dots) and unstable (green dots) regions, and represent unstable solutions which neither disappear (yellow dots) nor turn into homogeneous solution (green dots), but take some other shape.
It has been observed that variety of exotic structures can occur at the boundary between stable and unstable regions in the case of two dimensional CGLE with conservative nonlinearity (Soto-Crespo et al. 2008). Similar to that, at some points located at the boundary between red and blue region in Fig. 2, DS can change their initial Gaussian shape and turn into a stable non-Gaussian structure. Figure 5 shows initial and final distribution of two examples of numerically obtained stable structures. In case of lower e we observe stable DS in a form of a bell Fig. 5(a). Higher e results in a flat-top solution. To confirm stability of the observed structures the numerical simulation is performed until z = 10 4 .
Besides the structures shown in Fig. 5, at the boundary we can also observe a breather structure like the one in Fig. 6. The parameters (T and P) of the obtained breather experience oscillations with constant amplitude and frequency during the evolution Fig. 6(a) which confirms the long-life existence of the breather. Figure 6(b) shows contour plot of the breather's field distribution as a function of z.

Conclusion
In this work, the stability of DS in the resonant interaction of laser radiation with nonlinear dispersive medium is analyzed analytically and numerically. As a result of the variational approach, the domain of existence of DS is analytically obtained. Analytical solutions have been tested by exhaustive numerical simulations which have resulted in defining the area of stability and confirmed the existence of stable solutions for both normal and anomalous dispersion. "Exotic" solutions have been observed at the boundary between stable and