Regular, beating and dilogarithmic breathers in biased photorefractive crystals

The propagation of light beams in photovoltaic pyroelectric photorefractive crystals is modelled by a specific generalization of the nonlinear Schrödinger equation (GNLSE). We use the variational approximation (VA) to predict the propagation of solitary-wave inputs in the crystal, finding that the VA equations involve the dilogarithm special function. The VA predicts that solitons and breathers exist, and the Vakhitov-Kolokolov criterion predicts that the solitons are stable solutions. Direct simulations of the underlying GNLSE corroborates the existence of such stable modes. The numerical solutions produce both regular breathers and ones featuring beats (long-period modulations of fast oscillations). In the latter case, the Fourier transform of amplitude oscillations reveals a nearly discrete spectrum characterizing the beats dynamics. Numerical solutions of another type demonstrate spontaneous splitting of the input pulse in two or several secondary ones.

The propagation of breathers in photorefractive crystals has received less attention, in comparison to breathers which appear in other systems [12][13][14][15][16]).In the present work, we focus on the breather dynamics.
The propagation of breathers and solitons in photorefractive crystals, apart from their possible applications [17], is particularly interesting from a theoretical point of view, as it is modelled by equations similar to the nonlinear Schrödinger (NLS) equation.In the present work we study the propagation of light beams in a biased photovoltaic photorefractive crystal featuring the pyroelectric effect.It is known that this medium maintains spatial solitons, which are governed by a generalized NLS equation (GNLSE), that includes the usual diffractive term (the second derivative with respect to the transversal coordinate), Kerr-type and saturable nonlinearities, produced by the pyroelectric and photovoltaic effects, respectively, and a nonlinear term induced by the external field [18].In this work we introduce a saturation parameter, which is a measure of the combined effect of the three distinct nonlinearities.The corresponding GNLSE is not integrable, unlike the classical NLS equation [19,20], but it admits a Lagrangian structure.The latter property suggests constructing solutions by means of the Anderson´s variational approximation (VA) [21,22], which is used, in various forms, in studies of optical solitons [23][24][25][26][27][28][29][30][31], see also reviews in Refs.[32] and [33].The existence of the Lagrangian makes it also possible to apply the Noether´s theorem to identify dynamical invariants of GNLSE [33].
The rest of the paper is structured as follows.In Sec. 2 we derive the GNLSE modelling the propagation of light in the setting under the consideration.In Sec. 3 we introduce the Lagrangian for the GNLSE and elaborate the VA, using two different ansätze (trial forms of the approximate solution).We thus derive two sets of Euler-Lagrange equations corresponding to these ansätze.They involve effective Lagrangians which include the dilogarithm function Li2(z), where z is the propagation distance (evolution variable).With the help of Noether´s theorem, we obtain dynamical invariants (conserved quantities) of the GNLSE.Moreover, the Euler-Lagrange equations obtained via the VA permit us to apply the Vakhitov-Kolokolov criterion in order to show that the soliton solutions of this equation are stable.In Sec. 4 we demonstrate that the VA produces both, stationary solitons and breathers.The existence of the breathers is confirmed by direct numerical simulations of the GNLSE.The simulations also reveal robust breathers with intrinsic vibrations featuring beats, as well as solutions exhibiting spontaneous splitting of the initial pulse in two or several secondary ones.The paper is concluded by Sec. 5.

