text:
# Verifying the Kugo–Ojima Confinement Criterion in Landau Gauge Yang–Mills theory

\[

## Abstract

Expanding the Landau gauge gluon and ghost two–point functions in a power series we investigate their infrared behavior. The corresponding powers are constrained through the ghost Dyson–Schwinger equation by exploiting multiplicative renormalizability. Without recourse to any specific truncation we demonstrate that the infrared powers of the gluon and ghost propagators are uniquely related to each other. Constraints for these powers are derived, and the resulting infrared enhancement of the ghost propagator signals that the Kugo–Ojima confinement criterion is fulfilled in Landau gauge Yang–Mills theory.

preprint:   UNITU-THEP-08/01 hep-ph/0102332

\]

Modifications to the standard framework of quantum field theory are necessary to accommodate the confinement phenomenon in QCD. Here either relaxing the principle of locality or abandoning the positivity of the representation space (or both) is required. As quantum gauge field theories in covariant gauges demand indefinite metric spaces, and as locality is of fundamental importance, it is suggestive to relate confinement to the violation of positivity in Yang–Mills theories: No colored states should be present in the positive definite space of physical states defined by some suitable condition maintaining physical $`S`$-matrix unitarity. Exploiting the BRS symmetry of the gauge fixed action a nilpotent BRS-charge $`𝒬_{}`$ in an indefinite metric space $`𝒱`$ can be constructed. The subspace $`𝒱_\text{p}=\text{Ker}Q_B`$ is defined by those states which are annihilated by the BRS charge $`𝒬_{}`$. Positivity is then proved for physical states  (see also refs. ) which are given by the cohomology $`(Q_B,𝒱)=\text{Ker}Q_B/\text{Im}Q_B`$, the covariant space of equivalence classes of BRS-closed modulo BRS-exact states.

Longitudinal and timelike gluons form BRS quartets together with ghosts and are thus unphysical. At the same time the global symmetry $`J_{\mu ,\nu }^a`$ corresponding to gauge transformations generated by $`\theta ^a(x)=a_\mu ^ax^\mu `$ is spontaneously broken quite analogous to the displacement symmetry in QED. The identification of BRS-singlet states with color singlets is non-perturbatively possible if this global symmetry is restored . The condition for a dynamical restoration of this symmetry is the Kugo–Ojima confinement criterion which is a requirement on the infrared behavior of a four-point Green’s function. In Landau gauge, a sufficient condition for this criterion is that the nonperturbative ghost propagator is more singular than a massless pole in the infrared :

$$D_G^{ab}(p)=\delta ^{ab}\frac{G(p^2)}{p^2},\text{with}G(p^2)\stackrel{p^20}{}\mathrm{}.$$
(1)
This mechanism is correlated to the derivation  of the Oehme–Zimmermann superconvergence relations  from Ward–Takahashi identities. These superconvergence relations formalize a long known contradiction between asymptotic freedom and the positivity of the spectral density for transverse gluons in the covariant gauge.

The violation of positivity of transverse gluons is unambiguously established by a variety of independent nonperturbative studies of the gluon propagator . Furthermore, recent lattice calculations investigate the Kugo–Ojima criterion directly . As infrared singularities are anticipated, also non-perturbative continuum methods, besides lattice calculations, are necessary to investigate this and related pictures of confinement in more detail.

The Green’s functions of a quantum field theory obey the Dyson–Schwinger equations. For realistic theories this infinite tower of integral equations cannot be solved directly, and truncations are necessary to obtain (numerical) solutions. The ones for QCD propagators can then be used to formulate a successful hadron physics phenomenology, see e.g. the recent reviews .

The results of two recent truncation schemes  suggest that the Kugo–Ojima confinement criterion for the ghost propagator and the violation of positivity for transverse gluon correlations are fulfilled. Obviously, the question arises whether these results are sensitive to the employed truncation schemes.

In this letter we will exploit the general structure of the simplest QCD Dyson–Schwinger equation, the one for the ghost two–point function, to relate the infrared behavior of gluon and ghost propagators. Thereby the Kugo–Ojima confinement criterion will be verified by purely analytical methods. Using the rather unique properties of the ghost-gluon vertex in the Landau gauge and assuming the applicability of a powerlaw ansatz in the infrared the ghost Dyson–Schwinger equation admits a completely general qualitative analysis. Instead of truncating the system, the general forms for the unknown functions constrained by various consistency arguments, notably multiplicative renormalizability, are input into the equation under an asymptotic expansion in the infrared momentum scale, c.f. the method described in ref. .

We start by providing some basic notations. The general form of the gluon propagator in Landau gauge is purely transverse and is given by

$$D_{\mu \nu }^{ab}(p)=\delta ^{ab}\left(g_{\mu \nu }\frac{p_\mu p_\nu }{p^2}\right)\frac{Z(p^2)}{p^2},$$
(2)
which defines the gluon renormalization function $`Z(p^2)`$. The ghost renormalization function $`G(p^2)`$ is established in eq. (1). We define the ghost-gluon vertex

$$\mathrm{\Gamma }_\mu ^{abc}(p_1,p_2,p_3)=igf^{abc}\mathrm{\Gamma }_\mu (p_1,p_2,p_3)$$
(3)
with antighost ($`p_1`$), ghost ($`p_2`$) and gluon ($`p_3`$) momentum incoming, i.e. $`p_1+p_2+p_3=0`$. At tree-level the dressing functions reduce to

$$Z(p^2)=G(p^2)=1,\mathrm{\Gamma }_\mu (p_1,p_2,p_3)=p_{1\mu }.$$
(4)
These renormalized dressing functions are related to the bare ones by the usual renormalization constants $`Z_3`$ (gluon wave function), $`\stackrel{~}{Z}_3`$ (ghost wave function) and $`\stackrel{~}{Z}_1`$ (ghost-gluon vertex). In Landau gauge one has $`\stackrel{~}{Z}_1=1`$ , i.e., the ghost-gluon vertex has not to be renormalized .

The renormalized ghost dressing function $`G(p^2)`$ is a dimensionless function of a single variable. In pure Yang–Mills theory there are only two scales present: the external scale $`p^2`$ and the renormalization scale $`\mu ^2`$. Therefore $`G`$ is a function of the ratio $`p^2/\mu ^2`$ only. As the external scale vanishes this function can be written as an asymptotic expansion. Renormalizability implies that it must also depend on the renormalized coupling $`g`$. This motivates the power series ansatz

$$G(p^2;\mu ^2)=\underset{n}{}d_n\left(\frac{p^2}{\mu ^2}\right)^{\delta _n},$$
(5)
where the coefficients $`d_n`$ and the exponents $`\delta _n`$ are, at least in principle, functions of $`g`$. We will see below that the exponents $`\delta _n`$ will be independent of $`g`$. Furthermore, as will be shown later on, the series (5) has a lowest power $`\delta _0`$ such that as the external scale $`p^2`$ vanishes the ghost dressing function $`G(p^2)`$ is well approximated by a powerlaw. Multiplicative renormalizability allows to derive the renormalization group (RG) equation

$$\left(\mu \frac{}{\mu }+\beta (g)\frac{}{g}2\gamma _\mathrm{G}\right)G(p^2;\mu ^2)=0$$
(6)
with $`\beta (g)`$ being the beta function, and

$$\gamma _\mathrm{G}=\mu \frac{\mathrm{log}\stackrel{~}{Z}_3^{\frac{1}{2}}}{\mu }$$
(7)
is related to the ghost anomalous dimension. Inserting the series (5) into eq. (6), and exploiting that $`p^2`$ is a free (though small) parameter one obtains

$`d_n\delta _n+\beta (g)\left({\displaystyle \frac{d_n}{g}}+d_n\mathrm{log}({\displaystyle \frac{p^2}{\mu ^2}}){\displaystyle \frac{\delta _n}{g}}\right)+d_n\gamma _\mathrm{G}=0.`$ (8)
Since the logarithm has no matching factor the exponents $`\delta _n`$ must be independent of the renormalized coupling ($`\delta _n/g=0`$). Using $`\beta (g)=g(2\gamma _\mathrm{G}+\gamma _\mathrm{A})`$ with $`\gamma _\mathrm{A}=\mu \mathrm{log}Z_3^{\frac{1}{2}}/\mu `$ the resulting differential equation for the coefficients $`d_n`$ has the general solution

$$d_n=\mathrm{const}.g^{\frac{2(\delta _n+\gamma _\mathrm{G})}{2\gamma _\mathrm{G}+\gamma _\mathrm{A}}}.$$
(9)
Expanding the renormalized gluon dressing function

$$Z(p^2;\mu ^2)=\underset{n}{}e_n\left(\frac{p^2}{\mu ^2}\right)^{ϵ_n}$$
(10)
one arrives in a completely analogous way at

$$e_n=\mathrm{const}.g^{\frac{2(ϵ_n+\gamma _\mathrm{G})}{2\gamma _\mathrm{G}+\gamma _\mathrm{A}}}.$$
(11)
In the renormalization process we have to specify the corresponding normalization conditions. We require that the unrenormalized functions reduce to a constant at asymptotically large momenta. Exploiting additionally multiplicative renormalizability and using eq. (5) one obtains the formal relation

$$\stackrel{~}{Z}_3(\mu ^2,\mathrm{\Lambda }^2)\left(\underset{n}{}d_n\left(\frac{\mathrm{\Lambda }^2}{\mu ^2}\right)^{\delta _n}\right)^1.$$
(12)
This equation is only well-defined if the corresponding series makes sense, e.g. if it is resumable by some method. This is the main assumption throughout the presented derivation. One consistency check is provided by inserting the general solution for the coefficients $`d_n`$ (9) and rewriting the renormalized coupling $`g`$ in terms of the scale $`\mu `$. This yields

