Exact solutions for periodic and solitary matter waves in nonlinear lattices

We produce three vast classes of exact periodic and soliton solutions to the one-dimensional Gross-Pitaevskii equation (GPE) with the pseudopotential in the form of a nonlinear lattice (NL), induced by a spatially periodic modulation of the local nonlinearity. It is well known that NLs in Bose-Einstein condensates (BECs) may be created by means of the Feshbach-resonance technique. The model may also include linear potentials with the same periodicity. The NL modulation function, the linear potential (if any), and the corresponding exact solutions are expressed in terms of the Jacobi's elliptic functions of three types, cn, dn, and sn, which give rise to the three different classes of the solutions. The potentials and associated solutions are parameterized by two free constants and an additional sign parameter in the absence of the linear potential. In the presence of the latter, the solution families feature two additional free parameters. The families include both sign-constant and sign-changing NLs. Density maxima of the solutions may coincide with either minima or maxima of the periodic pseudopotential. The solutions reduce to solitons in the limit of the infinite period. The stability of the solutions is tested via direct simulations of the GPE. As a result, stability regions are identified for the periodic solutions and solitons. The periodic patterns of the cn type, and the respective limit-form solutions in the form of bright solitons, may be stable both in the absence and presence of the linear potential. The stability of the two other solution classes, of the dn and sn types, is only possible with the linear potential.

solutions, affiliated with each of these functions; accordingly, the families are referred to as ones of the "cn", "sn", and "dn" types. In the case when only the NL pseudopotential is included, while the linear potential is absent, each family is controlled by three free parameters, two continuous and one which determines the sign of the nonlinearity. If the linear-lattice potential is included too, each family features four free continuous parameters. In the limit of the infinite period, the solutions of the cn and dn types go over into bright solitons, while the respective limit of the sn-type solutions is a dark soliton. The stability of all these solutions, both periodic and solitary (except for the dark solitons), is tested in this work by means of systematic direct simulations of the perturbed evolution.
A noteworthy feature of the exact solutions of all the three types is that they can be found in cases when the local nonlinearity periodically changes its sign, thus forming NLs with the spatially periodic alternation of self-focusing and defocusing layers. It is relevant to note that exact solutions which were recently reported in a model combining nonlinear and linear periodic potentials [4] cannot be obtained in such a case, as they are generated by a gauge transformation of the NLSE with constant coefficients.
The normalized form and physical meaning of the GPE for the atomic wave function, Ψ(x, t), and of the NLSE for the local amplitude of the electromagnetic wave in optics, are well known [36,50] (in the latter case, t is the propagation distance of the beam launched into the layered medium in the longitudinal direction, see Refs. [8,17,26,44,59]): Here g(x) and V (x) represent the nonlinear and linear lattices, respectively, with g(x) > 0/g(x) < 0 corresponding to the local attraction/repulsion between atoms in the BEC. As usual, stationary solutions to Eq. (1) with chemical potential µ are sought for as Ψ (x, t) = ψ(x) exp(−iµt), where real function ψ(x) obeys equation The corresponding expression for the energy of the matter-wave configuration is The analysis presented in this work is focused, first, on looking for special forms of periodic functions g(x) and V (x), for which Eq. (2) admits exact solutions in terms of elliptic functions. Then, the stability of the exact solutions is tested by means of simulations of the perturbed evolution of those solutions within the framework of Eq. (1). The paper is organized as follows. In Section 2, we introduce the class of exact solution of the cn type, first in the model without the linear potential, and then in the its more general version, including the potential. Results of tests of the stability of this class of exact periodic solutions are reported in Section 3. The cn patterns may be stable both without the linear potential, and in the presence of the potential. Solutions of the two other classes, of the dn and sn types, together with the results of tests of their stability, are presented in Sections 4 and 5, respectively. Both these species of the patterns require the presence of the linear potential for their stability. Section 6 reports results for bright solitons, which can be obtained as ultimate-form solutions from families of the cn and dn types. The bright solitons may be stable in the absence and presence of the linear potential alike. Dark solitons, which are limit-form solutions of the sn type, are briefly considered too in Section 6, without the stability analysis, which will be reported elsewhere. Conclusions drawn from this work are formulated in Section 7.

