Received 4 May 2016; accepted 10 June 2016; published 13 June 2016
Dissipative systems in nonlinear optics have a dynamic involving many different phenomena such as nonlinear gain, the saturable losses, the dispersion and others effects. The interaction between these different physical manifestations leads to a rich variety of structures  . Dissipative solitons are stable elementary localized structures that appear in dissipative nonlinear medium and depend on the balance between nonlinearity and dispersion/diffraction, and the balance between gain and loss. They have several properties such as their internal energy exchange mechanism which is similar to that of biological structures or the possibility to reveal pulsating dynamics  . These properties make them attractive for research. The concept of dissipative solitons has been widely studied in several fields of nonlinear science  such as nonlinear dissipative optics. Among the important applications there are passively mode-locked laser systems and optical transmission lines.
For dissipative systems in nonlinear optics, stable solitons can arise in one, two, and three dimensions  . For example, in case of the one-dimensional (1D), new types of localized waves  -  have been experimentally  observed in lasers cavities.
With a good choice of the parameters of the system, stable solitons over very large distances of propagation are perfectly obtained. However, a variation of the parameters of the system changes the types of solutions via bifurcations. Dissipative soliton stability depends crucially on the energy balance and exists as long as there is a continuous energy supply to the system. Its shape, amplitude, and velocity are all fixed and defined by the parameters of the system  rather than by the initial condition.
Properties and conditions of the existence of solitons have been extensively studied in (1D) one-dimension and (2D) two-dimension. The (3D) three-dimensional case is still largely in its infancy and requires an extremely lengthy and costly procedure. Indeed solving numerically a (3D) equation for a given set of parameters and a given initial condition can take up to several days in a standard PC. In this context, it is important to develop theoretical tools that can perceive soliton solutions more efficiently and envisage their domains of existence.
It has been recently demonstrated in previous works that the collective variable approach is a useful tool which can reduce significantly the computation time. This approach is useful for predicting approximately the domains of existence of stable light bullets in the parameter space of the (3D) complex cubic-quintic Ginzburg- Landau equation  .
Here, previous studies are focusing on this type of (3D) optical pulse in order to obtain the properties and conditions of their existence. According to the specific choice of the trial function, different internal dynamics of the pulsating dissipative soliton are shown.
2. Complex Cubic-Quintic Ginzburg-Landau Equation (CGLE)
In this study, the propagation of pulse in a system described by an extended complex Cubic-quintic Ginzburg- Landau equation model is considered. This equation (CGLE) is one of the universal equations used to describe dissipative systems. Many nonequilibrium phenomena, such as the generation of spatio-temporal dissipative structure in lasers  , soliton propagation in optical fiber systems with linear and nonlinear gain and spectral filtering (such as communication links with lumped fast saturable absorbers or fiber lasers with additive-pulse mode-locking or nonlinear polarization rotation) may all be described by the CGLE  . The quintic dissipative term in CGLE is essential to provide the stability of the optical pulse  . Moreover this equation could be applied to the modeling of a wide-aperture laser cavity in the short pulse regime of operation. The model includes the effects of two-dimensional transverse diffraction of the beam, longitudinal dispersion of the pulse, and its evolution along the cavity. Dissipative terms describe the gain and loss of the pulse in the cavity. The propagation equation reads:
Equation (1) is written in normalized form. The optical envelope ψ is a complex function of four real variables, where t is the retarded time in the frame moving with the pulse, z is the propagation distance, and x and y are the two transverse coordinates.
The left-hand-side contains the conservative terms, namely which is for the anomalous (normal) dispersion propagation regime and ν which represents, if negative, the saturation coefficient of the Kerr nonlinearity. In the following, the dispersion is anomalous, and ν is kept relatively small. The right-hand-side of Equation (1) includes all dissipative terms: δ, ε, β and μ are the coefficients for linear loss (if negative), nonlinear gain (if positive), spectral filtering (if positive) and saturation of the nonlinear gain (if negative), respectively. This distributed equation finds application in modeling for instance a wide-aperture active optical cavity in the regime of short pulse operation. The model includes the effects of two-dimensional transverse diffraction of the beam, longitudinal dispersion of the pulse and its evolution along the cavity. Dissipative terms describe gain and loss of the pulse in the cavity. Higher-order dissipative terms are responsible for the nonlinear transmission characteristics of the cavity which allows, for example, passive mode-locking. This equation is a natural extension of the one-dimensional CGLE.
Now, to best of our Knowledge, there is no analytical solution for the (3D) complex Cubic-quintic Ginzburg- Landau equation. Indeed, research conducted to date use namely direct numerical simulations of the CGLE. This procedure is extremely and costly  .
The lack of general analytical solutions for the complex Cubic-quintic Ginzburg-Landau equation CGLE leads us to the necessity of using simple approaches to explore the existence of certain class of solutions. A few approximate, semi-analytical methods based on various physical backgrounds were developed and applied to study nonlinear pulse propagation. Several of them make use of a trial function and its associated finite-dimen- sional dynamical system.
Our study is therefore to develop theoretical tools that can perceive soliton solutions more efficiently and envisage their domains of existence. Using collective variable approach allow to rich variety of solutions that include stationary and pulsating dissipative solitons.
3. Collective Variable Approach
The dynamics of ultra-short light pulses in a fiber system can be described by the Equation (1). But the field describes not only the pulse as a collective entity (localized in time and space) but also all other localized or non-localized excitations, such as noise or radiation, which are always more or less present in the real system. The pulse may not only be able to translate as a whole entity but it may also execute more or less complex internal vibrations depending on the type of the perturbations in the system. In this context it is useful, then, to somehow simplify the characterization of the pulse by use of a low dimensional equivalent mechanical system based on a finite number of degrees of freedom  . Each degree of freedom can then be described by means of a coordinate called the collective variable. The mean idea in the collective variable approach is to associate collective variables with the pulse’s parameters of interest for which equations of motion may be derived. One may introduce N collective variables, z dependent, say Xi with, in a way such that each of them can correctly describe a fundamental parameter of the pulse (amplitude, width, chirp, …)  . To this end, one can decompose the field in the following way:
where f, the ansatz function or trial function is a function of the collective variables. And is a residual field that represents all other excitations in the system (noise, radiation, dressing field, etc.)  . The precise form of the ansatz function that introduces the collective variables in the theory is rather crucial, especially when approximations are made.
The choice of the trial function is important for the success of the technique. It is impossible to get any idea of the exact profile without solving numerically the CGLE (1). Hence, the choice of a particular trial function has a certain degree of arbitrariness, stressing furthermore the approximate nature of the collective variable approach. After choosing the ansatz function one can pursue the process of characterization of the pulse by neglecting the residual field. This approximation is called the bare approximation  . In this way one can consider the fact that the pulse propagation can be completely characterized that the ansatz function. This approximation depends very strongly on the choice of the trial function. Consequently, when approximations are made, such as the use of the variational approach, the precise shape of the trial function that introduces the collective variables in the analysis of the dynamics of soliton becomes rather crucial.
4. Variational Equations
In this study the propagation of solitons in a system described by an extended (3D) complex Cubic-quintic Ginzburg-Landau equation model is considered. This model includes cubic and quintic nonlinearities of dispersive and dissipative types. Transverse operators are added to take account of the spatial diffraction in the wave paraxial approximation.
To obtain a better understanding of dynamic processes which influence the behavior of pulse during its propagation, it is considered that the real field represents the ansatz function which is endowed with the usual physical meaning of the pulse parameters.
Thus, it is possible to reduce the equation of impulsive field to an ordinary differential equation (ODE) describing the evolution of the parameters of the soliton during propagation. The advantage is that, the ordinary differential equation can be solved numerically with relative ease.
However, the precise shape of the trial function is crucial to have solutions with the desired properties. For the complex systems, this choice is preceded by careful analysis and based on a comparison of some analytical approximate solutions before an exact numerical solving of the CGLE can be made for comparison.
In previous work  , it has been demonstrated with different trial functions that the collective variable approach is a useful tool for predicting approximately the domains of existence of stable light bullets in the parameter space of the complex cubic-quintic. These predictions were confirmed by the direct numerical solutions of the CGLE, qualitatively as well as quantitatively, that give us confidence in our collective variable approach.
On the basis of preliminary results of previous works  and to describe a richer variety of dynamical behaviors, a Gaussian trial function which admits asymmetric pulse shapes in the transverse plane is used and is given by
A, wt, wx, wy, ct, cx, cy, and p represent the collective variables. With t, x and y, the temporal and transverses variables along x and y axis respectively. A stands for soliton amplitude, wt the temporal width along t, wx the transverse width along x axis and wy the transverse width along y axis. ct represents the temporal chirp parameter, cx the transverse chirp along x, cy the transverse chirp along y and p is the global phase that evolves along with propagation. When a stationary regime is reached, the phase becomes a linear function of the propagation distance z.
The choice of the trial function is done according to the master Equation (1) and the type of solutions pursued. Consequently, the precise shape of the trial function that introduces the collective variable in the analysis of the dynamics of soliton becomes rather crucial. After choosing the trial function one can pursue the process of characterization of the pulse in two completely different ways depending on the level of accuracy desired. First, one can make use of the exact pulse field to obtain the pulse parameters. The second approach of characterization would be to carry out a variational analysis neglecting the residual field. The approximation of neglecting the residual field is called the bare approximation  , just like those in previous works  . This approach leads generally to results which are less accurate than the first one, but its application to the characterization of the propagation of soliton would be too times consuming.
After choosing the trial function one can carry out variational analysis neglecting the residual field (the bare approximation). Applying the bare approximation to the 3D CGLE, that is, substituting the field by the given trial function and projecting the resulting equations in the following direction
The collective variables evolve according to the following set of eight coupled ordinary differential equation are easily obtained:
Obviously the ordinary differential equations in Equation (4) depend on the choice the Equation (3).
These equations give no explicit information with regard to the different solutions of the Equation (1) and their stability. They give us the first idea on the dynamic of the light pulse. The variational equations allow seeing the influence of each Equation (1) parameters on the various physical parameters of the soliton. They are usually functions of time that evolve subject to the constraints of the system and finally converge to fixed point or a limit cycle.
A meticulous analysis of the variational equations show that the evolution of the amplitude (A) is dominated by the linear loss (δ), the nonlinear gain (ε) and its saturation (μ), as well as that the terms of spectral filtering (β) and dispersion term (D). This confirms quite well that the perfect balance between losses and gains is required to maintain the shape and stability of the soliton. The temporal (wt) and spatial widths (wx, wy) also depend on the nonlinear gain (ε) and its saturation (μ). As expected, the terms of spectral filtering (β) and dispersion term (D) affect the temporal width. As well the spatial parameters cx, cy and temporal parameters ct are influenced in the same way by the Kerr term saturation of the optical nonlinearity (ν), but the temporal term is also affected by the terms of spectral filtering (β) and dispersion term (D). Finally, not any parameters of the soliton are influenced by (p), the global phase.
One advantage of the collective variable approach is that the total energy can also expressed as function of the trial function parameters. Here it is interesting to gain insight from this simple and useful quantity, which is defined as
This above equation point out that the total energy is strongly ruled by the amplitude (A) the temporal (wt) and spatial widths (wx, wy).
Hence, the significance of our collective variable approach helps to analyze the variational equations and the influence of various parameters. Following this in-depth analysis, the first major step in our study is to provide a mapping of the regions of existence of stable solutions in the parameter space of the (3D) CGLE.
5. (3D) Pulsating Dissipative Soliton
From the analytical results of the variational Equation (4), the fixed points are carefully analyzed and their stability studied. The fixed points (FPs) of the system are found by imposing the left-hand side of Equation (4) to be zero. The threshold of existence of FPs can be estimated by the relation  . If, both stable and unstable fixed points appear. The stable fixed points correspond to stationary solutions. The unstable fixed points correspond to the pulsating and non-stationary solutions.
A major goal of our study is to provide a quick approximate mapping of the regions of existence of stable and unstable solutions in the parameter space of the Equation (1). Stationary and pulsating solutions exist in regions defined by the space of the parameters of the equation, but it is extremely difficult to map these regions by varying all these parameters. To remedy this, it is convenient to set all parameters of the equation except two chosen as variable parameters. Thus the different solutions obtained are defined in this plan. Chosen variable parameters are the dispersion D and the cubic nonlinear gain ε. So for each set of parameters, the existence of the fixed point and its stability is studied.
By investigating the parameter regions situated in the neighborhood of the parameters μ = −0.1, δ = −0.4, β = 0.1, ν = −0.08 and γ = 1 and according to our previous study   , a rich variety of solutions of Equation (1) including stationary and pulsating dissipative solitons is observed in the plane (D, ε).
For a given set of CGLE parameters (μ, δ, β, ν, et γ), with an initial pulse
For each value pair (D, ε), the fourth-order Runge-Kutta algorithm provides the fixed point whose stability is analyzed.
The mapping Figure 1 brings to light this rigorous analysis. The dispersion coefficient D and the non-linear gain coefficient ε are the parameters that can be modified easily in most experiments. It is the case of experiments using lasers passively locked mode.
Figure 1. Cartography of stationary (color domain) and pulsating (dotted line domain) dissipative soliton found from collective variable approach in the (D, ε) plane using the Gaussian trial function. Other CGLE parameters appear inside the figure.
Regions of existence of stable light bullets following the same technique that we used for dissipative solitons in the (3 + 1)D dimensional case  have been determined.
Inside this cartography, two different regions according to the values of the dispersion and the nonlinear gain are observed. The coloured area contains the stable fixed points, which are the basin of shows. And near this points, all the others points converge. The stable fixed points regions represents the domain of stationary solitons of the Equation (1) found from collective variable approach. In this domain, all the solitons parameters (amplitude, width, chirp…) stay stationary throughout propagation. Above the stationary domain, instable fixed points which can be dived in two categories: the limit-cycle attractor and the instable solutions could be observed. Our main interest is to study the dynamic of the pulse in the limit-cycle attractor area (dotted line domain). Indeed close to boundaries of the existence domains of stationary solitons solutions, there is more often an intermediate region in which pulsating solutions can be found. The existence of pulsating solutions is indeed a general feature of most nonlinear dissipative systems. This behavior of pulsating soliton can be attributed of a limit-cycle attractor; it then possesses inherent stability the same way as stationary stable solutions do. Above the pulsating domain, instable solitons which are little physical interest in this current study are noticed. The cartography (Figure 1) is quite like those achieved in others parameters space (v, ε), (μ, ε), (β, ε) and (μ, v)   . Thus according the value of the cubic nonlinear gain ε, the stability diagram shows types of dynamic. In the two-dimensional plane of parameters, the regions of existence of stable stationary solitons presented above are contained between an upper and a lower values of ε. Below the lower limit, the energy pumped into the system is not enough to support the solitons. In this case the localized solution dissipates and eventually vanishes. On the contrary, when is above the upper value, the energy supply inflates the soliton to the extent that it grows indefinitely in size. However, this process does not occur in one single step. There is a small intermediate region of ε values where stationary solutions are transformed into pulsating ones before the continuous inflating begins at even higher ε. On can notice that in the pulsating domain, for each a paired value of (D, ε), there is one and only one type of pulsating soliton with its own characteristics (amplitude, widths, energy, …). From Figure 2 and Figure 3, pulsating solutions in the form of stable limit cycles are observed. These solutions can be considered as attractors. These figures are obtained thanks to our collective variable approach for μ = −0.1, δ = −0.4, ε = 0.59, β = 0.1, D = 6, ν = −0.08, and γ = 1.
Figure 2 (below the enlarged views) shows the evolution on the total energy of the soliton for a given set of parameters corresponding in the pulsating domain (dotted line domain in Figure 1). Dynamics begins with a short transitional phase, which is followed by a permanent dynamics with an increase of the energy are noticed. The energy remains in this stable state when the permanent dynamic is reached. The total energy gives us the main information about the soliton dynamics. It’s not conserved but evolves in accordance with the so-called balance equation. When a stationary solution is reached, the total energy converges to a constant value. However,
Figure 2. Evolution of the total pulse energy, enlarged view of this pulsation is plotted bottom with μ = −0.1, ε = 0.59, D = 6, δ = −0.4, β = 0.1, and γ = 1.
Figure 3. Radially asymmetric pulsations obtained with the use of the collective variable approach and the trial function. The evolution of the following collective variables with respect to z is shown: (c1) the amplitude, (c2) the temporal width, (c3) the width along the x axis, and (c4) the width along the y axis. The onset of stable harmonic pulsations appears clearly. Enlarged views of those pulsations are plotted to the right in (d1)-(d4). The pulsations along the x and y axes are out of phase and share the same amplitude μ = −0.1, ε = 0.59, D = 6, δ = −0.4, β = 0.1, and γ = 1.
the soliton, is a pulsating one, the total energy is an oscillating function of z. In Figure 3, all the others parameters (amplitude, widths, …) of the soliton follow the same dynamics. The dynamics are illustrated in this picture [d(i) are the enlarged views of the c(i)]. It shows the evolution of the soliton amplitude (Figure 3(c1)), temporal (Figure 3(c2)) width and the transverse width along x (Figure 3(c3)), and y (Figure 3(c4)), axis. After a brief transitional period characterized by small oscillations due to the adjustment of the initial condition, and after the onset of oscillations, the pulsating dynamics becomes steady. It is characterized by large oscillations between the two constant limits. The pulse breathing in the transverse and in the temporal domain can be seen from the evolution plots Figure 3(d3) and Figure 3(d4). The amplitude oscillates between 2.82 and 2.92, the temporal width between 2.78 and 2.85 and the transverse widths between 0.9 and 1.5. It could be clearly seen that the difference oscillation between the temporal and the transverse width. As well, on the same propagation distance z, the amplitude of the soliton evolves much faster compared with its spatial widths. The same holds true for the temporal width. The close-up view of oscillations (Figure 3(d3) and Figure 3(d4)) shows a nearly harmonic evolution of the collective variables. The x and y oscillations are out of phase and have the same magnitude. This shows periodic out of phase consecutive contractions of the soliton in the x and y directions. At the same time, the total energy is not a constant and the amplitude and temporal width oscillate at a doubled frequency.
5.1. Bifurcation Diagram
In Figure 1, two different kind of stable 3D solutions in the specific domain appear. When a pulsating solution can be attributed to the existence of a limit-cycle attractor, it then possesses inherent stability the same way as stationary stable solutions do. For a given set of the parameters and varying the nonlinear gain parameter, the soliton move from a stationary solution to pulsating one. So the total energy leads to a pulsation mode whose spectrum contains two (three) main frequencies. To show pulsations with two oscillation periods, the system has to undergo two bifurcations from the stationary solution. And for the three oscillation periods, the system has to undergo four bifurcations from the two bifurcations. The picture Figure 4 shows a bifurcation diagram obtained when varying the nonlinear gain parameter ε from 0.575 to 0.595, while keeping the rest of parameters (written in Figure 4) fixed. The curve represents a local maximum or minimum of the total energy. For the nonlinear gain value inferior to 0.583 the energy has a single value (maximum equal to minimum). This means that stationary solution is reach. From the up value to 0.583, a bifurcation occurs, which is related to the onset of pulsations with single frequency. In the case discussed above, energy pulsations remain very close to harmonic pulsation.
5.2. The Effect of the Saturation of the Kerr Nonlinearity
The influence of the saturation coefficient of the Kerr nonlinearity ν on the oscillations is investigated. To this end, parameters a μ = −0.1, ε = 0.6, D = 6, δ = −0.4, β = 0.1, γ = 1 are set with different values of ν (−0.07, −0.08, −0.09). The dynamics of the energy with respect to the propagation distance z are followed. Figure 5 stands for the result of this detailed analysis [b(i) are the enlarged views of the a(i)]. Obviously, the value of saturation of kerr nonlinearity for which the soliton oscillates is chosen, but this dynamics changes when ν is modified. So for different values of the saturation nonlinear parameter ν, the other parameters (written below the picture Figure 5) are fixed. When ν decrease (from −0.07 to −0.09), the amplitude of the pulsations increase, e.g. the energy augment. In Figure 5(b1) and Figure 5(b2) the same dynamics with stable oscillation is obtained, but in Figure 5(b3) the motion changes. More precisely, pulsations that consist of two oscillation periods instead of one are found. With precision, the asymptotic evolution of the total energy leads to a pulsation whose spectrum contains two main frequencies, associated to intensities (Figure 5(b3)) of the same magnitude. It could also be noticed that the internal dynamic change including many oscillations in the same period. It emerges from this study that according to the value of saturation nonlinearity and to the choice of the trial function, it is possible to find the other types of pulsating dissipative soliton and the dynamic within the same period differs, from harmonic to non-harmonic pulsations.
Figure 4. Bifurcation diagram of the dissipative pulsating soliton.
Figure 5. Evolution of the total pulse energy for different values of the Kerr nonlinearity saturation, with μ = −0.1, ε = 0.6, D = 6, δ = −0.4, β = 0.1 and γ = 1. Enlarged views of those pulsations are plotted to the right in (b1)-(b3). The pulsation (b3) contains two main frequencies.
In this work, based on collective variable approach, we have expanded the studies of 3D dissipative pulsating solitons in the complex Ginzburg-Landau equation with the cubic-quintic nonlinearity. In particular, the regions of coexistence of stationary and pulsating dissipative soliton are obtained. The collective variable approach is very efficient for approximating stable pulsating solutions when a suitable trial function is chosen. A nearly harmonic evolution of the widths along x and y axis is shown and the x and y oscillations are out of phase with the same magnitude. The dynamics of soliton can be controlled by the choice of the system parameters. So according to the values of the nonlinear gain, the system has to undergo a bifurcation from the stationary solution, to obtain pulsations with two oscillation periods. A pulsating soliton whose spectrum contains two main frequencies (two oscillation periods instead of one), associated to intensities of the same magnitude could also be predicted. The latter effect may be used, in principle, to grow photonic channels and multichannel arrays in bulk optical media. The collective variable technique is incomparably quicker than direct numerical computations. Of course, it should be used at the final stage of studies to confirm, complement, or invalidate the collective variable approach predictions.
 Saarlos, W.V. and Hhenberg, P.C. (1992) Fronts, Pulse, Sources and Sinks in Generalized Ginzburg-Landau. Physica D, 56, 303-367.
 Soto-Crespo, J.M., Akhmediev, N. and Town, G. (2002) Continuous-Wave versus Pulse Regime in Passively Mode Locked Laser with a Fast Saturable Absorber. Journal of the Optical Society of America B, 19, 234-242.
 Akhmediev, N. and Ankiewicz, A. (2005) Dissipative Solitons. Springer, Heidelberg.
 Rosanov, N.N. (2002) Spatial Hysteresis and Optical Patterns. Springer, Berlin.
 Sakaguchi, H. and Malomed, B.A. (2001) Instabilities and Splitting of Pulses in Coupled Ginzburg-Landau Equations. Physica D, 154, 229-239.
 Akhmediev, N., Soto-Crespo, J.M. and Town, G. (2001) Pulsating Solitons, Chaotic Solitons, Period Doubling, and Pulse Coexistence in Mode-Locked Lasers: CGLE Approach. Physical Review E, 63, Article ID: 056602.
 Deissler, R.J. and Brand, H.R. (1994) Periodic, Quasiperiodic, and Chaotic Localized Solutions of the Quintic Complex Ginzburg-Landau Equation. Physical Review Letters, 72, 478-481.
 Soto-Crespo, J.M., Grelu, M.Ph. and Akhmediev, N. (2004) Bifurcations and Multiple-Period Soliton Pulsations in a Passively Mode-Locked Fiber Laser. Physical Review E, 70, Article ID: 066612.
 Akhmediev, N. and Ankiewicz, A. (2003) Solitons around Us: Integrable, Hamiltonian and Dissipative Systems. In: Porsezian, K. and Kurakose, V.C., Eds., Optical Solitons: Theoretical and Experimental Challenges, Springer, Heidelberg.
 Kamagaté, A., Grelu, Ph., Tchofo-Dinda, P., Soto-Crespo, J.M. and Akhmediev, N. (2009) Stationary and Pulsating Dissipative Light Bullets from a Collective Variable Approach. Physical Review E, 79, Article ID: 026609.
 Haus, H.A., Fujimoto, J.G. and Ippen, E.P. (1995) Structures for Additive Pulse Mode Locking. Journal of the Optical Society of America B, 8, 2068-2076.
 Soto-Crespo, J.M., Akhmediev, N.N., Afanasjev, V.V. and Wabnitz, S. (1997) Pulse Solutions of the Cubic-Quintic Complex Ginzburg-Landau Equation in the Case of Normal Dispersion. Physical Review E, 55, 4783-4796.
 Moores, J.D. (1993) On the Ginzburg-Landau Laser Mode-Locking Model with Fifth-Odersaturable Absorber Term. Optics Communications, 96, 65-70.
 Soto-Crespo, J.M., Akhmediev, N. and Grelu, P. (2006) Optical Bullets and Double Bullet Complexes in Dissipative Systems. Physical Review E, 74, Article ID: 046612.
 Boesch, R., Stancioff, P. and Willis, C.R. (1988) Hamiltonian Equations for Multiple-Collective-Variable Theories of Nonlinear Klein-Gordon Equations: A Projection-Operator Approach. Physical Review B, 38, 6713-6735.
 Tchofo-Dinda, P., Moubissi, A.B. and Nakkeeran, K. (2001) Collective Variable Theory for Optical Solitons in Fibers. Physical Review E, 64, Article ID: 016608.
 Soto-Crespo, J.M., Akhmediev, N. and Ankiewicz, A. (2000) Pulsating, Creeping, and Erupting Solitons in Dissipative Systems. Physical Review Letters, 85, 2937.