of Eq. (3)) [56]. Particularly, one should expect that in the case of a purely harmonic trap the soliton oscillates with a frequency ω$_{ds}$ > Ω $_{/}$√$_{2, which can be found numerically in}$ the absence of the optical lattice (note that Ω $_{/}$√$_{2 is the characteristic value of the soliton}$ oscillation frequency in a 1D harmonic trap of strength Ω in the TF limit [21–29]). For a more detailed discussion see Sec. IV A. Moreover, for sufficiently small optical lattice wavenumbers (such as the chosen value of k = π/ 5 . 37 µm − $^{1}$), and large chemical potentials (TF limit) the effective potential of Eq.(13) can be simplified to the following form: We complete this Section by mentioning that the case of two (or more) spatially well sep- arated solitons can also be treated analytically, taking into account the repulsive interaction potential between two solitons. In the case of two spatially separated solitons located at z$_{1}$ and z$_{2}$ , on top of a constant background with density n$_{0}$ , the interaction potential takes the form [15], Let us now proceed by investigating the existence, the stability and possible bifurcations of the lowest macroscopically excited states of Eq. (3). We fix the parameter values as follows: Ω = 0 . 1, k = π/ 5 . 37 µ m − 1 and V$_{0}$ = 1 . 16 × 10 − 12 eV; for such a choice, the four lowest eigenvalues of the operator 1 + ˆ H (corresponding to the linear problem) are found to be: and the energy of the first excited state is almost degenerate with that of the ground state, i.e., ω$_{0}$ ∼ ω$_{1}$ . Fig. 1 shows the number of atoms as a function of the chemical potential for different branches of solutions. The insets show the spatial density profiles for µ = 3 × 10 − 12 eV, with the units of the horizontal and vertical axes being given by µm and 1 /µm , respectively. Branch (1) corresponds to the states with the largest number of atoms, for fixed chemical potential. The respective states are symmetric and have no nodes. Continuation of this ω$_{0}$ = 2 . 138 × 10 − $^{12}$eV , ω$_{1}$ = 2 . 142 × 10 − $^{12}$eV , ω$_{2}$ = 2 . 634 × 10 − $^{12}$eV , ω$_{3}$ = 2 . 700 × 10 − $^{12}$eV , (16) V $^{int}$( z$_{1}$, z$_{2}$ ) = n$_{0}$ sinh $^{2}$( √ n$_{0}$ ( z$_{2}$ − z$_{1}$ )) . (15) V$_{eff}$ ( z ) = ω 2 $_{ds}$z 2 + 1 2 V$_{0}$ cos(2 kz ) . (14) 9 III. BIFURCATION AND BDG ANALYSIS