A. The model without the linear potential
We first consider the case of the NL pseudopotential alone, with V (x) = 0 in Eq. (1). Using the techniques elaborated in Refs. [19-22, 40, 48, 60], which are based on manipulations with the Jacobi's elliptic functions, it is possible to identify general forms of the NL structural function, i.e., coefficient g(x) in Eq. (2), which admit exact solutions to Eq. (2) with V = 0, and the form of the solution, ψ(x), itself. The first type of the solutions is based on the elliptic cosine (cn), the corresponding expressions for g(x) and ψ(x) being The spatial periods of solution (5) and NL potential (4) are, respectively, L = 4K(k)/r (6) and L/2, where k is the modulus of cn, which takes values 0 < k ≤ 1, K(k) is the complete elliptic integral of the first kind (the case of k = 1 corresponds to L = ∞, i.e., localized solutions). Generally, expressions (4) and (5) contain a set of four free parameters, viz., g 0 , r, k, b, while the remaining coefficient (modulation depth), g 1 , together with the amplitude and chemical potential of the solution, are found to be Using the obvious scaling invariance and assuming g 0 = 0, we fix |g 0 | ≡ 1 (if g 0 = 0, the above solution becomes trivial). Throughout this paper, g 0 = ±1 is kept as the sign parameter. Thus, k and r determine the period of the NL (L), g 0 is the overall sign of the nonlinearity, and b controls the depth of the spatial modulation of the periodic pseudopotential. An obvious constraint on b, which may be both positive and negative, is b > −1, which is necessary to avoid singularities in Eqs. (4) and (5). Further, the remaining scaling invariance of Eqs. (1) and (2) makes it possible to set r ≡ 1, without the loss of generality, which is fixed below, unless the linear potential is included. Thus, in the absence of the linear potential, the family of the exact solutions of the cn type depends on two free coefficients, k and b, and the sign parameter, g 0 . Comparing expressions (4) and (5), it is easy to see that maxima of the density of the periodic solutions, ψ 2 (x), which are always collocated with points where cn 2 (rx) = 1, coincide with maxima or minima of g(x)-i.e., with, respectively, minima or maxima of the NL pseudopotential-in cases when, severally, g 1 > g 0 b or g 1 < g 0 b. Then, substituting expression (7) for g 1 , we conclude that, for b < 0, the density maxima always coincide with maxima of g(x), while for b > 0 this is true under condition Otherwise, density minima (ψ = 0) are located at maxima of g(x). In fact, taking into account condition A 2 0 > 0, as it follows from Eq. (8), one can conclude that condition (10) may only hold in the case of g 0 = −1. Further, it is seen from Eq. (3) that the ground state, which realizes the minimum of the energy, must have density maxima collocated with minima of the pseudopotential, i.e., maxima of g(x) (of course, this condition is only necessary but not sufficient for the exact solution to represent the ground state). It is relevant to mention that various localized modes in the simplest model of the NL, represented just by two spots at which the self-attractive nonlinearity is concentrated, may be stable, even if they do not represent the ground state [43].
Lastly, we note that, as it follows from expression (4), Eq. (2) degenerates into the ordinary GPE, with g(x) ≡ ±1, in two cases: b = 0, or In the former case, expressions (5) and (9) go over into the usual cnoidal solution of the GPE with g(x) ≡ +1, i.e., provided that g 0 = +1 [for g 0 = −1 and b = 0, the solution does not exist, as Eq. (8) yields A 2 0 = −k 2 ]. The other degeneration condition, given by Eq. (11), if combined with Eq. (7), yields In this case, solution (5) is relevant for g 0 = −1, being [for g 0 = +1, one obtains A 2 0 < 0 from Eq. (8)]. It is easy to check that expression (14) is a straightforward exact solution to Eq. (2) with g(x) ≡ −1 and µ = (1/2) 1 + k 2 , in agreement with Eq. (9). B. The model with the linear potential A more general class of exact periodic solutions of the cn type can be found if the NL pseudopotential, introduced, as above, in the form of Eqs. (4) [but Eq. (7) should be dropped, see an explanation below], is combined with the following linear-lattice potential: The respective exact cnoidal solutions keeps the same general form as in Eq. (5), but with the amplitude and chemical potential given by This solution family depends on four continuous real parameters, viz., g 1 , r, k, b, and sign parameter g 0 . Unlike the case of V = 0, coefficient r cannot be scaled out, unless b = 0, and the value of g 1 is now another free parameter, rather than the one given by expression (7); in fact, Eq. (7) follows from Eq. (16) if one sets V 0 = b, which reduces linear potential (15) to a trivial form, V (x) ≡ 1.
C. Further analysis of the exact solutions (without the linear potential) The above solutions are meaningful provided that Eqs. (8) and (17) yield A 2 0 > 0. Because the exact solution is defined in a vast parameter space, several cases should be considered separately. In this subsection, we do that for the solutions obtained in the absence of the linear potential.
The first particular case corresponds to g 0 = +1, b > 0, see Eq. (4). Then, as follows from Eq. (8), condition A 2 0 > 0 requires For all b > 0, Eq. (19) is compatible with constraint k < 1. On the other hand, Eq. (19) is incompatible with condition (10), hence this subfamily of the exact solutions cannot represent the ground state (and, in fact, it cannot be stable, as shown below). As said above, an interesting possibility is to find exact solutions in the sign-changing NL, which makes the model drastically different from the usual NLSE. As follows from Eqs. (4) and (7), the NL-modulation coefficient, g(x), periodically changes its sign for g 1 < −1, i.e., In view of underlying condition (19), the denominator on the left-hand side of Eq. (20) is negative, hence this inequality, together with Eq. (19), identifies a parameter region in which the exact solutions are available for the sign-changing NL: It is easy to check that, for b > 0, interval (21) always exists, and, as a whole, it lies within the region of 0 < k < 1.
The case of g 0 = +1 and −1 < b < 0 in the nonlinearity-modulation function (4) can be considered too (recall the solutions with b < 0 may represent the ground state). The inspection of Eq. (8) shows that exact solution (5) exists for all k 2 < 1 at 0 < −b < 1/2, and, at 1/2 < −b < 2/3, it exists in interval cf. Eq. (19), there being no solutions at −1 < b < −2/3. Further considerations of Eqs. (4) and (7) demonstrate that the NL cannot change its sign in this case, always corresponding to the local self-attraction. The other generic subfamily of the exact solutions is obtained for g 0 = −1. First, considering this case with b > 0 in Eq. (4), condition A 2 0 > 0 imposes a constraint on the elliptic modulus, which always complies with k < 1. Further straightforward considerations demonstrate that constraint (23) does not admit sign-changing NLs, i.e., in the case of g 0 = −1 and b > 0, the exact solutions pertain to the spatially modulated nonlinearity which remains self-repulsive at all x.
Another possibility is to take g 0 = −1 and −1 < b < 0 in Eq. (4). A simple algebra demonstrates that the corresponding condition , combined with k < 1, is satisfied with any k for 2/3 < −b < 1, and is never satisfied for 0 ≤ −b < 1/2. In the intermediate case of 1/2 < −b < 2/3, the exact solutions exist in interval cf. Eq. (22). The consideration of expression (7) demonstrates that condition g 1 > 1 automatically holds in the present case, hence NL modulation function (4) is always of the sign-changing type. Thus the exact solutions in this case may, in principle, represent the ground state of the BEC with nonzero mean density, loaded into the sign-changing NL. In fact, it will be demonstrated below that solely in this case the stationary periodic solutions of the cn type may be stable (without the addition of the linear potential). It is relevant to stress that the existence of the exact cn-type solutions with g 0 = −1 is a nontrivial finding in the following sense: according to Eqs. (5), (7), (8), and (9), in the limit of b = 0 the same type of the solution, but taken with g 0 = +1, goes over into the commonly known exact unstable cnoidal solution for the GPE with g(x) ≡ +1 (the uniform self-attraction), given by expressions (12). Thus, the cn solutions for g 0 = +1 may be regarded as a continuation of this simple unstable solution, which helps to explain the fact that they are always unstable too, as shown below. On the other hand, in the case of g 0 = −1 the existence range of the exact solutions is separated, as demonstrated above, by interval 0 i.e., this subfamily of the cn solutions has no counterpart among solutions of the GPE with the constant nonlinearity coefficient, and it does not originate from that limit-in particular, it may be stable for this reason.