$$\stackrel{~}{Z}_3(\mu ^2,\mathrm{\Lambda }^2)\left(\frac{\mathrm{\Lambda }^2}{\mu ^2}\right)^{\gamma _\mathrm{G}}$$
(13)
which is precisely the correct dependence to agree with the definition of $`\gamma _\mathrm{G}`$. A completely analogous argument applies to the gluon renormalization constant $`Z_3`$.

The tensor structure of the ghost-gluon vertex reads

$$\mathrm{\Gamma }_\mu (p_1,p_2,p_3)=p_{1\mu }A(p_i^2)+p_{3\mu }B(p_i^2).$$
(14)
In Landau gauge the ghost-gluon vertex becomes bare when the ghost momentum $`p_2`$ vanishes :

$$\mathrm{\Gamma }_\mu (p_1,p_2,p_3)=p_{1\mu }\left(AB\right)p_{2\mu }B\stackrel{p_20}{}p_{1\mu }.$$
(15)
Thus, as $`p_20`$ the functions $`A`$ and $`B`$ cannot be singular. This is in accordance with the fact that the renormalization constant $`\stackrel{~}{Z}_1=1`$.

The RG equation for the vertex function can be split into two parts, one for each of the tensor components. It is only necessary to consider the equation for $`A`$, and since $`\stackrel{~}{Z}_1=1`$, it is identical to the RG equation for the coupling. Writing the vertex function as a series has the complication of several momentum scales. For our considerations it is sufficient to simply extract the scale and coupling dependence

$$A(p_i^2;\mu ^2)=\underset{k}{}a_k(p_1^2,p_2^2)\lambda _k\left(\frac{p_3^2}{\mu ^2}\right)^{\kappa _k},$$
(16)
with

$$\lambda _k=\left(g^2\right)^{\frac{\kappa _k}{\gamma _\mathrm{A}+2\gamma _\mathrm{G}}}.$$
(17)
We must renormalize the ghost-gluon vertex at the point $`p_2=0`$ since this corresponds to a zero finite renormalization from $`\stackrel{~}{Z}_1=1`$. Knowing the $`p_20`$ limit of the full vertex does not constrain $`A`$. (When $`p_2=0`$ one only has $`AB=1`$.) With this given limit and the additional requirement that the ghost equation should provide a consistent solution one infers that

$`A(p_1^2=p_3^2=\mathrm{\Lambda }^2,p_2^2=0;\mu ^2)={\displaystyle \underset{k}{}}a_k\lambda _k\left({\displaystyle \frac{\mathrm{\Lambda }^2}{\mu ^2}}\right)^{\kappa _k}`$ (18)
is a finite positive constant. Asymptotic freedom requires that $`\gamma _\mathrm{A}+2\gamma _\mathrm{G}>0`$. In the unrenormalized theory the vertex should be bare when the unrenormalized coupling vanishes because $`\stackrel{~}{Z}_1=1`$. Therefore $`\kappa _k0`$ with the highest value $`\kappa _0=0`$.

The ultraviolet region of the ghost Dyson–Schwinger equation corresponds to a small (albeit non-vanishing) ghost momentum $`p_2`$ and all other momenta large. For the following discussions it is helpful to write

$`A(p_i^2;\mu ^2)_{p_20}={\displaystyle \underset{k,m}{}}b_k(p_1^2)\lambda _k\left({\displaystyle \frac{p_3^2}{\mu ^2}}\right)^{\kappa _k}\left({\displaystyle \frac{p_2^2}{p_3^2}}\right)^{\rho _m}`$ (19)
where the $`\rho _m0`$ and $`\rho _0=0`$.

The full Euclidean ghost Dyson–Schwinger equation can be written after having performed the two trivial angular integrals as

$`G^1(x)=\stackrel{~}{Z}_3`$ (20)
$`{\displaystyle \frac{C_Ag^2}{16\pi ^2}}{\displaystyle _0^{\mathrm{\Lambda }^2}}𝑑yZ(y)`$ $`{\displaystyle \frac{2}{\pi }}{\displaystyle _0^\pi }𝑑\theta \mathrm{sin}^4\theta {\displaystyle \frac{G(z)}{z}}A(z,x,y),`$ (21)
where $`\theta `$ is the angle between the two internal loop momenta corresponding to $`y=q^2`$ and $`z=(pq)^2`$ with $`x=p^2`$. $`C_A=N_c`$ is the resulting color factor. When using the expansions (5) and (10) the angular integration splits the radial integral into two parts, the lower part $`I_{gh}`$ ($`0x`$) and the upper part $`J_{gh}`$ ($`x\mathrm{\Lambda }^2`$). The expression derived from the upper limit will be some combination of factors with the ratio $`x/\mathrm{\Lambda }^2`$ to some power. The effect of an exact angular integration is to generate a hypergeometric function whose argument is $`x/y`$ and which can be expanded in a convergent series with positive integer powers of $`x/y`$ starting with unity. Thus, the $`z`$-dependence of the vertex function is unimportant.

The most important aspect is that the ghost-gluon vertex reduces to its tree-level form as the in-ghost momentum vanishes. In the ultraviolet part of the integral this is precisely the case: the in-ghost momentum is the infrared external scale $`x`$. Therefore the vertex function can be written in the form of eq. (19).

Using the series for $`G`$ and $`Z`$, it is possible to write the ghost equation with the most general vertex function as

$`\left({\displaystyle \underset{n}{}}d_n\left({\displaystyle \frac{x}{\mu ^2}}\right)^{\delta _n}\right)^1=\stackrel{~}{Z}_3`$ (22)
$`{\displaystyle \frac{C_Ag^2}{16\pi ^2}}{\displaystyle \underset{j,k}{}}e_jd_k\left({\displaystyle \frac{x}{\mu ^2}}\right)^{ϵ_j+\delta _k}`$ $`\left(I_{gh}(ϵ_j,\delta _k)+J_{gh}(ϵ_j,\delta _k)\right).`$ (23)
where $`I_{gh}`$ and $`J_{gh}`$ represent the lower and upper regions of the integral, respectively. The integrals $`I_{gh}`$ must contain factors characterized by the $`\kappa _k`$ but there are no restrictions on its form except that the vertex must have positive powers of the coupling. As the integral $`I_{gh}`$ is very difficult to compute we will not be able to proceed by studying its general properties.

It is frequently asserted that only the infrared integrals of Dyson–Schwinger equations can contribute to the renormalized propagator functions on the left hand side because this is the only fully renormalizable possibility. Indeed, one can show that upper integral $`J_{gh}`$ will not contribute directly to the propagator functions . It is the renormalization constant $`\stackrel{~}{Z}_3`$ that provides the key to the infrared behavior of the propagators, and only the ultraviolet integral can yield the ultraviolet divergent part of the ghost equation. First we note that $`J_{gh}(ϵ_j,\delta _k)`$ can be expanded as

$`J_{gh}={\displaystyle \underset{l,m,q}{}}J_{j,k,l,m,q}\left({\displaystyle \frac{\mathrm{\Lambda }^2}{x}}\right)^{ϵ_j+\delta _k\rho _mq}\lambda _l\left({\displaystyle \frac{\mathrm{\Lambda }^2}{\mu ^2}}\right)^{\kappa _l}`$ $`,`$ (24)
where $`\rho _m0`$ represents the kinematical content of the vertex function, and $`q=0,1,\mathrm{}`$ denotes the order of the angular integration series. As $`\stackrel{~}{Z}_3`$ has to be independent of $`x`$ one deduces from eqs. (22,24) that only $`\rho _0=q=0`$ will contribute. All other terms involve $`x/\mathrm{\Lambda }^2`$ to a positive power and subsequently vanish:

$`\stackrel{~}{Z}_3(\mu ^2,\mathrm{\Lambda }^2)`$ $`=`$ (25)
$`{\displaystyle \frac{C_Ag^2}{16\pi ^2}}`$ $`{\displaystyle \underset{j,k,l}{}}e_jd_kJ_{j,k,l,0,0}\lambda _l\left({\displaystyle \frac{\mathrm{\Lambda }^2}{\mu ^2}}\right)^{ϵ_j+\delta _k+\kappa _l}.`$ (26)
Using that expression (18) is a constant, employing the expansion (12) and multiplying through with the left hand side yields

$$\frac{C_Ag^2}{16\pi ^2}\underset{i,j,k}{}e_jd_id_kJ_{0,0,0}(ϵ_j,\delta _k)\left(\frac{\mathrm{\Lambda }^2}{\mu ^2}\right)^{ϵ_j+\delta _i+\delta _k}=\mathrm{const}.$$
(27)
Since $`ϵ_0`$ and $`\delta _0`$ are the unique lowest powers, and the term with $`i=j=k=0`$ has nothing to compete with it must be responsible for the non-vanishing part of the left hand side of this equation. Therefore it must be independent of $`\mathrm{\Lambda }^2/\mu ^2`$, and we obtain

$$ϵ_0=2\delta _0,$$
(28)
i.e. we have shown that the leading infrared powers of the gluon and the ghost propagator are uniquely related. Note that this result is valid even if corresponding series are not convergent. We have simply matched the coefficients of an asymptotic series. Therefore, our underlying assumption is that the employed series have to exist in the sense of being, at least, asymptotic ones.

Eq. (27) constrains the values of the related leading powers. To do this, one needs to explicitly evaluate the integrals $`J_{j,k,0,0,0}`$ which are readily calculated from the full ghost equation with a bare vertex stripping off the terms that depend on the external scale $`x`$. Concentrating on the leading powers we obtain:

$$\frac{3}{4}\frac{C_Ag^2}{16\pi ^2}e_0d_0^2\frac{1}{\delta _0}>0.$$
(29)
To yield a positive ghost renormalization function (i.e. to keep the negative-norm property of the ghost states) the corresponding leading order coefficients $`e_0`$ and $`d_0`$ must be positive. As the coupling is positive this yields:

$$\delta _0<0.$$
(30)
This relation implies that the Kugo–Ojima confinement criterion, see eq. (1), is fulfilled.

A lower bound on $`\delta _0`$ can be derived from the lower bound of the integral in the ghost equation. As this bound is considerably weaker as the one following from a generalized masslessness condition (see eq. (192) in ref. ) $`\delta _0>1`$, we will not present its derivation here. This implies $`2>ϵ_0>0`$: The gluon propagator is suppressed in the infrared as compared to its tree-level form.

In summary, we have investigated the general structure of the ghost Dyson–Schwinger equation in Landau gauge. A general series form for the propagator and vertex functions has been employed. No specialized ansatz for the ghost-gluon vertex has been used. Care has been taken not to use approximations that destroy the generality of the results. Multiplicative renormalizability, and especially the non-renormalization of the ghost-gluon vertex, have placed the constraints to relate uniquely the leading infrared behavior of the gluon and the ghost propagator.

We have found an infrared enhanced ghost propagator thereby verifying the Kugo–Ojima confinement criterion  directly. Due to the relation (28) this corresponds to an infrared suppressed gluon propagator. For a non-integer leading infrared power this implies that positivity is violated and the Oehme–Zimmermann superconvergence relation  is fulfilled non-perturbatively. Considering only the ghost equation as in this letter an infrared constant gluon propagator cannot be excluded. To investigate this question one has to consider additionally the gluon Dyson–Schwinger equation .

Acknowledgements
We want to express our gratitude to M.R. Pennington who has accompanied the whole research presented here with helpful discussions. We thank C. Fischer, S. M. Schmidt and L. von Smekal for their valuable remarks.

P. W. acknowledges the University of Durham and the CSSM, University of Adelaide, where major parts of this research have been done.

This research has been supported by DFG under contract Al 279/3-3.


format: text

text:
# Stability of trapped Bose-Einstein condensates To appear in Physical Review A

## I Introduction

In trapped Bose-Einstein condensates(BEC), it is well-known the occurrence of collapse when the two-body interaction is attractive and if the number of atoms $`N`$ exceeds a critical value $`N_c`$, as in the case of atomic condensates with $`{}_{}{}^{7}\text{Li}`$ . In this case, experiments with attractive two-body interaction have been performed , with results consistent with the limitation in the number of atoms and with the grow and collapse scenario. As the non-linear terms presented in the non-linear Schrödinger equation (NLSE) (cubic, quintic, and so on) are due to the first terms of an expansion of the effective many-body interaction in the mean-field approximation , a study of stability of the equation that describes BEC is necessary, considering the observed few-body interactions that occur in the atomic gas. When the overall sign of the effective many-body interaction is negative, we have the conditions for the occurrence of collapse, with the stability of the condensate being restricted by the number of atoms.

We should stress that the criteria for stability usually involves one parameter: the number of atoms. Here, in order to introduce a criterion for stability of a trapped condensate with two and three-body effective interactions, we introduce another parameter: a chirp, that is related to the initial focusing (defocusing) of the cloud. Without trap, a chirp parameter was previously introduced in Ref.  to study the stability of the NLSE; and with trap, the stability has been studied in Refs. . One should also note that an oscillating condensate always has a chirped wave function, where the chirp is proportional to the time variation of the width, which is a periodic function of time. Thus, a real condensate, for example in optical traps, always has nonzero initial chirp, that should be taken into account.

Collapses for two and three dimension ($`2D`$ and $`3D`$) geometries of the NLSE with trap potential has been investigated in Refs. , where the corresponding critical numbers were obtained . These works are extensions of the Zakharov and co-workers theory of collapse  to the case when an additional trap potential is included. The trap potential helps to prevent the collapsing process, which cannot be avoided if the number of atoms exceeds the critical limit. In Refs. , considering the moments method in the analysis of the equation for the average squared value of the width, $`a^2`$, it was emphasized the influence of the corresponding initial condition on the collapse.

Variational approaches have also been used by many authors for the analysis of stability of NLSE, and proved to be useful in many aspects (see Ref. , for example). The variational approach (VA), using a simple Gaussian ansatz, gives a reasonable description of the conditions for the collapse and also an approximate value for the critical number of atoms. It is interesting to apply the time dependent variational approach to describe the dynamics of BEC with attractive two-body interactions and in particular to study the influence of the initial conditions (using a chirp constant, $`b`$, in the wave-function) on the stability of BEC. This can be important for the condensation in optical traps and in strongly inhomogeneous traps. Also recently, in Ref. , the stability analysis was performed, using both VA and exact calculations, considering that the condensate can unstabilize by tunneling, due to quantum fluctuations.

When the two-body interaction term is switched off, for example by tunning an external magnetic field , the three-body interaction term can play an important role. Thus, it is interesting to study the collapse conditions in the model with trap potential and an attractive three-body interaction (quintic) term. Recently, in Ref. , it was investigated the collapse in one dimensional ($`1D`$) model with attractive three-body interaction. The trap was narrow and induced the radiative loss of atoms, reducing the number below the critical value, which made an arrest of the collapse possible. The two-dimensional ($`2D`$) case, with cubic nonlinearity (two-body interactions), may be useful in the modulation theory around the Townes soliton . Then it is possible to estimate the role played by the continuum.

In the present work, we first consider a review of the stability criteria, without application of initial chirp. Next, we analyze the role of an initial chirp in the wave-function of stable condensate with attractive two-body nonlinear term, which can be useful for actual experimental analysis in BEC. A generalization of this study of stability is considered with the inclusion of three-body interaction term in the NLSE. Such study can be relevant in the perspective of atomic systems with enhanced three-body effects, that can happen as the two-body scattering length is altered .

When a positive three-body interaction term is included in a trapped three-dimensional ($`3D`$) NLSE with negative two-body interaction, it is already known the possibility of extending the region of stability of the equation to larger number of atoms, with the occurrence of two stable phases in the condensate . In Ref. , it was calculated the frequency of the collective excitations in denser (liquid) and dilute (gas) phases of the condensate. Here, we also investigate a bifurcation phenomena related with the time variation of the width parameter of the theory, by switching the oscillations from one phase to the other. The present study is done with the aid of a time-dependent variational approach for the collapse of the atomic cloud in $`3D`$ model with harmonic trap potential. We compare the predictions of the VA with exact numerical simulations of partial differential equation (PDE). As detailed in our conclusions, the Gaussian VA gives a good estimate for the observables in the region of stability, and it starts to deviate from the exact results when the system is close to the collapse conditions and particularly for the unstable solutions (maxima) of the total energy. For this regions, where the Gaussian VA fails, one should improve the ansatz or take the results as a qualitative picture to guide the exact numerical calculations.

In cases without traps, the NLSE was previously studied by Vakhitov and Kolokolov (VK) , where a criterion for stability was settled. This and other criteria for stability in nonlinear systems have been studied and extensively used by many authors (see  and references therein). The VK criterion was recently detailed in Ref. . When a trap potential was added to a NLSE with a negative cubic term, an alternative stability criterion was derived in Ref. , after considering similarities between BEC atomic systems and compact object like neutron stars. Recently, the VK criterion for trapped BEC with cubic term was formally demonstrated in Ref. .

In recent numerical studies , one can also observe that the VK criterion is not generally valid; it cannot be extended, for example, to the case of a trapped system with both attractive cubic and repulsive quintic terms. The solutions can be stable irrespectively to the sign of the derivative of the eigenvalues of the NLSE with respect to the number of atoms, $`d\mu /dN`$. When $`g_3`$ is positive, and for large enough number os atoms, the NLSE is stable and $`d\mu /dN`$ can be positive (see, for example, Ref. ), implying that the criterion cannot be extended to trapped NLSE with two and three-body terms. Here, the nonvalidity of the VK criterion is clearly verified in case we have a positive harmonic trap, with two and three-body nonlinear terms with opposite signs.

We extend the study of the stability of the NLSE with cubic (two-body) and quintic (three-body) terms, through the analysis of chirp response, and through the frequency of collective excitations. By applying the perturbation techniques to the theory of nonlinear oscillations , we also derive the frequency of weakly nonlinear collective excitations. This study extend numerical calculations performed in Ref.  for the frequencies of collective excitations.

One should also note that such study is of actual interest, as recently in Ref.  is was reported observations of nonlinear oscillations in BEC of a gas with rubidium atoms. The present study can be of interest not only to atomic BEC but also to other branches of physics and mathematics, where non-linear effective interactions are added to a trapped potential in the time dependent Schrödinger equation, such as in optics and soliton physics.

In Section II, it is given a description of the model, using a variational approach for a system with $`D`$ dimensions. In Section III, we derive an expression for the chirp parameter, that is considered in our analysis of stability. The maximum initial chirp to keep the system oscillating in the same phase, around a minimum of the energy, without collapsing or without phase-transition, is obtained in this section by considering both variational and exact numerical calculations. In Section IV, we follow the study of Section III, and also obtain the frequencies for the linear and nonlinear oscillations. In this section we consider the presence of both, cubic and quintic terms, in several different configurations. The VK criterion of stability is discussed in section V. Finally, in Section VI, we present our concluding remarks.

## II Description of the model

In this section we will consider the dynamics described by the NLSE with harmonic potential and with cubic and quintic terms. Using a variational approach for a system with $`D`$ dimensions, we derive from the Lagrangian the expression for the anharmonic potential considering the mechanical analogy.

The time-dependent Gross-Pitaevskii-Ginzburg (GPG) equation that will be considered in the present approach, is given by

