Variational Approximations of Bifurcations of Asymmetric Solitons in Cubic-Quintic Nonlinear Schroedinger Lattices

Using a variational approximation we study discrete solitons of a nonlinear Schroedinger lattice with a cubic-quintic nonlinearity. Using an ansatz with six parameters we are able to approximate bifurcations of asymmetric solutions connecting site-centered and bond-centered solutions and resulting in the exchange of their stability. We show that the numerically exact and variational approximations are quite close for solitons of small powers.

1. Introduction. The variational approximation (VA) has long been used as a semi-analytic technique to approximate solitary wave solutions of nonlinear evolution equations with an underlying Hamiltonian structure [13]. There have been a number of papers exploring the VA with four parameters as a relevant approximation of localized modes in discrete nonlinear Schrödinger (DNLS) equations [6,14,19]. Kaup [10] extended the variational approximation with six parameters that allowed him to construct not only site-centered solutions (also called on-site solitons) from [14] but also the bond-centered solutions (solitons centered at a midpoint between two adjacent sites also known as inter-site solitons).
Site-centered and bond-centered solitons were recently considered in the context of the DNLS equations with competing cubic focusing and quintic defocusing nonlinearities both in the space of one [3] and two [4] lattice dimensions. It was found that the two branches exchange their stability while continued with respect to the underlying parameters. A salient feature of this stability exchange is that the two branches of site-centered and bond-centered solitions do not intersect directly but are connected by an intermediate branch of asymmetric solitons. It was argued in [4] that the discrete solitons have enhanced mobility near the regimes of stability inversion. These properties were originally discovered in the DNLS equations with a saturable nonlinearity both in the space of one [8] and two [20] dimensions as well as for the DNLS equations with on-site and next-site cubic nonlinearities in the space of one dimension [18]. It was recently shown for the same model in the space of two dimensions [17] that stability inversion and asymmetric solutions may not lead to enhanced mobility of discrete solitons if the bifurcation points of the solution branches are widely separated in the parameter space.
It is the purpose of this work to apply Kaup's variational method with six parameters from [10] to explain bifurcations of asymmetric solutions and stability exchange of site-centered and bond-centered discrete solitons in the context of the one-dimensional cubic-quintic DNLS equation. In this sense, our work is a complement to the previous paper [3] where discrete solitons were constructed numerically using a dynamical reduction. Dark solitons and staggered solutions of the cubicquintic DNLS equation were recently studied numerically in Refs. [15] and [16] respectively.
We consider a discrete nonlinear Schrödinger equation with a cubic-quintic nonlinearity in the form, where ψ n (t) : R + → C and (C, B, Q) are real-valued parameters. Nonlinear Schrödinger lattices have proved to be relevant models in a variety of contexts (see reviews in [7,11]), including the description of optical pulses in one-dimensional waveguide arrays [5]. In this application, the quantity |ψ n | 2 represents the intensity of the electric field of waveguide n, C > 0 represents coupling strength between adjacent waveguides and (B, Q) measure the nonlinearity strength. A large portion of the literature is dedicated to Eq. (1) with Q = 0, which would correspond to a medium with a Kerr nonlinearity. Recent experimental results [2,9,21] have shown that the nonlinear response of some materials is better fit with an additional competing quintic nonlinearity, i.e. B, Q > 0. This lends relevance to studying the cubic-quintic DNLS equation in the form (1). The cubic-quintic DNLS equation (1) has two conserved quantities, namely the power, and the Hamiltonian, Steady-state solutions have the form ψ n = u n e −iµt , n ∈ Z, where µ ∈ R and u n ∈ R are found from the stationary DNLS equation, We seek localized solutions of the stationary DNLS equation (4) which in turn correspond to discrete solitons. Spectral stability of the steady-state solutions is studied with the linearization ansatz, ψ n (t) = u n + (v n + iw n )e λt + (v * n + iw * n )e λ * t e −iµt , n ∈ Z, which leads to the spectral problem, n v n + 5Qu 4 n v n = −λw n , −µw n − C(w n+1 + w n−1 − 2w n ) − Bu 2 n w n + Qu 4 n w n = λv n , n ∈ Z. (5) We are looking for nonzero solutions of the linearized system in L 2 (Z, C 2 ) called eigenvectors. The steady-state solution is called unstable if there exists at least one eigenvector for which Re(λ) > 0. Otherwise, the solution is called spectrally stable. Steady-state solutions are considered in Section 2. Spectral stability of these solutions is studied in Section 3. For both problems, we compare results of the variational approximations and the direct numerical approximations.
(6) To find approximate solutions for discrete solitons, we pose a trial function in the form, where each of the six parameters are dependent on t. Substituting Eq. (7) into the Lagrangian (6) and evaluating the sums yields the effective Lagrangian, where χ = 2n 0 − 1 and, Since the center n 0 of the ansatz (7) can be arbitrarily chosen on [0, 1] module the discrete group of translations in n ∈ Z, we shall consider χ on [−1, 1] (this restriction was used already in the derivation of (8)). The solution with χ = 0 is centered between lattice sites and hence is called bond-centered. The solution with χ = ±1 is centered on a lattice site and hence is called site-centered. Solutions for χ ∈ (−1, 0) ∪ (0, 1) are called asymmetric.
According to the variational principle, the effective Lagrangian L eff achieves critical values at the Euler-Lagrange equations, where p j represents a parameter of ansatz (7). Varying α yields the conservation law, Before writing the remaining equations, it will be more convenient to make use of the fact that we seek steady-state solutions. This corresponds to where µ is the parameter of the steady-state solution. With this assumption, the equations corresponding to variation of k and β are identically satisfied. Varying χ and η and making use of Eq. (11) leads to the following two equations respectively, and, Note these equations with Q = 0 correspond to those in Ref. [10]. When χ = 0 (bond-centered solutions), Eq. (12) is identically satisfied and Eq. (13) becomes, which is easily solved for A 2 , giving an existence condition in terms of the parameter η. There exists exactly two bond-centered solutions (denoted by B 1 and B 2 ) for η ∈ (0, η cr ), which disappear as a result of the saddle-node bifurcation at η = η cr . See Fig. 1, where solutions of Eqs. (12) and (13) are plotted for C = 0.1, B = 2, and Q = 1 with η cr ≈ 1.56. We note, for larger C steady-state solutions become smoother (they approach the continuous counterpart) and so the ansatz (7), which is based on an exponential cusp, becomes irrelevant. For χ = 0, the term in parenthesis of Eq. (12) can be used as condition for A 2 , which can then be substituted into (13) yielding a root finding problem in (η, χ) parameter space. We find exactly three pairs of solutions for C = 0.1, B = 2, and Q = 1, which are denoted by S 1 ,S 2 , and A 1 in Fig. 1. Branches (S 1 , S 2 ) and (S 2 , A 1 ) are connected by means of the saddle-node bifurcations. Branch S 1 arises from χ = 0 as a result of the supercritical pitchfork bifurcation, while branches S 2 and A 1 arise from χ = 0 as a result of the subcritical pitchfork bifurcation.
The solution branches S 1 and S 2 approach χ = 1 rapidly, so that they correspond to the site-centered solitons slightly disturbed by the variational approximation. The solution branch A 1 is a true asymmetric solution that connects the bondcentered and site-centered solitons. Only one branch S 1 of site-centered solutions and only one branch B 1 of bond-centered solutions exist in the cubic case Q = 0, where these two branches extend for any η > 0.
To test the validity of the variational approximations we compare them against direct numerical solutions of the stationary DNLS equation (4). The "numerically exact" solutions are obtained using a Newton method with the VA solutions as initial seeds. The profiles of the discrete solitons that are obtained via the VA and the corresponding numerical solutions are shown in Fig. 2.  There are two principle branches, one corresponding to site-centered solutions (labeled S j , j = 1..4) and one for bond-centered solutions (labeled B j , j = 1..3). The two principle branches are connected via asymmetric solutions (labeled A j , j = 1, 2) at points of stability change. Only the branches S 1 , S 2 , A 1 , B 1 , and B 2 , are captured by the VA. Differences between the VA and numerical solutions are more visible in the zooms labeled 'a', 'b' and 'c'. The branch labeledB 3 corresponds to the VA that fails to capture B 3 .
It is more instructive for comparison to plot the solution branches in the (µ, M ) plane, where M is defined in Eq. (10) and µ = − dα dt , see Fig. 3. The predicted stability is also depicted (see Sec. 3 for details on stability computations).
Besides the five solutions predicted within the variational approximation, there exists additional branches S 3 , S 4 , B 3 , A 2 of solutions of Eq. (4) that appear on Fig. 3. See Fig. 4 for profiles of these "wide" solutions. The pattern of Fig. 3 suggests existence of an infinite number of additional branches of discrete solitons for large power M with a characteristic snaking behavior. The snaking behavior was also found in the higher dimensional cubic-quintic DNLS equation [4] and the Swift-Hohenberg equation [1,12]. It seems the role a higher order dispersion plays The bifurcations shown in Fig. 3 look as if the asymmetric solutions connect exactly at the turning points of bond-centered or site-centered solutions. This is not the case. Rather, there exists a pair of saddle-node and pitchfork bifurcations and the stability exchange occurs at the pitchfork bifurcation where the asymmetric solutions emanate. Interestingly, the VA is able to accurately capture this subtle bifurcation scenario. In Fig. 5 a plot of such a bifurcation pair is shown. This particular pair represents a saddle-node bifurcation of a short (B 1 ) and tall (B 2 ) bond-centered solution and a pitchfork of two asymmetric solutions (A 1 ) and the short bond-centered solution (B 1 ). The agreement between the variational and numerical approximations is impressive. The stability is also correctly predicted.
On the other hand, bifurcations including the site-centered solutions are not captured as well. This should not be surprising since site-centered solutions are never truly represented by the VA (see Fig. 1). We point out other minor discrepancies between the variational approximations and the numerical results in the zooms of Fig. 3. Zoom (a) shows where the VA solutionB 3 departs from the corresponding numerical solution B 3 . This is expected as the ansatz (7) is only valid for "narrow" solutions, i.e. those that have at most two initially excited sites. Zoom (b) shows that the asymmetric solution A 1 is connected to the site-centered solution S 2 (opposed to S 3 ) and is underestimated. The wider short site-centered solution S 3 is not captured at all by the VA (as expected). Finally, zoom (c) shows that the VA falsely predicts collision of the bond-centered B 1 and site-centered S 1 solutions for a non-zero value of M , similar to the cubic case discussed by Kaup in [10].
3. Variational approximations of spectral stability. In order to determine stability within the VA we return to the variational equations (9) but this time without the assumption that β = k = dχ dt = dη dt = 0. Let x = (β, χ, η, k) T represent the four parameters of ansatz (7) after Eqs. (10) and (11) are used. To perform a linear stability analysis we substitute, into the four variational equations, where the steady-state solution is defined by x 0 = (0, χ 0 , η 0 , 0) T and (χ 0 , η 0 ) satisfy Eqs. (12) and (13). Keeping only the terms linear in ǫ leads to the generalized eigenvalue problem, where the entries of the 4 × 4 matrices A and B are given by, a 11 = 0, a 21 = 0, a 32 = 0, a 42 = 0, a 33 = 0, a 43 = 0, Here we have defined, and each entry is evaluated at the fixed point x 0 . For the form of the perturbation chosen, an eigenvalue with Re(λ) > 0 indicates instability of the steady-state solution x 0 . Since A and B are real-valued matrices, both λ and λ * are eigenvalues. Since A and B has a clear block structure, if λ is an eigenvalue with the eigenvector y, then −λ is an eigenvalue with the eigenvector S y, where S = diag(1, −1, −1, 1). Therefore, eigenvalues of (15) occur as pairs of real or imaginary eigenvalues or as quartets of complex eigenvalues. A typical example of eigenvalues of (15) is shown on the left panel of Fig. 6 for the unstable solution B 1 with a pair of real and a pair of imaginary eigenvalues.
The stability of the VA solutions corresponding to Fig. 3 was computed from eigenvalues of the generalized problem (15). The spectrum for each of the variational solutions shown in Fig. 3 is plotted in the left panels of Fig. 7. Only the bondcentered solution B 1 has a branch that is predicted to be unstable. The bottom left panel of Fig. 7 is a zoom of small eigenvalues in the top left panel. Where the asymmetric solution A 1 meets the bond-centered solutions B 1 and B 2 corresponds to stability exchange, whereas the connection to the site-centered solution S 2 occurs at λ 2 < 0 and hence no stability change takes place. Note that each branch of VA solutions has exactly two pairs of real or imaginary eigenvalues of the generalized problem (15).
A linear stability analysis of the "numerically exact" solutions was also carried out from the spectral stability problem (5) using the procedure described in Ref. [3]. The numerical approximation of the spectrum for the unstable bond-centered solution B 1 is shown on the right panel of Fig. 6. Note that there is always a double zero eigenvalue in the stability problem (5) which is related to the gauge symmetry. This double zero eigenvalue corresponds to the eigenvector (v n , w n ) = (0, u n )  Figure 6. Eigenvalues corresponding to the bond-centered solution B 1 in the variational system (left) and the full system (right). and the generalized eigenvector (v n , w n ) = (∂ µ u n , 0) of system (5). This double zero eigenvalue is captured by the variational approximation and it results in the conservation law (10) and the variational equation (11). On the other hand, the VA gives two pairs of eigenvalues and only the real pair persists in the numerical approximation. The purely imaginary pair of eigenvalues does not exist in the numerical solutions, where it is replaced by the continuous spectral band on the two symmetric intervals of the purely imaginary axis.
Isolated, non-zero eigenvalues for solutions B 1 , B 2 , S 2 , S 3 , and A 1 are shown in the right panels of Fig. 7. We note that no isolated eigenvalues of S 1 exist. The overall "look" of the top panels of Fig. 7 are quite similar, but more importantly, the stability is correctly predicted. Note that the right panels of Fig. 7 contains fewer eigenvalues than the left panels of Fig. 7 because some of the purely imaginary eigenvalues of the VA solutions approximate the continuous spectrum of the numerical solutions. An interesting difference is seen in the bottom right panel of Fig. 7. In the full problem, the eigenvalue pair for the asymmetric solution A 1 vanishes when A 1 intersects with the site-centered solution S 2 .
In conclusion, we have extended the results of Ref. [10] to show that not only does the VA faithfully represent the fundamental localized modes, but is also able to correctly predict the corresponding stability for small coupling constant C and power M . We showed this in the context of the cubic-quintic DNLS equation, which exhibits a family of discrete solitons, five of which were accurately captured by the variational approximation. It would be interesting to derive the variational equations in the context of time-dependent perturbations, although, the resulting equations would be far more complex and may undermine the utility of the method. It would also be relevant to extend the results of [4] and apply this analysis to the higher dimensional cubic-quintic DNLS equation.