III. THE STABILITY OF THE CN-TYPE SOLUTIONS
The test of stability of the exact solutions is a crucially important part of the analysis. In this work, we do not aim to compute stability eigenvalues for small perturbations around stationary solutions. Instead, we focus on systematic direct simulations of perturbed solutions. Although less rigorous mathematically, this approach is closer to the description of the physical situations, where finite random perturbations are always present. The numerical integration was performed in the spatial domain with the length corresponding to 100 periods of the structure, see Eq. (6), and periodic boundary conditions. The simulations were run up to time t = 5000, in terms of the present notation. For most patterns, which have the spatial period L ≤ 5 (see below), this implies t > 100 characteristic dispersion times, estimated as T disp ∼ L 2 . In fact, in all cases the instability, if any, manifests itself at a much shorter time scale, t < 500. The simulations were extended to t = 5000 to make sure that solutions identified as a stable ones do not develop a delayed instability (still longer simulations, carried out for selected cases, completely corroborate the stability). Initial random perturbations, at the 1% amplitude level, were always sufficient to unambiguously categorize stable and unstable solutions (the application of stronger perturbations did not change the results).

A. The solutions without the linear potential
We start the presentation of the stability results by considering the cn-type solutions in the absence of the linear potential, V (x) = 0 in Eqs. (1) and (2). First of all, in all cases with g 0 = +1 and b > 0 (which suggests that the nonlinearity is self-attractive, on the average), the simulations reveal that the cn-type solutions are unstable, see a typical example in Fig. 1. In accordance with results of the above analysis, density maxima of the unperturbed solution in Fig. 1(a) coincide with minima of g(x) in Fig. 1(b), i.e., maxima of the respective pseudopotential, which readily explains the instability. Note that, although only six spatial periods of the stationary solution are shown in the figure, it is actually a representative fragment of the picture which was generated in the domain covering 100 periods, as said above. The same pertains to other figures displayed below.
A generic feature observed in Fig. 1 is that the instability sets in after evolution time t ≤ 10T disp . The instability onset may be delayed for the unperturbed solutions with a small amplitude, which is an obvious manifestation of the nonlinear nature of the instability.
In the same case of g 0 = +1 but with b < 0, unambiguously stable stationary solutions were not found either. The difference of this case from the one with b > 0 is that the respective NL function g(x) is a sign-constant one, i.e., the nonlinearity is locally attractive everywhere (which makes the modulational instability of stationary periodic solutions quite plausible).
Proceeding to the case of g 0 = −1, i.e., the model with the nonlinearity which tends to be self-repulsive on the average, no stable solutions have been found for b > 0 in Eq. (4). However, stable solutions are possible if g 0 = −1 is combined with b < 0 and appropriate values of k, which should be relatively close to k = 1 [recall that the respective exact solutions exist in the interval of 1/2 < −b < −1, and they all correspond to a sign-changing function g(x)]. An example of a stable pattern is shown in Fig. 2. In fact, this figure displays a marginally stable solution, in the sense that it is never destroyed by the originally imposed random perturbation, but the disturbance does not dissipate either, remaining trapped in the solution. Note that this example shows density maxima collocated with maxima of g(x), i.e., minima of the corresponding NL pseudopotential, which facilitates the stabilization of the pattern, and  gives it a chance to be a ground state.
Nevertheless, in the case of g 0 = −1 and b < 0 the cn solutions may be unstable if the elliptic modulus is taken too far from k = 1. For instance, Fig. 3 displays an example of a quick onset of the instability, found at the same value of b = −0.7 as in Fig. 2, but with k = 0.95 replaced by k = 0.70.
Collecting data of the systematic numerical simulations for the presently considered subfamily of the exact cn-type solutions, with g 0 = −1, 1/2 < −b < 1 and various values of k, makes it possible to identify a stability area in the plane of (b, k), as shown in Fig. 4 As said above, the stability is limited to relatively long-wave patterns, with k sufficiently close to 1, viz., k ≥ 0.9, the respective values of period (6) being L ≥ 10.12. In accordance with the above analysis, the existence region of the solutions is limited to 1/2 < −b < 1. It may be relevant to notice that the stability boundary in Fig. 4 is not monotonous, which may be understood as a consequence of the fact that the solution depends on b in a complex manner, as seen from Eqs. (7) and (8). The linear potential given by expressions (15) and (16) affects not only the shape of the solutions but also their stability. First of all, in the cases of g 0 = +1 with either sign of b, and g 0 = −1 with b > 0, well-defined stable solutions have not been found, as well as in the same cases in the absence of the linear potential, see above ("welldefined" implies that we do not count solutions with a very small amplitude, that may seem stable as the respective nonlinearity is extremely weak). A typical example of a retarded onset of the instability of a small-amplitude solution is displayed in Fig. 5, for g 0 = +1 and b = +0.5.
Similar to the situation reported above for the solutions with V (x) = 0, a vast subfamily of stable solutions is found in the case of g 0 = −1 and b < 0. Actually, the stability region is located at b > −0.73, this border being practically independent of elliptical modulus k (at b = −0.72, the exact solutions are found to be stable in interval 0.01 ≤ k ≤ 0.99, while at b = −0.74 they are unstable in the same interval). The shape of the stable solutions is quasi-harmonic at small and moderate values of k (see an example in Fig. 6), while featuring sharp peaks at k close to 1 (Fig. 7). Note that the amplitude of the solutions shown in these typical examples is not especially small, which confirms that the observed stability is not a trivial consequence of the weak nonlinearity. It is also worthy to note that density maxima of the solutions coincide with maxima of the NL modulation function, g(x), and minima of V (x), i.e., the density maxima are collocated with minima of both the nonlinear pseudopotential and linear potential. This circumstance obviously facilitates the stabilization of the periodic patterns, giving them a chance to realize the ground state.

