small-amplitude approximation) [61]. Nevertheless, exact analytical dark soliton solutions of Eq. (3) can be found in the weakly-interacting limit (4 | ψ | 2 ≪ 1), and in the absence of the external potential: in this case, Eq. (3) is reduced to the completely integrable cubic NLS model, which possesses single- and multiple-dark-soliton solutions on top of a background with constant density n = n$_{0}$ = µ (with µ being the chemical potential). A single dark soliton solution has the form [62], where η = √ n$_{0}$B [ z − z$_{0}$ ( t )], z$_{0}$ ( t ) is the soliton center, the parameter B ≡ √$_{1}$ − ν 2 sets the soliton depth given by √ n$_{0}$B , while the parameter ν sets the soliton velocity, given by dz$_{0}$/dt = √ n$_{0}$ν . Note that for ν = 0 the dark soliton becomes a black soliton (alias a stationary kink), with a π phase jump across its density minimum. Aside from the single- dark-soliton, multiple dark soliton solutions of the cubic NLS equation are also available [62–64]. In the simplest case of a two-soliton solution, with the two solitons moving with opposite velocities, ν$_{1}$ = − ν$_{2}$ = ν , the wave function can be expressed as [64] (see also Ref. [65]): where F = 2( n$_{0}$ − 2 n$_{min}$ ) cosh( T ) − 2 n$_{0}$ν cosh( Z ) + i sinh( T ), G = 2 √ n$_{0}$ cosh( T ) + 2 √ n$_{min}$ cosh( Z ), while Z = 2 √ n$_{0}$Bz , T = 2 √ n$_{min}$ ( n$_{0}$ − n$_{min}$ ) t , and n$_{min}$ = n$_{0}$ − n$_{0}$B 2 = n$_{0}$ν 2 is the minimum density (i.e., the density at the center of each soliton). Generally, the single dark soliton, as well as all higher-order dark soliton states, can be obtained in a stationary form from the non-interacting (linear) limit of Eq. (3), corresponding to N → 0 [66, 67]. In this limit, Eq. (3) is reduced to a linear Schr¨ odinger equation for a confined single-particle state, namely the equation of the quantum harmonic oscillator, together with the contribution of the optical lattice. However, due to the presence of the harmonic trap, the problem is characterized by discrete energy levels and corresponding localized eigenmodes. In the weakly-interacting limit, where Eq. (3) is reduced to the cubic NLS equation, all these eigenmodes exist for the nonlinear problem as well [66, 67], describing an analytical continuation of the linear modes to a set of nonlinear stationary states. Recent analysis and numerical results [68] (see also Ref. [15]) suggest that there are no solutions of Eq. (3) without a linear counterpart, at least in the case of a purely harmonic potential. However, as we will show below, the presence of the optical lattice in the model results in 7 ψ ( z, t ) = F ( z, t ) G ( z, t ) exp( − iµt ) , (9) ψ ( z, t ) = √ n$_{0}$ [ iν + B tanh( η )] exp( − iµt ) , (8)