The fundamental equation (GNLSE)
We consider a non-centrosymmetric photorefractive crystal, whose principal axes coincide with the axes (x,y,z) of a rectangular coordinate system, the x axis designating the crystal´s optical axis (the c-axis).A bias voltage is applied between the crystal´s boundaries in the x direction.A polarized laser beam, with the electric field directed along x, propagates in the z direction.Once a stationary state is reached, the beam´s electric field can be written as: with  =  0   = 2  / 0 , where   is the unperturbed refractive index along the extraordinary c-axis, and  0 is the carrier wavelength in vacuum.Field  ⃗ (, ) satisfies the Helmholtz equation, where   ′ is the full refractive index along the c-axis, (  ′ ) 2 =   2 −   4     [34,35], where   is the electro-optic coefficient, and   ⃗⃗⃗⃗⃗⃗ =    ̂ is the space-charge field induced in the crystal.Substituting the form of  ⃗ , defined as per Eq. ( 1), in Eq. ( 2), and using the paraxial approximation, we obtain an equation of the NLS type: In a biased photovoltaic pyroelectric photorefractive crystal, the space-charge field has three contributions [6]: where  =   || 2 /(2 0 ) is the beam´s intensity (with  0 = √ 0 / 0 ),   is the dark irradiance,  ∞ = ( → ±∞), and  0 = ( → ±∞).The background value  0 is determined by the external voltage,   depends on the polarization of light and photovoltaic coefficient,   is the pyroelectric field, and coefficient   is determined by properties of the crystal.The pyroelectric contribution to   (the last term in Eq. ( 5)) is adopted according to the approximation introduced in [36].
Finally, we consider the illumination of the crystal corresponding to  ∞ =  = 0, hence Eq. ( 6) reduces to the GNLSE: This equation predicts the propagation of stationary spatial solitons in the photorefractive waveguide [6].
Stationary soliton solutions to Eq. ( 7) with a real propagation constant  are looked for as: where () is a real function.In particular, the main characteristic of the soliton is its integral optical power, see also Eq. ( 12) below.