$$2\mathrm{i}\psi _t=\mathrm{\Delta }\psi +V\psi +\lambda _2|\psi |^2\psi +\lambda _3|\psi |^4\psi .$$
(1)
In this equation and in the following, the explicit space and time dependences of the variables and parameters are implicit, unless it is necessary or convenient for clarity. The time derivatives will be denoted by indices $`t`$. In Eq. (1) we are assuming dimensionless variables: the unit of energy is $`\mathrm{}\omega /2`$; the unit of length is $`\sqrt{\mathrm{}/(m\omega )}`$; and the unit of time is $`1/\omega `$. $`VV(\stackrel{}{r})`$ is a static trap potential, that we assume the harmonic oscillator with spherical symmetry and given by $`V=r^2`$, in the present units. $`\lambda _2`$ and $`\lambda _3`$ are the parameters of the two- and three-body interactions, which in general can be complex quantities. The imaginary parts of $`\lambda _2`$ and $`\lambda _3`$ describe, respectively, the effects of inelastic two and three-body collisions on the dynamics of BEC. In the present paper, we are not considering dissipative terms, and such cubic and quintic parameters are real. In order to compare with the formalism given in Ref. , $`\lambda _2`$ is proportional to the two-body scattering length $`a_{sc}`$ and given by $`\lambda _28\pi a_{sc}`$.

The chemical potential $`\mu `$ is given by the eigenvalue solutions of Eq. (1), with $`\psi (r,t)=e^{\mathrm{i}(\mu t/2)}\phi (r)`$:

$$\mu \phi =\mathrm{\Delta }\phi +r^2\phi +\lambda _2|\phi |^2\phi +\lambda _3|\phi |^4\phi .$$
(2)
The Lagrangian density corresponding to Eq. (1) is given by

$``$ $`=`$ $`2\mathrm{Im}(\psi _t^{}\psi )|\psi |^2r^2|\psi |^2{\displaystyle \frac{\lambda _2}{2}}|\psi |^4{\displaystyle \frac{\lambda _3}{3}}|\psi |^6.`$ (3)
### A Gaussian variational approach

To analyze the dynamics of the BEC under two and three-body interactions, it is convenient to follow the variational approach developed in Refs. . This approach was successfully employed, recently, in , for BEC with two-body interaction. We choose the simple Gaussian ansatz

$$\psi (r,t)=A(t)\mathrm{exp}\left(\frac{r^2}{2a^2(t)}+\mathrm{i}\frac{b(t)r^2}{2}+\mathrm{i}\varphi (t)\right),$$
(4)
where $`A(t)`$ is the amplitude, $`a(t)`$ is the width, $`\varphi (t)`$ is the linear phase of the condensate, and $`b(t)`$ is the “chirp” parameter previously discussed.

Without dissipative terms, the normalization of the wave-function is conserved and given by the number of particles $`N`$. The mean-square-radius and the normalization, for a system with $`D`$ dimensions, are

$$r^2=\frac{D}{2}a^2\mathrm{and}N=A^2(\sqrt{\pi }a)^D=\mathrm{constant}.$$
(5)
For a system with $`D`$ dimensions and radial symmetry, the averaged Lagrangian expression is given by

$`L=(\sqrt{\pi }a)^DA^2[2\varphi _t+{\displaystyle \frac{D}{2}}a^2(b_t+{\displaystyle \frac{1}{a^4}}+b^2+1)+`$ (6)
$`+{\displaystyle \frac{\lambda _2A^2}{2\sqrt{2}^D}}+{\displaystyle \frac{\lambda _3A^4}{3\sqrt{3}^D}}],`$ (7)
where

$$\delta L𝑑t=0,L=(r,t)𝑑\stackrel{}{r}.$$
(8)
Analyzing the corresponding Euler-Lagrange equations,

$$\frac{L}{\eta _i}\frac{d}{dt}\frac{L}{\dot{\eta }_i}=0,$$
(9)
where $`\eta _i`$ refer to the variational parameters ($`A(t)`$, $`a(t)`$, $`b(t)`$ and $`\varphi (t)`$), we obtain

$$a_t=ab;$$
(10)
and

$$b_t=\frac{1}{a^4}b^21+\frac{\lambda _2Na^{(D+2)}}{2(2\pi )^{(D/2)}}+\frac{2\lambda _3N^2a^{(2D+2)}}{3\left(\pi \sqrt{3}\right)^D}.$$
(11)
Equation (10) expresses the chirp parameter $`b`$ via initial focusing (defocusing), $`(a_t)`$, of the wave-function. Combining the Eqs. (10) and (11), we have

$$a_{tt}=\frac{1}{a^3}a+\frac{P}{a^{(D+1)}}+\frac{Q}{a^{(2D+1)}},$$
(12)
$`P`$ $``$ $`{\displaystyle \frac{\lambda _2N}{2\sqrt{(2\pi )^D}}},`$ (13)
$`Q`$ $``$ $`{\displaystyle \frac{2\lambda _3N^2}{3\left(\pi \sqrt{3}\right)^D}}={\displaystyle \frac{(2)^{(D+3)}}{(3)^{(D/2+1)}}}{\displaystyle \frac{\lambda _3}{\lambda _2^2}}P^2.`$ (14)
Comparing the result of $`a_{tt}`$ with $`D=2`$, from Eq. (12), with Ref. , we observe that the resulting equation obtained by the variational approach coincides with the one obtained by the moments method . The equation is also close to the one obtained by the modulation method for the $`2D`$ Townes soliton . In case of $`D=3`$, expressing the above definitions by the parameters given in Ref. , where $`\lambda _2=8\pi a_{sc}`$ in our dimensionless units ($`a_{sc}`$ is the two-body scattering length), and $`N=n/\left(2\sqrt{2}|a_{sc}|\right)`$:

$$|P|=\frac{n}{2\sqrt{\pi }}\mathrm{and}Q=\frac{8n^2}{9\pi \sqrt{3}}g_3.$$
(15)
The anharmonic potential is derived using the mechanical analogy. By using explicitly the dimensions in (5),

$$R\sqrt{r^2}=\sqrt{\frac{D}{2}}a\left(\frac{\mathrm{}}{m\omega }\right)^{1/2},$$
(16)
$$m\frac{d^2}{dt^2}R=\sqrt{\frac{D}{2}}a_{tt}\left[m\omega ^2\left(\frac{\mathrm{}}{m\omega }\right)^{1/2}\right],$$
(17)
a result that can be identified with

$$\frac{U_P(R)}{R}=\left(\frac{\mathrm{}\omega }{2}\right)\sqrt{\frac{2}{D}}\left(\frac{m\omega }{\mathrm{}}\right)^{1/2}\frac{U}{a},$$
(18)
where

$$U_P(R)U(a)\frac{\mathrm{}\omega }{2}.$$
(19)
The total energy for the BEC system with $`N`$ particles is given by

$$E_T=NU_P(R)=NU(a)\frac{\mathrm{}\omega }{2}E(a)\frac{\mathrm{}\omega }{2}.$$
(20)
So, for the dimensionless quantities we obtain

$$\frac{U}{a}=(D)a_{tt}.$$
(21)
With (12), we have the anharmonic potential, in dimensionless units:

$$U(a)=\frac{D}{2}\left(a^2+\frac{1}{a^2}\right)+\frac{P}{a^D}+\frac{Q}{2a^{(2D)}}.$$
(22)
Correspondingly, from Eq. (2), the chemical potential is given by

$$\mu (a)=\frac{D}{2}\left(a^2+\frac{1}{a^2}\right)+\frac{2P}{a^D}+\frac{3Q}{2a^{(2D)}}.$$
(23)
The asymptotic limit $`a\mathrm{}`$, for $`U(a)`$, is the same of the oscillator \[$`U(\mathrm{})\mathrm{}`$\], and is not determined by the parameters $`P`$ and $`Q`$. So, the only limiting condition that is strongly affected by these parameters is the limit $`a=0`$. When $`Q0`$, the three-body term is dominant in this limit: $`U(0)Q/(2a^{2D}).`$ Now, two cases have to be considered: $`Q<0`$ give us $`U(0)\mathrm{}`$; and $`Q>0`$ will give us $`U(0)+\mathrm{}`$. As the other end is fixed by the oscillator condition \[$`U(\mathrm{})+\mathrm{}`$\], in the first case ($`Q<0`$) we cannot obtain more than one maximum and one minimum for finite $`a`$, that will depend on the sign and relative value of the parameter $`P`$. The case $`Q>0`$ can be very interesting, if $`P<0`$: in both ends $`U(a)`$ goes to $`+\mathrm{}`$, such that in between we can have two minima and one maximum for the total energy. The two minima will represent two possible phases, and a rich dynamics can be described by the NLSE. From this analysis, we should note that phase transition is possible only in case $`Q>0`$ with $`P<0`$. When $`Q0`$, there is only a single phase.

### B Critical parameters for stability

The relevant extremes of $`U(a)`$ are given by the real and positive roots of

$$a_s^2\frac{1}{a_s^2}\frac{P}{a_s^D}\frac{Q}{a_s^{2D}}=0.$$
(24)
In such extreme positions, we obtain

$$U(a_s)=Da_s^2+\frac{P}{2a_s^D}\left(2D\right)+\frac{Q}{2a_s^{2D}}\left(1D\right),$$
(25)
$$\mu (a_s)=Da_s^2+\frac{P}{2a_s^D}\left(4D\right)+\frac{Q}{2a_s^{2D}}\left(3D\right).$$
(26)
In Table I we resume the particular cases, from the above three equations, where exact analytical solutions can be found for the extremes of $`U`$ and the corresponding critical parameters for stability. We note that, once a solution for a system with $`D`$ dimensions and $`P=0`$ is known, the same solution can be used for the roots of a system with $`2D`$ dimensions and $`Q=0`$, by just exchanging $`P`$ and $`Q`$. When $`D=2`$, the result $`P_c=1`$ for the particular case when the quintic term is zero ($`Q=0`$) is well known, and corresponds to $`\nu =2`$ in Ref. . The solutions will be unstable for $`P=|P|`$ with $`|P|>1`$.