A. Stationary solutions
Another generic class of exact solutions, which, as well as the solutions of the cn type, are represented by even functions of x, can be obtained from cn by the continuation to values of k > 1, using the well-known transformation formula, cn (x, k) = dn (kx, 1/k) .
Thus, in the absence of the linear potential, V (x) = 0, expressions (4)-(9) generate the following class of exact solutions to Eq. (2):  Figs. 1 -3, but here we display an example of the perturbed evolution of a weakly unstable cn-type solution, found in the presence of the linear periodic potential, as given by Eqs. (15) and (16). The shape of the linear potential is displayed in panel (c). The parameters are g0 = +1, b = +0.5, r = 1, g1 = −2, and k = 0.60.
where, as above, we assume normalization g 0 = ±1, the elliptic modulus takes values 0 < k < 1, and an additional normalization may be imposed, r = 1. As concerns the transition to the ordinary GPE, with g = const, Eq. (26) demonstrates that this is possible in two cases, b = 0 or under condition coinciding with Eq. (11), see above. The latter condition, if applied to Eq. (28), leads to k 2 = (b + 1) /b. Unlike the similar condition obtained for the cn-type solutions in the form of Eq. (13), the present one, obviously, cannot give appropriate values of the elliptical modulus, 0 < k 2 < 1, if b takes values for which solution (27) is nonsingular, i.e., b > −1. As for the former limit case, b = 0, expressions (27) and (29) Similar to the class of the cn-based solutions, in the presence of the linear potential the NL-modulation function is given by the same general expression as above, Eq. (26), while Eq. (28) should be dropped (actually, the latter equation is the condition for reducing V 0 to zero). Accordingly, the solution family depends on four free parameters, g 1 , r, k, b, and the sign parameter, g 0 .