The Lagrangian structure
In terms of the appropriate Lagrangian density, Eq. ( 7) is written as the variational Euler-Lagrange equation: According to the Noether's theorem, the invariance of the action integral, corresponding to the Lagrangian density (10), with respect to translations along directions ζ and s, and an arbitrary shift of the phase of the complex amplitude U, implies the conservation of three dynamical invariants, viz., the Hamiltonian, momentum, and optical power (norm) [cf.Eq. ( 9)]: with densities: The Lagrangian representation of Eq. ( 5) suggests the possibility of applying the Anderson´s VA [21,22] to predict the propagation of a solitary-wave input in the photorefractive crystal.
To apply the VA, one needs to choose an adequate variational ansatz for a soliton (strictly speaking, it is a solitary wave, rather than a soliton in an integrable system).The traditional Gaussian, similar to the one used in Refs.[21,[28][29][30][31] for NLS equations with polynomial nonlinearities, would not allow the necessary analytical calculation of the resulting integral of the nonpolynomial Lagrangian density (10).Instead, two ansätze lead to VA in a tractable analytical form.The first one is based on the hyperbolic secant: where A, a, b, and c are, respectively, the soliton´s amplitude, width, phase, and chirp.Substituting this ansatz in the Lagrangian density (8), we calculate the corresponding averaged Lagrangian as: The result being: where  =  2 , and we define: with the dilogarithm Li 2 (also called the Spence´s function), [37], which is a special case of the polylogarithm Li  for  = 2.Note that the arguments of Li 2 in Eq. ( 19) are negative, in agreement with the fact that the integral definition ( 20) is well defined for t < 1.Note that the dilogarithm appeared only in few previous works dealing with optical solitons [38,39].Therefore, it is worthy to note that this special function naturally arises in the study of solitons in photorefractive media, in the present context.The averaged Lagrangian (18), gives rise to Euler-Lagrange equations for the variational parameters (), (), () and (), which are defined in ansatz (16): (this equation implies, as usual, the power conservation), where: The primary use of the Euler-Lagrange equations ( 21)-( 24) is to predict the evolution of the solitary wave.These equations also make it possible to apply the Vakhitov-Kolokolov (VK) criterion [40][41][42][43][44] to predict the stability of the soliton solutions of the system ( 21)- (24).
For stationary solutions (alias fixed points) of the variational equations we set Eq. ( 22), which implies, as usual, that the chirp must vanish,  = 0.Then, setting Eq. ( 24), leads to a certain VA-predicted relation between the amplitude and the width in the stationary solution: where the form of the function () can be obtained from Eq. ( 24).Further, Eq. ( 23) produces the phase of the stationary solution as: where the propagation constant can be written as: or, alternatively: taking into regard that Eq. ( 23) produces k as a function of A and a, and a can be expressed in terms of A as per Eq. ( 26).
To apply the VK criterion, one needs to obtain the sign of the derivative /, where P is the optical power defined in Eq. ( 9).Using the form of the ansatz (16), and removing the width a in favor of the amplitude A pursuant to Eq. ( 26), we obtain: Then, taking into account Eq. ( 29), we obtain: where ´() = /.
As an example, Fig. 1 displays the VK derivative /, as obtained from Eq. ( 31) with the following (physically relevant) values of the coefficients in Eq. ( 7 16), as produced by Eq. ( 31) with the coefficients fixed as per Eq. ( 32).
It is seen in Fig. 1 that the condition   > 0 holds for all values of k, and hence the VApredicted family of stationary soliton solutions may be completely stable.In the following section we present some numerical solutions corroborating this prediction.An alternative ansatz relies upon the square root of sech, instead of sech in Eq. ( 16), to approximate the soliton's shape: where the parameters , , , and  have the same meaning as in Eq. ( 16).Substituting this ansatz in the Lagrangian density (10) and calculating the integral (17), we obtain the averaged Lagrangian: where  =  2 as above, and: Note that the averaged Lagrangian (34) also includes the dilogarithm Li 2 .The occurrence of this special function in both averaged Lagrangians, ( 18) and (34), is due to the presence of the logarithmic term in the Lagrangian density (10).Consequently, it may be expected that Li2 will appear whenever the VA is applied to other models which contain a logarithmic term in the Lagrangian.
Having the Euler-Lagrange equations ( 36)- (39) we can apply the VK criterion to predict the stability of the soliton solutions.Following a procedure similar to that which was employed above to produce Fig. 1, and fixing the same parameters as in Eq. ( 32), we obtain a plot for the function / which has a form which is almost identical to the one shown in Fig. 1.Therefore, ansatz (33) also predicts VK-stable solitons as stationary solutions of Eq. (7).
To close this section we would like to emphasize that the VA is a powerful tool in the study of solitons and breathers because it allows one to transform a complex nonlinear partial differential equation (NLPDE) into a system of first-order ordinary differential equations (ODEs), which is definitely much easier to solve than the underlying NLPDE.The Euler-Lagrange equations provided by the VA may also permit us to obtain, in a compact way, approximate results that would be quite difficult to obtain if we dealt directly with the original NLPDE.For example, in this section we were able to estimate the stability of the solitons of Eq. ( 7) by means of the VK criterion, because the variational equations ( 22)-( 26) permitted us to calculate the sign of the derivative dP/dk.

Variational results produced by the sech ansatz (16)
We start the analysis by considering the same coefficients , β, θ and γ from Eq. ( 32), that we used above to generate Fig. 1.Actually, these coefficients adequately represent the generic case.
First, Fig. 2 demonstrates the soliton's shape, as produced by ansatz ( 16) with parameters obtained from the numerical solution of the Euler-Lagrange equations ( 21)- (24), with the initial conditions corresponding to a fixed point of the equations, which is  = √10 and  = 7/100, the respective ansatz being (0, ) = √10 sech (100 /7) . ( The stationary propagation observed in Fig. 2 corroborates that this fixed point is a stable constant solution of equations ( 21)-( 24).(47).The counterpart of this picture produced by full simulations of GNLSE ( 7) is presented below in Fig. 16.
Another relevant example of breathers is produced by Eqs. ( 21)-( 24) with a parameter set which differs from ( 42) by a smaller value of coefficient  in front of the cubic term:  = 31.5965, = 1.0149,  = −13.5328, = 1, (44) and initial values of the variational parameters corresponding to the following input: cf. Eq. ( 43).Similar to the example displayed in Figs. 3 and 4, the amplitude and width of expression (45) do not correspond to a fixed point of Eqs. ( 22) and (24).In this case, the solution of the Euler-Lagrange equations ( 21)-( 24) produces the results presented in Figs. 5  and 6.They exhibit oscillations of the breather with period  ≈ 0.411, which is larger by a factor ≃ 3.34 than in Figs. 3 and 4.
Fig. 5.The same as in Fig. 3, but for parameters (44) and the input corresponding to Eq. ( 45).
Lastly, Figs. 7 and 8 demonstrate an example of a robust breather featuring much slower oscillations, with period  ≈ 1.888, which is larger by a factor ≃ 4.6 than in Figs. 5 and 6.This solution of Eqs. ( 21)-( 24) was obtained for a parameter set:  = −6.7658, = 14.078,  = −13.5153, = 1 (46) and the same initial input (41) as used above, which in this case does not represent a fixed point of Eqs. ( 22) and ( 24) either.Fig. 7.The same as in Fig. 3, but for parameters (46) and the input corresponding to Eq. ( 41).The values of coefficients (46)   Next, we produce VA predictions based on the alternative analytically tractable ansatz (33), with the evolution of the variational parameters governed by the Euler-Lagrange equations ( 36)- (39).To begin with, in Fig. 9 we produce a stationary soliton generated by initial values of the variational parameters corresponding to the input: (0, ) = √10 sech(13.11) (48) whose amplitude and width correspond to a fixed point of Eqs. ( 37) and ( 39), and the respective set of coefficients in Eq. ( 7) being: Next, to produce a breather shown in Figs. 10 and 11, we take the same input ( 48), but with different coefficients in Eqs. ( 7) and ( 36)-( 39 In this case, parameters of input (48) do not correspond to a fixed point of Eqs. ( 37) and ( 39), and consequently the simulations of the variational equations ( 36)-( 39) produce a breather with a period of stable oscillations λ ≈ 0.24.
Lastly, a stable breather with a very large oscillation period,  ≈ 5.8, which is ≈24 times larger than in Figs. 12 and  in Eqs. ( 36)-( 39).This slowly oscillating breather (once again, initiated by an input which does not correspond to the fixed point of Eqs. ( 37) and ( 39)), is displayed in Figs. 14 and 15.
Fig. 14.Breather produced by ansatz (33).In this case, the variational parameters evolve according to the Euler-Lagrange equations ( 36)-( 39), with the coefficients shown in Eq. (52) and the initial condition given in Eq. ( 51).The values of coefficients ( 52  The values of coefficients ,  and  used to obtain the solutions shown in Figs.9-15 are translated, using the material characteristics given in book [18], into the values of the physical parameters of the photorefractive crystal given in Eqs.(47), along with   = −2 / and  0 = −3.99518/.

Direct simulations of the evolution of solitons and breathers
The predictions of VA reported above have been verified by means of the full numerical solution of Eq. ( 7).These numerical solutions were obtained by using a 4 th order Runge-Kutta algorithm to advance along ζ, and approximating the second derivative   with finite differences including first and second neighbours.To begin with, we solved Eq. ( 7) using the same coefficients [Eq.( 32)] and initial conditions [Eq. ( 41)] that we used to produce the variational solution shown in Fig. 2. The corresponding numerical solution of Eq. ( 7) is displayed in Fig. 16.It is seen that the full simulations produce a stable soliton in perfect agreement with the VA prediction.The surfaces seen in Figs. 2 and 16 look quite similar, but in order to appreciate to what extent these surfaces are really similar, in Fig. 17 7) with the same coefficients [see Eq. ( 32)] and the same input [Eq.( 41)] which were used above to produce the VA solution displayed in Fig. 2. Next, in Figs.18 and 19 we show numerical solutions of Eq. ( 7) which corroborate the existence of robust breathers which are quite close to those predicted above in Figs. 3 and 5 for the same parameters and initial conditions.It is worthy to note that the VA-predicted and numerically found oscillation periods are very close.Fig. 18. Results of the numerical simulations of Eq. ( 7) with the same coefficients [see Eq. ( 42)] and the same input [Eq.( 43)] which were used above to produce the VA solution displayed in Fig. 3. Fig. 19.Results of the numerical simulations of Eq. ( 7) with the same coefficients [see Eq. ( 44)] and the same input [Eq.(45)] which were used above to produce the VA solution displayed in Fig. 5.
To appreciate to what extent the variational solution shown in Fig. 5 is similar to the direct numerical solution shown in Fig. 19, in Fig. 20 we present the profiles |(,  = 1)| obtained from the variational solution shown Fig. 5 [with a thin line], and the profile obtained from the direct numerical solution of Eq. ( 7) [with black dots].It is seen that the agreement between both profiles is very good, thus confirming that the variational method indeed produces solutions which are close to the direct numerical solutions of Eq. ( 7).The breathing solutions shown in Figs.18 and 19 were produced by using initial conditions which were very close to exact breathers.A natural question is what will happen if one uses an input which is not too close to an exact breather.To address this issue, in Fig. 21 we show the numerical solution of Eq. ( 7) obtained with the input: and the coefficients given in Eq. ( 46).It is observed that the input evolves into a breather, even though the initial amplitude of the pulse is significantly higher than the breather´s amplitude.This result implies that even if the breather is significantly perturbed (by increasing its amplitude and reducing its width), the pulse returns to the breather´s shape, indicating that this breather is a robust solution.
All the breather solutions produced by the simulations of Eq. ( 7) demonstrates their robustness against decay into small-amplitude radiation waves.Strictly speaking, the emission of radiation, in the form of slow decay of an excited state of a fundamental soliton [45] is generated by any equation of the NLS type (a well-known fact is that the integrable NLS equation with the cubic self-focusing gives rise to exact breather solutions in the form of N-solitons, which are produced by the input in the form of a fundamental soliton multiplied by an arbitrary integer N = 2, 3, 4, … [46]).Nevertheless, in our simulations the decay is a very weak effect, which remains virtually invisible.7) for the coefficients shown in (46), and the initial condition (53).We can use the numerical solutions of Eq. ( 7) to corroborate that the VK criterion predicts the stability of the soliton solutions of Eq. (7).To this end, we looked for solitons corresponding to the same set of coefficients {, , , } as chosen in Eq. ( 49), and four different values of the propagation constant: k = 21.99,31.42,62.83 and 100.53.In each case, we computed the integral power P [see Eq. ( 12)], as shown in Fig. 22.This figure shows that dependence () has a positive slope (i.e., / > 0), hence the VK criterion indicates that these solitons are stable solutions, which is corroborated by the simulations of Eq. (7).Fig. 22. Optical power P vs. the propagation constant k of the soliton solutions of Eq. ( 7) with the coefficients given in Eq. (49).
Figures (17)- (19) show that the separation  between the breathers´ peaks changes if we modify the coefficients of Eq. ( 7), and it is worthy to note that that the VA presented in Figs.
3, 5 and 7 predicts similar changes in λ.To evaluate to what extent the VA-predicted changes in λ are similar to those observed in the direct numerical solutions of Eq. ( 7), we have calculated how  depends on  according to the direct numerical solutions of Eq. ( 7) and, on the other hand, according to the variational equations ( 21)- (24).We focus here on the parameter , as this is the coefficient associated to the cubic nonlinearity of Eq. ( 7), which is responsible for the creation and propagation of solitons in NLS-like equations.
In Fig. 23 we display the form of function () according to the direct numerical solutions of Eq. ( 7) [black squares], and according to the variational equations [stars].All the results shown in this figure were obtained using the coefficients  = 31.60, = −13.53, = 1 in Eq. ( 7), and the initial condition defined in Eq. (45).
Figure 23 shows that both procedures [performed using direct numerical solutions of Eq. ( 7), or the variational equations], show that λ decreases as  increases.This behaviour can be understood if we observe that the definition of  [below Eq. ( 6)], combined with Eqs. ( 4) and ( 5), imply that  will increase if  increases.And if  (the refractive index change due to the photorefractive effect) increases, the light velocity in the crystal decreases.Consequently, if ∆ is the time interval between two successive maxima in the light intensity of a breather, the distance ∆, advanced by the light during this time interval, diminishes if the velocity of light () decreases.In other words,  = ∆ will decrease.This is the physical interpretation of Fig. 23.It should be mentioned, however, that a change in  is not the only possible cause for a change in λ.In particular, it is known that light absorption may also change the value of λ [47].7), and the stars were obtained from the variational solution, using the ansatz (16).The input is taken as per Eq. ( 45), and other coefficients in Eq. ( 7) are:  = 31.60, = −13.53and  = 1.
The numerical solutions of Eq. ( 7) presented in this subsection, as well as the variational solutions presented in the subsections 4.1 and 4.2, show that Eq. ( 7) permits the propagation of solitons and breathers.Consequently, a natural question is if there is a criterion that permits one to forecast if a given solitary-wave input will propagate as a soliton, or it will evolve as a breather.In the case of the cubic NLS equation the integrability of this equation permits one to answer this question.More precisely: if we consider the input in the form (,  = 0) =  sech(/), the Zakharov-Shabat equations predict the number of solitons that will be generated as the number of bound eigenstates in a linear problem (in a similar way as in the usual quantum mechanical problem with a given potential well).Namely, it is well-known [46] that, if the parameters of the input satisfy the inequality  < 1/2, the potential well has no bound states, and no solitons appear.In the interval 1/2 <  < 3/2 the eigenvalue problem has exactly one bound state, hence exactly one soliton appears.However, when  > 3/2, the second bound state appears, which implies the creation of the superposition of two solitons, i.e., precisely the simplest breather.In other words, in the integrable-NLS case, the character of the output (i.e., the generation of a soliton or a breather), is determined by the size of the product .For a nonintegrable equation such as Eq. ( 7), there are no exact results like that, although the situation remains qualitatively similar: for relatively small values of , the input is able to maintain the self-trapping of a single soliton, but for larger values of , the input cannot avoid the generation of at least one additional soliton, whose interaction with the original soliton builds a breather.It should be mentioned, however, than the determination of precise criteria which may be helpful to predict if a given solitary-wave initial condition for Eq. ( 7) will evolve as a soliton, or as a breather, is outside the scope of this article.