Actually, one of the interesting cases occurs when the cubic term is zero ($`P=0`$), for $`D=1`$ and $`D=3`$, in view of the possibility of altering the two-body scattering length by an external magnetic field . If this condition is realized ($`P=0`$), the effect of the three-body non-linear term in the mean-field approximation of the condensate will be enhanced.

In case such three-body (quintic) term is negative, as shown in Table I, the present variational ansatz gives analytical estimates for the critical parameters, and for physical quantities as the mean-square-radius, energy and chemical potentials. In Section IV we study in more detail this case.

Other particular cases occur when phase transitions are possible. As explained before, these situations can only happen if both parameters $`P`$ and $`Q`$ are non-zero and $`Q>0`$ with $`P<0`$. The absolute collapse is not possible in the condensate; but, one can obtain transitions between two-phases. This implies that, for a system with two-body attractive interaction and three-body repulsive one, as we increase the number of atoms in the condensate ($`|P|`$ increases), we can reach a critical limit where only one phase (stable) remains. In other words, for a particular value of $`P=P_c`$, the position $`a`$ of one of the minima of $`U`$ coincides with the maximum, and we have an inflection point of $`U`$.

As an example, in Table I we show one particular case where we have phase transition, and also analytical solution for one of the roots of Eq. (24), when $`D=3`$. We can observe that this condition is realized if $`0>PP_c=(1/8)^{(1/4)}`$. When $`Q=P^8`$, the critical limit occurs for $`P_c=(1/8)^{(1/4)}`$ and we have three-roots for $`a_s`$, given by $`a_{s1}=a_{s2}=|P_c|=0.5946`$ and $`a_{s3}=`$0.4696.

A typical plot of the potential $`U(a)`$ as a function of the width $`a`$ and the corresponding phase portrait defined in the plane ($`a_t,a`$) for negative $`\lambda _2`$ and positive $`\lambda _3`$ is given in Fig. 1 for the case $`D=3`$ of a stable BEC. The parameters for the attractive cubic and repulsive quintic terms used in the figure are such that $`n=`$1.948 and $`g_3=`$0.016, corresponding to $`P=`$0.5495 and $`Q=`$0.0099. The number of atoms corresponds to the situation in which we have two minima with the same value for $`U(a)`$. From the variational expressions given in this section and our mechanical analogy, one should observe that, in the pictorial example of Fig. 1, if we have a solution located at the right minimum (for example) it cannot migrate to the left minimum, unless the chirp parameter corresponds to an energy greater than the difference between the right minimum and the maximum that is in between.

## III Analysis of the oscillations of the BEC in three dimensions

The different dynamical regimes in the condensate oscillations can be described by the cross sections of the curve of $`U(a)`$ with the levels of the effective total energy. Following our mechanical analogy, we have

$$H(P_R,R)=\frac{P_R^2}{2m}+U_P(R),$$
(27)
where $`P_R`$ is the momentum conjugate to $`R`$, given in (16). In our dimensionless variables, the corresponding initial Hamiltonian is given by

$$H_0=\frac{3}{2}(a_t)^2|_0+U(a_0)=\left(\frac{3}{2}\right)a_0^2b_0^2+U(a_0),$$
(28)
where Eq. (10) was used. Thus, by varying the initial conditions $`a_0,b_0`$ we alter $`H_0`$.

This case was recently analyzed using the moments methods in Refs. , under some specific assumptions. The assumption $`a_t|_0=0`$ was considered in Refs. ; and $`a_t|_00`$ in Ref. . In Ref. , a generalization of the Weinstein criterion for the collapse  was obtained, in case of three dimensions, and when there is only cubic term in the NLSE.

It will be interesting to derive criteria for stability more generally applied, using a variational approach and by comparison with results obtained with numerical simulations. For this purpose, we should first obtain expressions for the mean-square-root that extremizes the total energy. In the variational approach, for $`D`$ dimensions, the equations for the total energy (22) and chemical potential are given by (22) and (23). So, in case of $`D=3`$, the equations corresponding to (24), (25) and (26) are

$$a_s^8a_s^4Pa_s^3Q=0,$$
(29)
$$U(a_s)=2a_s^2+\frac{1}{a_s^2}+\frac{P}{2a_s^3}=\frac{5a_s^2}{2}+\frac{1}{2a_s^2}\frac{Q}{2a_s^6},$$
(30)
$$\mu (a_s)=3a_s^2+\frac{P}{2a_s^3}=\frac{7a_s^2}{2}\frac{1}{2a_s^2}\frac{Q}{2a_s^6}.$$
(31)
Here we should note that an exact expression for the Hamiltonian, corresponding to the VA given by Eqs. (22), can be derived for the general case. This will be useful for the exact numerical calculations that we are going to perform in 3D. We consider the following scaling in the wave-function $`\psi `$ of the system :

$$\psi =\alpha ^{3/2}\chi (\stackrel{}{\xi }),\mathrm{with}\stackrel{}{\xi }\frac{\stackrel{}{r}}{\alpha }.$$
(32)
Then, for the total Hamiltonian, we obtain

$$H(\alpha )=\frac{X_s}{\alpha ^2}+\alpha ^2\xi ^2+\frac{\lambda _2}{2}\frac{Y_s}{\alpha ^3}+\frac{\lambda _3}{3}\frac{Z_s}{\alpha ^6},$$
(33)
where

$$X_s\chi ^{}\mathrm{\Delta }\chi d^3\xi ,\xi ^2=\xi ^2|\chi |^2d^3\xi ,$$
(34)
$$Y_s|\chi |^4d^3\xi ,Z_s|\chi |^6d^3\xi .$$
(35)
By taking $`\delta H(\alpha )/\delta \alpha |_{\alpha =1}=0`$, we obtain the characteristic equation

$$X_sr^2+\frac{3}{4}\lambda _2Y_s+\lambda _3Z_s=0.$$
(36)
With the above equation, for the energy we have

$$H=\frac{2X_s}{3}+\frac{4}{3}r^2\frac{\lambda _2}{4}Y_s.$$
(37)
Correspondingly, for the chemical potential, we obtain

$$\mu =2r^2+\frac{\lambda _2}{4}Y_s.$$
(38)
Interesting to observe that when the two-body interaction is zero, the exact result for the chemical potential is twice the mean square radius, as shown above \[in agreement with the variational result of Eq.(31), where $`a_s^2=2/3r^2`$\]. As specific cases will be considered in the next subsections, we will discuss this point in more detail. For now, we can use the specific example of $`Q=0`$, represented in Fig. 2, in order to illustrate the validity of a general criterion (also valid for $`Q0`$) that we consider. In this case, we have just one minimum and a maximum of $`U(a)`$.

A given chirp $`b(t)`$ in the wave-function corresponds to the square root of the kinetic energy, given by the Hamiltonian in Eq. (28). With a small initial chirp, we should have oscillations around the minimum of $`U(a)`$. The square of the chirp must be larger than the difference between the minimum and the maximum of the effective potential energy $`U(a)`$ for the system to become unstable. So, from Eq. (28) and given that $`H_0`$ is conserved, we can derive a minimum criterion for stability. The given kinetic energy (related to the chirp) must be smaller than the corresponding variation of the potential energy, in order to maintain stable oscillations around the minimum (where $`a=a_{s1}`$). The initial chirp $`b_0`$ must satisfies the condition

$`b_0^2a_{s1}^2`$ $``$ $`{\displaystyle \frac{2}{3}}\mathrm{\Delta }U(a),`$ (39)
$`b_0^2`$ $``$ $`{\displaystyle \frac{2}{3}}{\displaystyle \frac{U(a_{s2})U(a_{s1})}{a_{s1}^2}},`$ (40)
where $`a_{s2}`$ is the value of $`a`$ corresponding to the maximum of $`U(a)`$, and $`a_{s1}`$ corresponds to the minimum of $`U(a)`$. The maximum initial chirp is given by

$`b_{0,m}`$ $`=`$ $`\sqrt{{\displaystyle \frac{2}{3}}}{\displaystyle \frac{\sqrt{U(a_{s2})U(a_{s1})}}{a_{s1}}}.`$ (41)
From the above equation, we should also note that the critical values for the number of atoms, related to the parameters $`P`$ and $`Q`$, can be obtained by the condition that $`U(a_{s2})U(a_{s1})`$. We should also observe that the chirp parameter is introduced in the exact wave function by an exponential factor, as in the VA \[see (4)\].

The frequency of the linear collective oscillations, $`\omega _L`$, is also another relevant quantity to analyze the stability of the condensate. When the system becomes unstable, it reaches the value zero. From the radial variation of the force around a minimum $`a_{s1}`$, obtained from Eqs. (29) and (30), $`\omega _L`$, in units of the trap frequency $`\omega `$, is given by

$`\omega _L=\sqrt{{\displaystyle \frac{1}{3}}{\displaystyle \frac{d^2U}{da^2}}|_{a=a_{s1}}}`$ $`=`$ $`\sqrt{5{\displaystyle \frac{1}{a_{s1}^4}}+{\displaystyle \frac{3Q}{a_{s1}^8}}}`$ (42)
$`=`$ $`\sqrt{8{\displaystyle \frac{4}{a_{s1}^4}}{\displaystyle \frac{3P}{a_{s1}^5}}},`$ (43)
One of the requirements for instability of the BEC is that $`\omega _L`$ is zero, when the minimum and maximum of $`U`$ disappear. So, it is convenient to analyze the critical points when the system becomes unstable; at such points, the system collapses (in case we have only one minimum) or one of the phases (characterized by a minimum) disappears. From Eq. (30),