B. The stability of the dn-type solutions
In the absence of the linear potential, the exact solutions given by Eqs. (26)- (29) are found to be strongly unstable, see a typical example in Fig. 8. Also unstable are the exact solutions found in the presence of the linear potential, as per Eqs. (31)- (33), in the case of g 0 = −1, on the contrary to the solutions of the cn type, and also to the odd solutions of the sn type (see below), where the potential can readily stabilize the solutions with g 0 = −1.
On the other hand, in the case of g 0 = +1, when the solutions of both the cn (see above) and sn (see below) types cannot be stable, the dn species of the periodic patterns readily demonstrates its stability, in the presence of the linear potential, as shown in Fig. 9. Note that panels (b) and (c) of this figure imply the competition between the nonlinear pseudopotential and linear potential, as minima of g(x) coincide with minima of V (x), this sets of points also coinciding with density maxima of the solution, see panel (a). A similar competition between the nonlinear and linear effective potentials was recently considered in a model of a one-dimensional photonic crystal [42]. In fact, almost all the solutions belonging to this subfamily are stable, with the exception of ones with large b or k very close to 1.

A. Solutions without the linear potential
In addition to the two classes of even stationary solutions considered above, it is possible to find a class of the NL-modulation functions, g(x), and related potentials, V (x), which are even too, but they support exact stationary In the case of g 0 = +1, expression (38) cannot give A 2 0 > 0 for b > 0. Further, a straightforward analysis demonstrates that, in the same case but for b taken from 0 < −b < 1/3, the solutions exist in interval of values of the elliptic modulus. This entire interval satisfies constraint k < 1. For 1/3 < −b < 2/3, the solutions exist for all k 2 < 1. Finally, for 2/3 < −b < 1, the existence condition is A similar analysis performed for the case of g 0 = −1 demonstrates that, with b > 0, the solutions exist for all values of k 2 < 1. Proceeding to negative b in this case, it is easy to conclude that, for 1/3 < −b < 0, the solutions with A 2 0 > 0 exist for cf. existence condition (40) for g 0 = +1, which was obtained in the same interval of b. In the adjacent interval, 1/3 < −b < 2/3, the solutions do not exist at all. Finally, in the end portion of the allowed interval of b, viz., 2/3 < −b < 1, the sn-type solutions are available in the region of cf. region (41) in the case of g 0 = +1 [the whole interval (43) complies with constraint k < 1]. The transition of g(x) given by Eq. (35) to g = const, corresponding to the ordinary GPE, may be performed by taking b = 0, or by imposing the same condition (11) as above. In the former case, the respective solution, ψ(x) = k sn(x), satisfies the GPE with g(x) ≡ −1. In the latter case, the solution is relevant for g(x) ≡ +1, reducing to b = −k 2 , µ = (1/2) 1 − 2k 2 , and ψ(x) = k √ 1 − k 2 sn(x)/dn(x), cf. similar limit solution (14) of the cn type.