Breathers modulated by beats
The solutions produced by the simulations of Eq. ( 7) in Figs. 18, 19 and 21 may be considered as regular breathers, as they exhibit strictly periodic oscillations.On the other hand, nonintegrable systems (such as the NLS equation combining the cubic self-defocusing and harmonic-oscillator trapping potential [48]) may give rise, under special conditions, to dynamical states featuring quasi-periodic oscillations.In accordance with this expectation, Eq. ( 7) gives rise, in addition to the regular breathers presented above, to more sophisticated spatially self-trapped modes in the form of breathers featuring long-period beats, such as one displayed in Fig. 24.It corresponds to coefficients:  = −6.7658, = 67.6584, = −13.2643, = 1 (54) in Eq. ( 7), and the input: Fig. 24.Beating breather, obtained as a numerical solution of Eq. ( 7) with parameters (54) and input (55).The values of coefficients (54) correspond to the physical parameters given in Eqs.(47),   = 40 /,  0 = −3.92098/ and   = −2 /.
To identify the character of the beats exhibited by the breather displayed in Fig. 24, it is relevant to calculate the power spectrum of oscillations of their amplitude, | ̃()| 2 , where  ̃() is the Fourier transform of (,  = 0),  standing for the respective spatial frequency.The power spectrum for the amplitude oscillations of the breather shown in Fig. 24 is plotted in Fig. 25.In this case, the Fourier transform of |(,  = 0)| was computed using 2 14 values of this function of , corresponding to 16384 equally spaced points in the interval of 0 ≤  ≤ 7.5.As ζ is a dimensionless variable, the Fourier transform of |(,  = 0)| is a function of a dimensionless spatial frequency.It is observed that the discrete spectral components at frequencies ±1.09, ±2.19, ±3.28 represent anharmonic oscillations of the breather, while the dc component concentrated around the zero frequency accounts for longperiod beats.
The results presented in this section, and in the previous ones, show that breathers in photorefractive crystals constitute a subject deserving further investigation.In this connection, it is relevant to mention that the study of breathers in other systems, such as linear and nonlinear phononic crystals, and granular crystals [49][50][51], and the study of discrete breathers [52][53][54], are also interesting fields.