$`{\displaystyle \frac{1}{3}}{\displaystyle \frac{d^2U}{da^2}}|_{a=a_s}`$ $`=`$ $`{\displaystyle \frac{5a_s^8a_s^4+3Q}{a_s^8}}`$ (44)
$`=`$ $`{\displaystyle \frac{8a_s^54a_s3P}{a_s^5}}.`$ (45)
We can only have two minima (given by $`a_{s1}`$ and $`a_{s3}`$) and one maximum (at $`a_{s2}`$) for negative $`P`$ and positive $`Q`$. In this case, two critical limits are possible, corresponding to the situations in which the position of one of the minima, $`a_{s1}`$ or $`a_{s3}`$, is equal to the position of the maximum (one of the phases disappears). Note that in the critical limit we have $`a_{min}=a_{max}=a_c`$. In the next two subsections, we will consider the particular cases in which we have only one nonlinear term present in the NLSE: the cubic term or the quintic term. The discussion of the more general cases, where both terms are present in the NLSE, we leave for the section IV.

### A Case of $`Q=0`$, with attractive two-body term

The roots for the maximum ($`a_{s2}`$) and the minimum ($`a_{s1}`$) of $`U`$ can be obtained from (29):

$$a_s^5a_s+|P|=0.$$
(46)
Combining this equation with Eq. (45), we have

$`{\displaystyle \frac{1}{3}}{\displaystyle \frac{d^2U}{da^2}}|_{a=a_s}`$ $`=`$ $`5{\displaystyle \frac{1}{a_s^4}}.`$ (47)
In the critical limit ,

$`a_c`$ $`=`$ $`\left({\displaystyle \frac{1}{5}}\right)^{\frac{1}{4}},P_c={\displaystyle \frac{4}{5}}\left({\displaystyle \frac{1}{5}}\right)^{\frac{1}{4}}=0.53499`$ (48)
$`\mathrm{and}n_c`$ $`=`$ $`2\sqrt{\pi }|P_c|=1.8965,`$ (49)
where Eqs. (46) and (15) were used. With the above root of Eq. (46), the other roots can also be easily found numerically for any value of $`P`$. Given $`P`$, once we obtain the positions of the minimum ($`a_{s1}`$) and maximum ($`a_{s2}`$) of $`U`$, we can calculate the value of the maximum initial chirp for the collapse, using Eq. (41):

$$b_{0,m}=\sqrt{\frac{a_{s1}^2a_{s2}^2}{3a_{s1}^4a_{s2}^2}(15a_{s1}^2a_{s2}^2)}.$$
(50)
Our results, when we have only the cubic term in the nonlinear effective potential ($`Q=0`$), is shown in Figs. 2 and 3. In Fig. 2, for $`P=0.1693`$ (corresponding to $`n=0.6`$), we represent in the lower frame the phase space given by $`a_t`$ versus $`a`$. In the upper frame, we have the total energy, in units of $`N(3/2)\mathrm{}\omega `$, as a function of $`a`$.

In Fig. 3, we show in the lower frame our results for the maximum initial chirp, $`b_{0,m}`$, as a function of $`n`$, for stable solutions. $`n`$ is the reduced number of atoms, given by Eq. (15). In this figure, we also represent the total energy (upper frame) and the average of the square radius, $`r^2`$, as a function of $`n`$.

The static variational expression for the maximum initial chirp, was derived from Eq. (10) and energy conservation requirements, given by Eq. (40). The circles correspond to results obtained by the time-dependent VA, by increasing the initial conditions for the chirp till the system collapses. We also obtain the results for the maximum initial chirp using the same Eq. (40), but considering exact PDE results for the observables presented in it. The results, for the static case, are represented by a solid line; and, for the time-dependent case, by squares. We have observed small differences in the absolute values of $`b_{0,m}`$ when considering the trial $`b_0`$ as positive or negative. This fact can be qualitatively understood by a small dephasing of the wave that occurs near the wall of the oscillator, when it is bouncing back.

For the mean-square-radius, we observe that the VA results approaches the exact PDE results also in the unstable branch in the limit $`n0`$. These results are consistent with the ones obtained in Ref. , where it is shown that, for the trapped NLSE with cubic term, the unstable branch of the chemical potential diverges in the limit $`n0`$, with $`r^2`$ collapsing to zero (this limit is the non-trapped solution of the NLSE).

### B Case of $`P=0`$ with attractive three-body term

As in the previous subsection, for a stable condensate the maximum initial chirp is given by Eq. (41), where the roots are given by Eq. (29). The real and positive roots correspond to stable solutions, such that only one is meaningful when $`Q>0`$. When $`Q<0`$, we obtain two real and positive roots of $`a`$, in case $`0>Q1/4`$. These roots are

$$a_{s1/s2}a_{+/}=\left[\frac{1}{2}\pm \sqrt{\frac{1}{4}+Q}\right]^{\frac{1}{4}},$$
(51)
where the minimum of $`U`$ is given by $`a_+`$ ($`a_{s1}`$) and the maximum by $`a_{}`$ ($`a_{s2}`$). The critical limit happens for $`a_c^4=1/2`$, or $`Q_c=|Q_c|=1/4`$. From Eqs. (51) and (45),

$$\frac{1}{3}\frac{d^2U}{da^2}|_{a=a_s}=\pm \frac{2\sqrt{14|Q|}}{|Q|}\left(1\sqrt{14|Q|}\right),$$
(52)
and the frequency of the linear collective oscillations, $`\omega _L`$, near the minimum is given by

$$\omega _L=\sqrt{\frac{2\sqrt{14|Q|}}{|Q|}\left(1\sqrt{14|Q|}\right)}.$$
(53)
In Fig. 4, we present the main results for this subsection with $`P=0`$. As noted in Eqs.(31) and (38) (respectively, the variational and exact expressions), the chemical potential is just twice the mean-square-radius. In the plots, the VA results are shown in dashed lines and the exact PDE in solid lines. The variational critical limit, $`Q=Q_c=1/4`$ ($`a_c^4=1/2`$), corresponds to $`n_{3,c}=1.2371,`$ where $`n_3`$ is the reduced number of atoms, that is given by the normalization of the wave-function. In this case, Eq. (14) will give us the relations between $`|Q|`$, $`|\lambda _3|`$ and $`n_3`$:

$$n_3=\frac{\sqrt{|\lambda _3|}N}{2\pi }\mathrm{and}|Q|=\frac{8n_3^2}{9\pi \sqrt{3}}.$$
(54)
Given $`\lambda _3`$, we can obtain the corresponding critical number of atoms:

$$N_c=\frac{3\pi ^{3/2}3^{1/4}}{2\sqrt{2}\sqrt{|\lambda _3|}}=\frac{7.7729}{\sqrt{|\lambda _3|}},$$
(55)
a result consistent with no limit in the number of atoms if $`\lambda _30`$.

In Fig. 4, we also observe that $`\mu /n_3<0`$ (branch of minima) and $`\mu /n_3>0`$ (branch of maxima) correspond, respectively, to stable and unstable solutions. We note that localized stable structure can only exist for positive $`\mu `$. The results in this case, that the effective interaction contains only the trap and quintic term, are in agreement with the VK criterion . The observed exact relation between the chemical potential and the mean-square-radius ($`\mu =2r^2`$), in the present case that $`\lambda _2=0`$, was proved in Eq. (38).

In the branch with the unstable solutions, we observe numerically that in the limit $`n_30`$ the wave-function approach the form $`\psi \mathrm{exp}(r^2/2)/r`$; with $`\mu 1`$ and $`r^2`$1/2. This solution is proved analytically to be the irregular solution of the oscillator \[Recall that, when $`\lambda _2=0`$ and $`\lambda _3=0`$, in Eq.(2), the regular oscillator solution gives $`\mu =3`$; as for the irregular solution in the origin we have $`\mu =1`$ in our energy unit ($`\mathrm{}\omega /2`$)\]. So, for the unstable solutions, with small $`|Q|`$, the Gaussian ansatz (started by considering the exact regular solution of the oscillator in the limit $`n_3=0`$) fails. By keeping (artificially) the Gaussian shape for the unstable solutions, the radius is forced to be zero. This explains the discrepancy between the present VA and the PDE results in the unstable region. The more appropriate ansatz for the trapped NLSE (when $`\lambda _2`$ and/or $`\lambda _3`$ is nonzero) should be build considering the two (regular and irregular) solutions of the harmonic oscillator.

Even considering the discrepancy between VA and PDE results for the unstable (maxima) solutions, there is a reasonable agreement for the results obtained for the maximum initial chirp $`b_{0,m}`$, as shown in Fig. 5. In Fig. 6, we have the variational result for the potential energy as a function of $`a`$, and the corresponding phase space, as in Figs. 1 and 2. Static predictions on PDE are also in good agreement for the maximum initial chirp. In PDE calculations we have verified numerically small differences between negative and positive $`b_{0,m}`$, as in the case analyzed and explained in the previous subsection (related to Fig. 3). Using Eq. (41) with $`P=0`$, the maximum initial chirp is given by

$`b_{0,m}`$ $`=`$ $`\sqrt{{\displaystyle \frac{1}{3|Q|}}\left(12\sqrt{|Q|}\right)\left[\sqrt{14|Q|}1+2\sqrt{|Q|}\right]}`$ (56)
$``$ $`\sqrt{{\displaystyle \frac{2}{3\sqrt{|Q|}}}}\left(1{\displaystyle \frac{3\sqrt{|Q|}}{2}}\right)\mathrm{for}|Q|<<1.`$ (57)
The above result shows that, when the three-body parameter is negative and small, the maximum initial chirp is $`|Q|^{1/4}`$, or $`(1/\sqrt{n}_3)`$, as shown in Fig. 5.

## IV Dynamics of BEC with two- and three-body interactions

In this section, we analyze the cases where we have non-zero cubic ($`\lambda _2`$) and quintic ($`\lambda _3`$) interaction terms in the NLSE.