B. Exact sn solutions with a linear potential
Another family of the sn solutions can be found, adding to Eq. (2) the linear potential in the following form, Then, ψ(x) is again given Eq. (36), but Eqs. (38) and (39) are replaced by while Eq. (37) is dropped.
In the limit of b → 0 and g 1 → 0, the present solution makes sense for g 0 = −1. In this limit, we may again fix r = 1, hence Eqs. (36), (46) and (47) yield Expressions (48) provide for a solution to an equation with the constant nonlinearity coefficient, These particular results, obtained for g(x) ≡ −1, reproduce findings reported in Ref. [12].

C. The stability
In the absence of the linear potential, nearly all the solutions of the sn type are unstable, as illustrated by a typical example in Fig. 10. Note that, in this case, maxima of the density coincide with minima of g(x), thus pushing the solutions towards instability. The only stable solutions could be found for g 0 = −1 and very small |b|, when, as said above, the exact solutions are close to the commonly known stable exact solution, ψ(x) = k sn(x), to the GPE with g(x) ≡ −1.
While, unlike the cn solutions, all patterns of the sn type tend to be unstable without the linear potential, in the presence of the potential they may be readily stabilized, in the case of g 0 = −1 and b < 0, which resembles the respective property of the cn-type patterns, see Fig. 4. An example of the stable sn structure is shown, for this case, in Fig. 11. As is typical for other types of stable patterns, in this figure we observe that maxima of the density coincide with maxima of g(x) and minima of V (x), cf. Figs. 6 and 7.
A stability region, summarizing the results of many runs of the simulations, can also be identified in this case, as shown in Fig. 12, cf. the stability diagram shown in Fig. 4 for the cn solutions (without the linear potential). The dashed curve with arrows designates the existence border for the solutions of this type, as determined by condition A 2 > 0. It follows from Eq. (46), that, in the case shown in Fig. 12, this existence condition amounts to k 2 < |b|.