Splitting breathers
Another example of a breather with slowly modulated intrinsic oscillations is displayed in Fig. 26.It is obtained as a numerical solution of Eq. ( 7) with the same values of coefficients as in Eq. ( 53) and the input: (,  = 0) = √10 sech (8.70). ( Unlike the dynamical scenarios demonstrated above, this breather splits in two parts: a small solitary wave that moves to the left (towards s < 0), and the principal one, which moves to the right with a small velocity.While the small solitary wave which moves to the left is not seen in Fig. 26 because it is located behind the principal wave, it is visible under the different angle used in Fig. 27.Actually, splitting is one of generic outcomes of the evolution of chirped pulses in equations of the NLS type [55,56].Splitting of a single pulse into multiple solitary waves is possible too.A typical example of this type is displayed in Fig. 28, where the initial pulse is fragmented into five solitons with weak internal oscillations.In this case, the coefficients of Eq. ( 7) are the same as in Fig. 24, i.e., they are given as per Eq. ( 54), and the input is: (,  = 0) = √10 sech (4.50).
(57) Fig. 28.Splitting of the input pulse in five secondary ones, exhibited by the numerical solution Eq. ( 7) with the coefficients shown in Eq. ( 54), and the initial condition (57).These coefficients correspond to the physical parameters given in Eqs. ( 47 The two examples presented in Figs.27 and 28 show that the solutions of Eq. ( 7) corresponding to localized initial conditions may produce solutions which split into two, or more, solitary waves.This is not a surprise, as it is well-known that in the case of the integrable NLS equation, a sech input, depending on amplitude  in front of it [(,  = 0) =  sech ()], may generate a soliton, or produce fission into a set of several solitary waves.Without any additional perturbation, the multi-soliton state would form a breather, while any perturbation splits it into a set of separating fundamental solitons.The number of solitons, , is determined by what is usually called "area"  of the input, i.e., the integral of (,  = 0).This integral should not be confused with the norm, which is the integral of [(,  = 0)] 2 , and unlike the norm, the area is not a dynamical invariant of the NLS equation.In the case of a nonintegrable equation, such as Eq. ( 7), there is no exact criterion to predict the number of solitary waves generated by the splitting of the input.However, a qualitative similarity still holds: also in the case of Eq. ( 7) the number of fundamental solitons created by the splitting of a sech input increases as the area  of the input increases.We can observe that the value  57 = 2.21, corresponding to the input (57), is nearly twice as large in comparison to  56 = 1.14 for the pulse (56), and, as a consequence, the number of waves generated by the splitting increases from 2 to 5. It is interesting to observe that in the integrable case, if an initial pulse with area  splits into 2 solitary waves, a broader pulse with area 2 would split into 2 × 2 = 4 new waves.Therefore, the splitting of the pulse (57) into five new waves (instead of four), is an interesting manifestation of the nonintegrability of Eq. (7).On the other hand, an essential similarity with the integrable case, is that the initial pulse (57) splits into a set of solitary waves with unequal amplitudes.Therefore, the fission of solitary waves in the case of Eq. ( 7) presents some similarities to the results observed in the integrable cubic-NLS case, but the similarity is not exact due to the fact that Eq. ( 7) is not integrable.

Conclusion
The results presented in this communication show that solitons and breathers exist in biased photovoltaic pyroelectric photorefractive crystals, where the light propagation obeys Eq. (7).The application of the VA (variational approximation) predicts the existence of breathers with different shapes and oscillation periods.The VK (Vakhitov-Kolokolov) criterion correctly indicates that these solitons are stable solutions.The VA is elaborated on the basis of two different ansätze, which are defined in Eqs. ( 16) and (33).In both cases the averaged Lagrangians, given by Eqs. ( 18) and (34), respectively, involve the dilogarithm function Li 2 .Therefore, we refer to these modes as dilogarithmic breathers.The existence of these VApredicted breathers is then confirmed by full numerical solutions of Eq. ( 7).The numerical results demonstrate that, in addition to the regular (periodic) breathers, such as those shown in Figs. 18, 19 and 21, Eq. ( 7) also gives rise to robust breathers exhibiting beats (long-period modulations) in their intrinsic dynamics, such as the one shown in Fig. 24.The spectrum of oscillations of the beating breathers is nearly discrete, thus confirming that the oscillations are indeed subject to long-period quasi-periodic modulations.Another class of numerical solutions demonstrate spontaneous splitting of the input pulse into two or several secondary quasi-soliton pulses, as shown in Figs.27 and 28.Systematic simulations of the splitting dynamics is a subject for a separate work, which is currently in progress.
we show the form of the profile |(,  = 20)| given by the numerical solution of Eq. (7) [black squares], and the profile obtained from the variational solution [the thin continuous line].It is seen that the agreement between the two solutions is excellent.

Fig. 16 .
Fig. 16.Results of the numerical simulations of Eq. (7) with the same coefficients [see Eq. (32)] and the same input [Eq.(41)] which were used above to produce the VA solution

Fig. 27 .
Fig.27.The same as in Fig.26, but observed from a different direction.