The different possibilities, relative to the signs of the two and three-body interactions are studied in the next three sub-sections, where we exclude one of the cases in which the signs are the same. When $`Q`$ has the same sign of $`P`$, the basic physical picture is not essentially altered in comparison to the cases already discussed, where one of these quantities is zero. However, it is worthwhile to examine the changes in the collective excitations in the case that $`\lambda _2`$ and $`\lambda _3`$ are positive, as such observable can be useful to obtain informations about possible manifestation of three-body interactions.

### A The case of repulsive cubic and quintic interactions ($`\lambda _2>0`$ and $`\lambda _3>0`$)

From Eq. (43), one should observe that the frequency of the collective oscillations will decreases, as we increase $`P`$, or if the position of the minimum, $`a_{s1}`$, decreases. As all quantities present in Eq. (43) are positive, this frequency in the VA has an upper limited at the value $`2\sqrt{2}`$, such that a reasonable frequency of oscillations around the minimum will be smaller than that (see also Ref. ).

The effects of the unharmonicity of oscillations can be taken into account by the expansion of the effective potential near the bottom $`a_c`$ on the power of deviations $`y=aa_{s1}`$ . The equation for $`y(t)`$ is:

$$y_{tt}=B_1yB_2y^2B_3y^3$$
(58)
where

$`B_1=1+{\displaystyle \frac{3}{a_{s1}^4}}+{\displaystyle \frac{4P}{a_{s1}^5}}+{\displaystyle \frac{7Q}{a_{s1}^8}},`$ (59)
$`B_2={\displaystyle \frac{6}{a_{s1}^5}}{\displaystyle \frac{10P}{a_{s1}^6}}{\displaystyle \frac{28Q}{a_{s1}^9}},`$ (60)
$`B_3={\displaystyle \frac{10}{a_{s1}^6}}+{\displaystyle \frac{20P}{a_{s1}^7}}+{\displaystyle \frac{84Q}{a_{s1}^{10}}}.`$ (61)
Then, applying the perturbation techniques to the theory of nonlinear oscillations , we find the frequency of weakly nonlinear collective excitations:

$$\omega _{NL}=\sqrt{B_1}\left[1+\left(\frac{3B_3}{8B_1}\frac{5B_2^2}{12B_1^2}\right)\sigma ^2\right],$$
(62)
where $`\sigma `$ is the amplitude of the oscillations. The correction to linear frequency is proportional to the square of the amplitude of oscillations of the condensate.

The above estimate of the nonlinear oscillations in BEC, when both the two and three-body terms are positive, can be a relevant piece of information to determine a possible manifestation of three-body interaction in the condensates and its magnitude. We should observe that, recently, in Ref. , it was reported observation of nonlinear oscillations in BEC of a gas with rubidium atoms.

Let us estimate the amplitude of oscillations of the width of the atomic cloud. The points of maximum and minimum width are defined by the section points of the line $`H_0=const`$ with the potential curve $`U(a)`$. For large width, $`a>>1`$, we can obtain the estimate

$$\sigma _\pm =\sqrt{\frac{H_0}{3}\pm \sqrt{\left(\frac{H_0}{3}\right)^21}},$$
(63)
where $`\sigma _{+()}\sigma _{max(min)}`$. It is natural to assume that the asymptotic value of the width is the averaged value between $`\sigma _{max}`$ and $`\sigma _{min}`$.

### B Case of attractive cubic and repulsive quintic interactions ($`\lambda _2<0`$ and $`\lambda _3>0`$)

In this case, there is no collapse, but the condensate can have up to two distinct phases when $`|\lambda _2|`$ is smaller than a critical value. From the beginning, when $`\lambda _3=0`$, $`U`$ has one minimum at $`a=a_{s1}`$ and a maximum at $`a=a_{s2}`$. By considering a fixed positive (and small) $`\lambda _3`$, a second minimum of $`U`$ ($`a_{s3}`$) appears, corresponding to a denser phase in the condensate. As we increase the value of $`|P|`$, we can reach a critical value where the minimum of the normal phase disappears .

When $`H_0/3<\mathrm{\Delta }U=|U(a_{s2})U(a_s)|`$, with $`a_s=a_{s1}`$ or $`a_{s3}`$, we have small amplitude oscillations near the fixed point $`a_s`$. With $`H_0/3>\mathrm{\Delta }U`$, the character of the oscillations changes and we have large amplitude oscillations , due to the motion of the effective particle between the wall given by the repulsive three-body interaction and the other curve given by the quadratic potential. With $`a_s`$ close to $`a_{s3}`$ and $`H_0/3<\mathrm{\Delta }U_2=|U(a_{s2})U(a_{s3})|`$, the condensate oscillates with the amplitude restricted from below by $`a_{s2}`$.

The position of the minimum $`a_{s1}`$ of $`U(a)`$ is defined by the equilibrium between two and three-body interactions; the positions of the other extrema of $`U`$, the fixed points $`a_{s2}`$ and $`a_{s3}`$ are mainly defined by the contributions of the two- and three-body terms. These considerations allow us to obtain analytical expressions for such fixed points, that are given by the real and positive roots of the equation given by the first derivative of $`U(a)`$.

In order to analyze the frequencies of the linear oscillations \[given by Eq. (43)\], we return to Eq (45), where $`P=|P|`$ and $`Q=|Q|`$. The physically relevant situations occur when $`|P|`$ is a fraction of one and $`Q<<1`$. And of particular interest are the cases where two phases (corresponding to two minima) are possible, at $`a=a_{s1}`$ and $`a=a_{s3}`$.

Let us consider, for example, the particular case examined in Table I, that we have a one-point solution given by $`a_{s1}=|P|`$ and $`Q=a_{s1}^8=|P|^8`$. This example is chosen not only due to convenience but because it is not far from a more realistic situation. Usually we scale the normalization of the wave-function with $`n`$ that is directly proportional to $`|P|`$. So, the three-body parameter $`\lambda _3`$ is related to the particular point we are considering. The critical limit, in this case, is given by $`|P_{c1}|^4=(1/8)`$ ($`|P_{c1}|=`$0.5946) and $`Q_{c1}=1/64`$, as shown in Table I, corresponding, respectively, to $`n_{c1}=2.1078`$ and $`g_3=0.021`$ \[see Eq. (15)\]. At such critical limit the frequency of the linear oscillations around the minimum $`a_{s1}`$ goes to zero \[In a practical situation, before the system reaches this critical limit in the normal phase, it tunnels to the denser phase, when both energies are equal\]. With $`|P|`$ larger than $`|P_c|`$, there is only one phase (minimum) corresponding to a denser phase. As we reduce the value of $`|P|`$, we reach a second critical point at $`a_{s2}=a_{s3}=|P_{c3}|`$, where $`a_{s3}`$ is the position of the minimum of the denser phase. A given value of $`|P|`$ should be in between the limits $`|P_{c3}|<|P|<|P_{c1}|`$ for the existence of two phases.

So, for the normal phase, we note that the corresponding frequencies of the oscillations are approximately given by the harmonic oscillator, with $`\omega _{L1}2`$ (twice the trap frequency). This can be seen in the previous cases where one of the nonlinear terms is zero and the other is small.

Consider the case exemplified in Fig. 1, where $`P0.55`$ ($`n`$1.95) and $`Q`$ 0.01 ($`g_3=0.016`$), and the corresponding roots for the minima are $`a_{s3}`$0.386 and $`a_{s1}`$0.726. Using Eq. (43), we can observe that the frequency of the oscillations in the denser phase, $`\omega _{L2}=4.493`$, is much larger than the frequency of the normal phase, $`\omega _{L1}=1.339`$. If we take more realistic parameters as, for example, $`g_3=0.001`$ and $`n=1.5`$, using Eqs. (15) and (29), we obtain $`a_{s1}=0.840`$ (minimum of the ‘gas’ phase), $`a_{s2}=0.434`$ (maximum) and $`a_{s3}=0.105`$ (minimum of the ‘liquid’ phase). The corresponding frequencies of the linear oscillations of the two phases are $`\omega _{L1}=1.73`$ and $`\omega _{L2}=258`$. So, in general we observe that it should be expected

$$\omega _{L2}\omega _{L1}.$$
(64)
This unequality for the frequencies was also confirmed in numerical PDE calculations done in Ref. .

The amplitude of the nonlinear oscillations in the dense phase can be found from the observation that $`a_{s3}<<1`$. Then the terms with $`a^6`$ and $`a^3`$ are dominant in the effective potential energy, and an approximate solution can be found:

$$\sigma _{max,min}=\left(\frac{|P|}{2|H_0|}\pm \sqrt{\frac{P^2}{4H_0^2}\frac{Q}{2|H_0|}}\right)^{1/3}.$$
(65)
An interesting phenomenon occurs when the initial $`a_t|_0=a_0b_0`$ is large. Then, the character of cloud oscillations will change, from oscillations near $`a_{s1}`$ and $`a_{s3}`$ to large amplitude oscillations defined mainly by the quadratic potential (see Fig. 1). Let us estimate the criterion for the bifurcation phenomenon in the oscillations. Taking into account the expression for the energy $`H_0`$, the bifurcation point is given by

$$b_{0,m}\sqrt{\frac{2\mathrm{\Delta }U}{3}}\frac{1}{a_{s1}}.$$
(66)
The result is shown in Fig. 7, using our Gaussian VA, which gives us a qualitative picture of the exact results, that were already presented in Ref. .

In Fig. 8 we show full numerical PDE calculation for the system beginning in the gas phase, which corresponds to the right minimum of the Fig. 1. Our results show that for a small chirp $`(b_00.10)`$, the mean square radius oscillates with small amplitude, according to its collective frequency. Applying a stronger chirp, near $`b_00.11`$, we observe transitions back and forth between two phases, characterized by different amplitudes of the oscillations. As $`b_0`$ is large enough we observe that the oscillation patterns remains with almost fixed and large amplitude. A similar picture happens if we begin in the denser phase, with different initial conditions given by the chirp. These results were obtained by using the exact PDE calculation, and they are in qualitative agreement with the VA approach. The quantitative agreement between the VA and the exact PDE deviates considerably in the phase transition region. This should be expected, as we are far from the harmonic oscillator behavior, in a region where an improvement in the ansatz is necessary.