A. Stationary solutions for bright solitons
In the limit of k → 1, all the above types of the exact periodic solutions go over into localized structures (solitons) pinned by the corresponding localized pseudopotential, possibly in the combination with the localized (trapping) linear potential. In this limit, the solutions of both the cn and dn types, as given by Eqs. (5), (8), (9), or (27), (29), (30), reduce to a common expression for the bright soliton, with µ sol = −1/2. The nonlinearity-modulation profile which supports this solution, in the absence of the linear potential, is given by Eqs. (4) and (7), in which cn(x) is replaced by sech x, i.e., and coefficient g 1 is replaced by its limit form for k = 1, i.e., as follows form Eq. (7),  The soliton may be stable only if g(x) given by Eq. (51) has a maximum at x = 0 [for a minimum of g(x) at x = 0, an obvious variational consideration of energy (3) immediately predicts an instability of the soliton placed at the maximum of the respective pseudopotential]. The general condition for the existence of the maximum of g(x) at x = 0 was actually obtained above, g 1 > g 0 b. In the present case, with g 1 taken as per Eq. (52), it amounts to −1 < b < 0, which is, thus, a necessary condition for the stability of the soliton. Since solution (50) exists only for g 0 (1 + 2b) > 0, we finally conclude that the solitons may be stable only in the case of In fact, the nonlinearity is sign-changing in this case, because g(x = ±∞) = g 0 = −1, while g(x = 0) = (1/2)(2 − |b|)/(2|b| − 1) > 0. Lastly, the norm of soliton (50) is A more general bright-soliton solution can be obtained by setting k = 1 in the cn or dn solutions including the linear potential, i.e., Eqs. (4), (5), (15), (16), (17), and (18) [recall that r cannot be set equal to 1 in this case, and g 1 is a free parameter, while Eq. (7) is irrelevant]. This leads to µ = 1 − r 2 /2 and The existence condition for solution (55) is which is different from that at which soliton (50) exists. Further, the nonlinearity-modulation function supporting solution (55) is given by the same expression (51) as above, with the difference that g 1 is a free parameter in it, and the respective form of the linear potential is Potential (57) stabilizes the soliton if it represents a potential well at x = 0, which means V 0 < b. A straightforward analysis of expression (58), taking into regard Eq. (56), demonstrates that the latter condition is met for

B. The stability of bright solitons
Systematic simulations of Eq. (1) with V (x) = 0, starting from the initial condition corresponding to a perturbed bright soliton, demonstrate that the solitons are indeed stable in region (53), where this should be expected, as argued above. An example of the stable soliton is shown in Fig. 13. For values of b approaching −1, see Eq. (53), the initially imposed random perturbations gradually grow, which, however, is interpreted as the growing sensitivity of the solution to the random perturbations, rather than as a true instability. A technical definition of the "low sensitivity" to the random noise may be adopted in the form of the requirement that the relative amplitude of the perturbations must remain below 3%, after the evolution time t = 5000, if the initial perturbation amplitude was 1%, as fixed above. In particular, for the same parameters as in Fig. 13, r = 1 and g 1 = 2.4, the so defined sensitivity border is found at b = −0.75. For instance, the final perturbation amplitude is 2.82% and 3.36%, at b = −0.7 and −0.8, respectively.
In the presence of the linear potential, stable solitons can be found both for g 0 = −1, as above, and for g 0 = +1, see Figs. 14-16. For instance, fixing g 0 = +1, r = 1 and g 1 = −2, as in Figs. 14 and 15, we find that the solitons are stable for b < 6, and unstable for b ≥ 6 [the soliton is unstable at b = 6, as shown in Fig. 15, but still stable-at least, up to t = 5000 -at b = 5.9 (not shown here)]. In fact, the transition from the stability to instability, in the form of a spontaneous escape of the trapped soliton, as observed in Fig. 15, is explained by the competition between the self-repulsive localized nonlinear pseudopotential, and the attractive linear trapping potential, see panels (b) and (c) in Figs. 14 and 15. Note that parameters corresponding to both figures satisfy condition (59), under which the linear potential is the trapping one.
Lastly, the combination of g 0 = −1 with the linear potential makes both the nonlinear pseudopotential and linear potential attractive, hence the solitons are definitely stable in this case, see an example in Fig. 16. For the parameters chosen as in this figure, g 0 = −1, r = 1, g 1 = 2, the above-mentioned technical definition of the sensitivity of the trapped soliton to the random noise (the perturbation amplitude at the level of 3% by t = 5000) sets the "sensitivity border" at b = −0.65, the solitons being "hypersensitive" at −0.65 < b < −1.

C. Dark solitons
The limit form of the exact solutions of the sn type for k = 1 are dark solitons. In the general case, when the respective linear potential (44) is included, the dark-soliton solution can be obtained from Eqs. (36) and (46) in the  Fig. 14, but for b = 6. In this case, the nonlinear self-repulsive pseudopotential turns out to be stronger than the linear trapping potential, hence the bright soliton is unstable, being spontaneously expelled from the bound state.
following form: The nonlinearity-modulation function and linear potential which support this dark soliton are obtained from Eqs. (35) and (44) by setting k = 1 in them: g(x) = (g 0 + g 1 ) − g 1 sech 2 (rx) V (x) = − (U 0 ) dark−sol tanh 2 (rx) The stability of the dark solitons will be considered elsewhere.

VII. CONCLUSIONS
This work reports three large classes of exact periodic solutions for the one-dimensional GPE (Gross-Pitaevskii equation) with the pseudopotential represented by the periodic modulation of the nonlinearity coefficient, which can be created in the BEC by means of the Feshbach-resonance technique. In the general case, the model includes a periodic linear-lattice potential too. The modulation function of the corresponding NL (nonlinear lattice), the linear potential (if any), and the exact solutions are expressed in terms of the Jacobi's elliptic functions of three types-cn, dn, and sn, which give rise to the three different classes of the solutions. In the absence of the linear potential, the solutions depend on two continuous free parameters, b and k, and the sign parameter, g 0 . If the linear potential is included, full solution families feature two additional continuous parameters, g 1 and r. The stability of the exact solutions was tested by means of systematic direct simulations of the perturbed evolution, except for dark solitons, whose stability will be considered separately.
Subfamilies of the pseudopotentials and respective exact solutions have been identified with both sign-changing and sign-constant nonlinearity-modulation functions. Density maxima of the solutions may coincide with minima or maxima of the periodic pseudopotential; in the former case, the exact solutions may represent the ground state of the BEC. In the absence of the linear potential, only the solutions with the density maxima collocated with minima of the pseudopotential may be stable (although this condition alone is not sufficient for the stability, in most cases). On the other hand, the addition of the trapping linear potential may stabilize solutions whose density maxima coincide with maxima of the nonlinear pseudopotential.
In the limit of the infinite period, the exact solutions reduce to solitons. The solitons are stable in a large part of parameter regions where they exist. In particular, they may be stable if located at a maximum of the nonlinear pseudopotential, provided that the competing trapping linear potential is strong enough.
This work calls for continuation in other directions. In the models considered here (with both the periodic and localized pseudopotentials), it is interesting to find, by means of numerical methods, a full set of nonlinear states, which should reveal how the exact solutions are embedded into the full family, and exactly identify the respective ground state. Further, it is desirable to investigate the dynamical stability of the solutions in a rigorous form, through the computation of eigenvalues for modes of small perturbations. Another interesting issue is a possibility of finding similar families of exact solutions for systems of two or three coupled GPEs describing multi-component BECs. In fact, the methods used in earlier works [20,22,60] for obtaining families of exact solutions in coupled NLSEs can be used in the latter context.