### C Case of repulsive cubic and attractive quintic interactions ($`\lambda _2>0`$ and $`\lambda _3<0`$)

As both two and three-body nonlinear terms are non-zero and with opposite signs (with attractive three-body and repulsive two-body interactions), one could expect a behavior similar to the case that was analyzed in the previous subsection, and also represented in Fig. 1 (where the two-body interaction is attractive and the three-body is repulsive). But contrary to such expectation, this case shows a different behavior for small radius, and no phase transition is possible. We observe that in this case the system can collapse, as the behavior for small radius is dominated by the three-body term, that is negative. This is represented in Fig. 9, in a variational plot of $`U(a)`$ versus $`a`$, together with a corresponding plot for the phase space.

In Fig. 10, we present the breathing mode collective excitations for a few values of $`g_3`$ calculated in the variational approximation. The collective excitations show that even for small negative $`g_3`$ a limited number of atoms is allowed in the condensate. Now, only a region of stationary condensate - the denser phase - exists. Two points exist that extremize the energy: one, $`a_{s1}`$, is stable; and the other, $`a_{s2}`$, is unstable. Here, the situation is similar to the cases that we have a single attractive nonlinear term (cubic or quintic); and the conditions for collapse, in terms of the initial chirp $`b_0`$, is such that it must be larger than $`b_{0,m}`$ given in Eq. (41).

One should note that this particular case ($`\lambda _2>0`$ and $`\lambda _3<0`$) can be relevant in the analysis of experiments with BEC performed with atomic systems that have repulsive two-body interaction. No collapse is expected if a real three-body effect is not manifested, or in case the possible three-body effect is also repulsive. However, an attractive three-body effect will change this scenario, as the system must collapse for certain critical maximum number of atoms. In this perspective, the present analysis shows that experiments with BEC may be useful in detecting negative three-body forces. The corresponding maximum critical number of atoms for stability of the condensate can also be obtained from the present approach, by taking the limit $`b_{0,m}0`$ (or $`a_{s2}a_{s1}`$). For instance, when $`g_3=0.5`$, the critical number $`n`$ can be obtained from the rhs (positive $`\lambda _2`$) of Fig. 11.

We present in Fig. 11, for a fixed attractive three-body interaction, with $`g_3=0.5`$, the total energy $`E`$, the mean-square-radius $`r^2`$, and the chemical potential, as a function of the reduced number of atoms $`n`$ multiplied by the sign of the two-body interaction. We consider both negative ($`\lambda _2<0`$) and positive ($`\lambda _2>0`$) two-body interactions, and the value of $`g_3=0.5`$ was chosen due to numerical convenience. The variational approach (dashed line) gives a good description of the minima for small $`|n|`$, but fails particularly for the maxima solutions, as one can see in the figure when comparing with the exact results (solid lines). In agreement with the variational results of collective excitations, the number of atoms allowed is limited to a critical number $`n_c`$ 3.5 (the corresponding exact result is $``$2.93). For a more realistic value of the three-body parameter, with $`g_3=`$0.01, the exact critical number is $`n_c`$4200.

In the upper frame of the Fig. 11, it is shown a similar pattern for the minimum and the maximum of the energy, as already described with $`\lambda _30`$ when $`\lambda _2<0`$ and with $`\lambda _3<0`$ when $`\lambda _2=0`$. Such results also indicates that an initial wave function with enough larger chirp can make the condensate unstable and that the corresponding magnitude of $`b_{0,m}`$ can be inferred from the absolute difference between the maximum and minimum of the energy. This implies that, when we have positive two-body interaction ($`\lambda _2>0`$), the presence of a negative three-body term ($`\lambda _3<0`$) can be detected in principle even if the number of atoms has not achieved the critical number. The initial maximum chirp plays a relevant role in this case.

## V VK criterion

In the context of stability analysis of NLSE bound states, the orbital stability, or stability with respect to the form, was first settled by the so-called Vakhitov-Kolokolov criterion , that is detailed in a recent review in Ref. . A solitary wave is orbitally stable if the initial orbit, chosen near the ground-state orbit, implies that the orbit of the solution at any $`t>0`$ remains close to the ground state. The criterion was first demonstrated in Ref. , by minimizing the deviations between the orbital states and the ground-state, with respect to the initial position parameters. In our notation, the VK criterion for stability can be written as We should note that $`\mu `$ corresponds to $`\lambda `$ of Ref. .

$$\frac{N}{\mu }<0.$$
(67)
However, the original VK criterion and the above cited demonstration have only considered cases without trap, with nonlinearities expressed by $`f(|\psi |^2)`$. More recently, the general conditions that guarantees orbital stability of stationary trapped condensates, described by the NLSE with cubic term, was formally demonstrated by Bergé et al. , in agreement with results analyzed in Ref.  for the case that the cubic term is negative.

The validity of the VK criterion appears to be solidly consolidated by the above described formal demonstrations and by the specific cases that have so far been considered. These results appear to support an extended range of applicability of the VK criterion to the NLSE. In the several cases that we have considered in the present paper (by exact numerical procedure and also by the Gaussian variational approach), the applicability of the VK criterion is out of the question in the mentioned cases that have been examined by other authors. We have also confirmed the applicability of the criterion in the case that we have only the trap and the quintic term in the effective interaction, as shown in Fig. 4.

The validity of the VK criterion cannot be extended to trapped NLSE with two and three-body terms, when the two body term is negative, as one can observe from the results obtained in Ref.. This is an unexpected result, considering that the criterion is applicable for non-trapped systems (with cubic and/or quintic terms) and also for the trapped cases when we have only the cubic or the quintic terms (and also in the cases that both terms have the same sign). The VK criterion is also not applicable when $`\lambda _2>0`$, as observed in Fig. 11, in the right side of the plot. The criterion fails particularly in the region near the oscillator solution (where $`\mu =`$3), where we can observe that $`\mu /n>0`$ in a stable branch (upper) and also in the unstable branch (lower). (We recognize the stable branch as the one that corresponds to minima of the energy.) Therefore, in the present work, we confirm the conclusion about the limitation of the VK criterion: The nonvalidity of the VK criterion is verified numerically when we have a positive harmonic trap, with two and three-body nonlinear terms with opposite signs.

## VI Conclusions

The stability of a trapped condensate, with two- and three-body nonlinear terms, was studied in the present work, considering several aspects, as the initial conditions in the wave-function and the validity of the Vakhitov-Kolokolov criterion. For the initial conditions we considered a chirp parameter, which is related to the initial focusing (defocusing) of the cloud. A non-zero initial chirp $`b_0`$ introduces oscillations in the condensate near the minimum of the energy, such that it can lead a previously stable system to collapse, when $`b_0`$ reaches a maximum limit $`b_{0,m}`$, that corresponds to the energy difference between a minimum and a maximum of the total potential energy. So, in the presence of a fixed value of the chirp, the number of atoms of a stable condensate with attractive interaction is smaller than the corresponding value when the chirp is zero. In case that the potential has two minima (when cubic and quintic terms are present in the NLSE), related to two phases of the condensate, a bifurcation phenomenon is predicted in the dependence of the value of the initial chirp. A maximum initial chirp can affect the system in such a way that the oscillations can switch from the gas phase to the liquid phase and vice-versa.

The present study was performed by using exact numerical solutions of partial differential equation, as well as by a corresponding variational approach. The analytical predictions based on the time-dependent VA, using Gaussian ansatz, are qualitatively confirmed by the exact time-dependent numerical simulations. We would like to point out the advantages of the Gaussian VA when comparing with the moments method, which encounters some difficulties in deriving the collapse conditions even in the pure attractive two-body case, for some classes of the initial data (the sign of the time dependent function in the equation for the mean-square-radius, $`r^2`$, is not defined ).

The VA starts to deviate from the exact results when the system is close to the collapse conditions and particularly for the unstable solutions (maxima) of the total energy. For this regions, where the Gaussian VA fails, one should improve the ansatz or take the results as a qualitative picture to guide the exact numerical calculations. In this respect, one should note that, for the NLSE without the trap potential, it is well known that the time-dependent variational approach fails to describe the region near the collapse .

Considering that in this work we have shown the relevance of an initial chirp parameter to study the stability of the condensate, we should observe that a good estimate for such parameter relies in a good approach to determining both minimum (stable) and maximum (unstable) solutions for the total energy. The Gaussian VA has a good agreement (quantitatively) with exact results for the minima, but can only give a qualitative description for the maxima. Still we can observe that the VA calculations of the chirp parameter, when compared with exact numerical calculations, are in reasonable agreement.

We also would like to point out that the present approach can be extended to study the stability of chirped laser beams in inhomogeneous three dimensional media with Kerr nonlinearity.

Finally, the main results of the present work are: (i) A chirp parameter in the wave-function is shown to be useful to study the initial conditions for a stable condensate to remain stable. (ii) we show that the VK criterion cannot be extended to cases of harmonic trapped BEC when the nonlinear two and three-body terms have opposite signs; (iii) if an atomic system has repulsive two-body interaction, the collapse is possible if the effective three-body interaction is negative. In this perspective, one can use the observed critical number of atoms in order to determine the corresponding three-body parameter.

## Acknowledgments

We are grateful to Fundação de Amparo à Pesquisa do Estado de São Paulo (FAPESP-Brazil) for partial financial support. LT and TF also acknowledge partial support from Conselho Nacional de Desenvolvimento Científico e Tecnológico.