Equilibrium and oscillations in a turbulent Bose–Einstein condensate

Turbulence is an important and fascinating topic in both classical and quantum fluids. This includes superfluid helium [1] and Bose–Einstein condensates (BEC) of low density alkaline gases, where turbulence has been studied in recent experiments [2]. From observation and simulations [3], it emerges the view that a turbulent condensate includes vortices and phonon oscillations.

Vortices have been observed in [4, 5], and the corresponding quantum theory was explored by a large number of authors [6, 7]. Large vortices of the size of the condensate, reminding those observed in [4], have also been theoretically described as Rossby waves, similar to those occurring in planet atmospheres [8]. Vortex condensates are known to organize in Abrikosov lattices, which have normal modes of oscillation called Tkachenko modes [9, 10]. These modes have been first observed by [11]. A more general class of lattice oscillations, associated with Rossby–Tkachenco modes, has also been considered [12].

Phonons are elementary excitations in the medium, and it is commonly believed that they cannot carry any angular momentum in uniform condensates. However, recent work has shown that twisted phonons, with angular momentum, can exist in a non-rotating and unbounded condensate [13]. Similar concepts have also been proposed, such as twisted phonons in classical fluids [14, 15].

Here we propose a general approach to oscillating BEC, where the mutual influence of quasi-static Thomas–Fermi (TF) equilibrium, vortices and elementary excitations of the turbulent spectrum is considered. The term turbulence is sometimes misleading as it can refer to different physical situations. Alternatively, we could talk about excited condensates, but this could also be misleading. Here we call turbulent a condensate where many different oscillations coexist. And we focus on a steady state turbulence, where the instabilities have saturated and the amplitude of the oscillations is nearly constant.

The spectrum of oscillations is assumed arbitrary, and described by a superposition of Bogoliubov–de Gennes (BdG) modes, in different geometrical configurations. BdG modes can be seen as generalized elementary oscillations of the condensate, which reduce to the simple phonons in a uniform medium.

The cases of cylindrical and toroidal media with arbitrary density profiles will be considered. The following questions will be addressed: (a) what are the changes introduced by turbulence on the quasi-static density profiles; (b) what are the changes introduced by such modified profiles on the basic BdG modes; (c) can these BdG modes include the case of twisted phonons, previously derived with a fluid model?

The structure of the paper is the following. In section 2, we start from the Gross–Pitaevskii (GP) equation, which efficiently describes the behavior of BECs in the mean-field approximation, and derive two different equations for the fast and quasi-static wavefunction components. These will be the basic equations of our model. In this introductory section it will be also shown that, for steady state turbulence, a modified TF equilibrium will eventually emerge, where the condensate density will adjust, in order to keep the unperturbed value of the chemical potential. Modified vortex solutions can also be derived.

In section 3, we successively consider the BdG solutions for an infinite BEC, and for cylindrical and toroidal configurations with uniform density. We first consider phonon beams with a beam waist much smaller than the transverse dimensions of the condensate, therefore justifying the infinite medium assumption and the neglect of boundary conditions. We show that the resulting modes are described by Laguerre–Gauss (LG) orthogonal functions and possess angular momentum. These are the expected twisted phonons, previously described using quantum fluid equations [13]. The cylindrical case is formally very similar, apart from the fact that the transverse boundary conditions are included, and a Fourier–Bessel (FB) decomposition is used instead. As for the toroidal case, it can be reduced to a periodic cylinder, if the major radius of the torus is much larger than the minor radius, which is relevant for experiments.

In section 4, we consider the case of a non-homogenous BEC, under quite general density profile conditions. Modified BdG modes will be derived, extending the results of the previous section. A general form of dispersion relation is obtained, which includes the effects of both the non-uniformity of the medium and the BdG mode configuration. Finally, in section 5, we consider nonlinear mode solutions and show that dark envelope solitons can occur for LG and FB modes. This is a new kind of soliton solution, which are not to be confused with the well known (non-envelope) dark solitons [12]. The present work confirms and generalizes our previous work based on a fluid description [13]. Finally, in section 6, we state some conclusions.

We start with the GP equation, describing the behavior of a condensate in the mean field approximation, using the standard notation, as

$$\begin{eqnarray}\begin{array}{rcl} {\rm i}\hbar \frac{\partial }{\partial t}\psi & = & \left( {{H}_{0}}+g{{\left| \psi \right|}^{2}} \right)\psi , \\ {{H}_{0}} & = & -\frac{{{\hbar }^{2}}{{\nabla }^{2}}}{2M}+{{V}_{0}}({\bf r}). \\ \end{array}\end{eqnarray} \tag{ 1 }$$

Here, $\psi \equiv \psi ({\bf r},t)$ is the collective wavefunction of the condensed matter, ${{V}_{0}}({\bf r})$ is the confining potential, and g is the coupling parameter associated with binary atomic collisions in the low energy limit. We assume that this equation is satisfied by a generic solution of the form

$$\begin{eqnarray}&&\psi ({\bf r},t)=\left[ {{\Psi }_{0}}({\bf r})+\tilde{\psi }({\bf r},t) \right]{\rm exp} (-{\rm i}\mu t/\hbar ),\end{eqnarray} \tag{ 2 }$$

where ${{\Psi }_{0}}({\bf r})$ describes the static equilibrium profile, $\tilde{\psi }$ is a superposition of BdG modes to be specified later, and μ is the chemical potential. Replacing this in equation (1), and taking a time average over a period much longer than the typical BdG mode period, such that $\left\langle \tilde{\psi } \right\rangle =0$, we obtain

$$\begin{eqnarray}&&\mu {{\Psi }_{0}}=\left[ {{\mu }_{0}}+2\;g\left\langle |\tilde{\psi }{{|}^{2}} \right\rangle \right]{{\Psi }_{0}},\end{eqnarray} \tag{ 3 }$$

where ${{\mu }_{0}}$ is determined by the unperturbed GP equation

$$\begin{eqnarray}&&{{\mu }_{0}}{{\Psi }_{0}}=\left[ {{H}_{0}}+g{{n}_{0}}({\bf r}) \right]{{\Psi }_{0}}.\end{eqnarray} \tag{ 4 }$$

The chemical potential μ introduced in equation (3) is clearly distinct from the unperturbed chemical potential, ${{\mu }_{0}}$, which is defined for a condensate in equilibrium and in the absence of oscillations [6]. On the other hand, ${{n}_{0}}({\bf r})=|{{\Psi }_{0}}({\bf r}){{|}^{2}}$ is the density profile. Assuming that equation (3) is satisfied, we can then use equation (1) to establish an evolution equation for the turbulent part of the wavefunction, of the form

$$\begin{eqnarray}&&{\rm i}\hbar \frac{\partial }{\partial t}\tilde{\psi }=\left( {{H}_{r}}+g{{\left| \tilde{\psi } \right|}^{2}} \right)\tilde{\psi }+\mathcal{L}({{\Psi }_{0}},\tilde{\psi }).\end{eqnarray} \tag{ 5 }$$

We can see that this equation strongly resembles the initial GP equation (1), but with a different Hamiltonian H_r, instead of H₀, and an additional term $\mathcal{L}({{\Psi }_{0}},\tilde{\psi })$, describing the interaction between the slow and the oscillating parts of the wavefunction. They are defined as

$$\begin{eqnarray}\begin{array}{rcl} {{H}_{r}} & = & {{H}_{0}}+2\;g{{n}_{0}}({\bf r})-\mu , \\ \mathcal{L}({{\Psi }_{0}},\tilde{\psi }) & = & g\left( \Psi _{0}^{2}{{\tilde{\psi }}^{*}}+\Psi _{0}^{*}{{\tilde{\psi }}^{2}} \right). \\ \end{array}\end{eqnarray} \tag{ 6 }$$

Equations (2) and (5) will be used as the starting point of our present work.

Before discussing the different BdG modes which are the basic elements of the turbulent spectrum, let us consider the modified equilibrium. For that purpose, we assume an arbitrary spectrum of oscillations, leading to an average density perturbation $\tilde{n}({\bf r})=2\left\langle |\tilde{\psi }{{|}^{2}} \right\rangle$. We can assume that a modified TF equilibrium will eventually be established, such that

$$\begin{eqnarray}&&{{n}_{0}}({\bf r})=\frac{1}{g}\left[ \mu -{{V}_{0}}({\bf r}) \right]-\tilde{n}({\bf r}).\end{eqnarray} \tag{ 7 }$$

This equality is obtained by neglecting the kinetic term of the Hamiltonian H₀, as usual. By assuming a radial potential shape ${{V}_{0}}({\bf r})=M\omega _{0}^{2}{{r}^{2}}/2$, where ${{\omega }_{0}}$ is the frequency of the harmonic trap, and taking the plausible assumption that turbulence is absent at the boundary r = a, where ${{n}_{0}}=0$ (as assumed in figure 1), we conclude that $\mu ={{V}_{0}}(a)\equiv M\omega _{0}^{2}{{a}^{2}}/2$. This contrast with the unperturbed value ${{\mu }_{0}}={{V}_{0}}({{a}_{0}})$, where a₀ is the radius of the condensate in the absence of turbulence. We notice that, due to the existence of a spectrum of oscillations, the size of the condensate increases, $a\gt {{a}_{0}}$. Integration over the BEC volume, assuming that the total number of condensed atoms is the same, leads to

$$\begin{eqnarray}&&{{a}^{4}}=a_{0}^{4}+\frac{4g}{\pi LM\omega _{0}^{2}}\int \tilde{n}({\bf r}){\rm d}{\bf r}.\end{eqnarray} \tag{ 8 }$$

This expression is valid for a cylindrical BEC with length L. The resulting modified density profile is illustrated in figure 1, for the cylindrical case in the presence of two BdG modes (these modes will be described below). In this figure we compare the dynamically renormalized profile with the unperturbed profile, as it would occur in a cylindrical nonuniform condensate. The observed static undulations of the density profile, the reduction of its maximum value ${{n}_{0}}(0)$, and the increase of radial dimension are due to the existence of oscillations.

Zoom In Zoom Out Reset image size

From the above discussion we can also conclude that the vortex equation for a modified TF equilibrium, will reduce to the usual form $\mu {{\Psi }_{0}}=\left( {{H}_{0}}+g{{\left| {{\Psi }_{0}} \right|}^{2}} \right){{\Psi }_{0}}$, where μ is now replacing ${{\mu }_{0}}$. No changes in the local structure of the vortex solutions will therefore be expected, except for the fact that vortices will exist on a modified density profile.

We now focus on the BdG modes, starting with the simple case of a uniform BEC, where the unperturbed density n₀ is assumed constant. The discussion of the elementary excitations of a BEC can be done by linearizing equation (5) with respect to the perturbation associated with a given turbulence mode $\tilde{\psi }$. We then get

$$\begin{eqnarray}&&{\rm i}\hbar \frac{\partial }{\partial t}\tilde{\psi }={{H}_{r}}\tilde{\psi }+g\Psi _{0}^{2}{{\tilde{\psi }}^{*}}.\end{eqnarray} \tag{ 9 }$$

For a mode with frequency ω, we can assume a solution of the form

$$\begin{eqnarray}&&\tilde{\psi }({\bf r},t)=u({\bf r}){\rm exp} (-{\rm i}\omega t)+{{v}^{*}}({\bf r}){\rm exp} ({\rm i}\omega t).\end{eqnarray} \tag{ 10 }$$

Replacing this in equation (9), we obtain the two equalities

$$\begin{eqnarray}&&(\hbar \omega -{{H}_{r}})u={{G}_{0}}v,\quad (\hbar \omega +H_{r}^{*})v=-G_{0}^{*}u,\end{eqnarray} \tag{ 11 }$$

where we have used ${{G}_{0}}=g\Psi _{0}^{2}$. These are the so-called BdG equations for the condensate. At this point we assume a solution satisfying the equations ${{\nabla }^{2}}\;u\;=\;-{{k}^{2}}\;u\quad {\rm and}\quad {{\nabla }^{2}}v=-{{k}^{2}}v$, where k is a wavenumber, to be specified. The kinetic operator inside H_r is then reduced to ${{\hbar }^{2}}{{k}^{2}}/2M$, and the two equations in (10) become purely algebraic, and can be written in matricial form as ${\bf \bar{L}}\cdot {\bf a}=0$, with

$$\begin{eqnarray}&&{\bf a}=\left[ \begin{array}{ccccccccccccccc} u \\ v \\ \end{array} \right],\quad {\bf \bar{L}}=\left[ \begin{array}{ccccccccccccccc} (\hbar \omega -{{H}_{r}})\quad \quad -{{G}_{0}} \\ G_{0}^{*}\quad \quad \quad (\hbar \omega +H_{r}^{*}) \\ \end{array} \right].\end{eqnarray} \tag{ 12 }$$

Solution of this matrix equation implies that ${\rm Det}.({\bf \bar{L}})=0$, which can be written as

$$\begin{eqnarray}&&\hbar \omega =\sqrt{|{{H}_{r}}{{|}^{2}}-{{\left| {{G}_{0}} \right|}^{2}}}=\sqrt{{{\left| {{H}_{0}}+2\;g{{n}_{0}}-\mu \right|}^{2}}-{{(g{{n}_{0}})}^{2}}}.\end{eqnarray} \tag{ 13 }$$

Ignoring the confining potential ${{V}_{0}}=0$, and assuming that the chemical potential is $\mu =g{{n}_{0}}$, this reduces to the well known dispersion relation [6]

$$\begin{eqnarray}&&{{\omega }^{2}}=c_{s}^{2}{{k}^{2}}+\frac{{{\hbar }^{2}}{{k}^{4}}}{4{{M}^{2}}},\quad {{c}_{s}}=\sqrt{\frac{g{{n}_{0}}}{M}},\end{eqnarray} \tag{ 14 }$$

where c_s is the Bogoliubov sound velocity. In order to understand the generality of this dispersion relation, we consider the case of twisted phonons [13] where, instead of a plane wave, we have a wave beam propagating along the z-axis, as $u({\bf r})=\tilde{u}({\bf r}){\rm exp} ({\rm i}kz),\quad v({\bf r})=\tilde{v}({\bf r}){\rm exp} ({\rm i}kz)$. Here, $\tilde{u}$ and $\tilde{v}$ are slowly varying amplitudes, such that $\left| {{\partial }^{2}}\tilde{u}/\partial {{z}^{2}} \right|\ll \left| k\;(\partial \tilde{u}/\partial z) \right|$ and $\left| {{\partial }^{2}}\tilde{v}/\partial {{z}^{2}} \right|\ll \left| k\;(\partial \tilde{v}/\partial z) \right|$ . They satisfy the two paraxial equations

$$\begin{eqnarray}&&\left( \nabla _{\bot }^{2}+2{\rm i}k\frac{\partial }{\partial z} \right)u=0,\quad \left( \nabla _{\bot }^{2}+2{\rm i}k\frac{\partial }{\partial z} \right)v=0.\end{eqnarray} \tag{ 15 }$$

Solutions of these equations are well known, and can be written as a superposition of LG modes [16]

$$\begin{eqnarray}&&u({\bf r})={{\tilde{u}}_{pl}}{{F}_{pl}}({\bf r}),\quad v({\bf r})={{\tilde{v}}_{pl}}{{F}_{pl}}({\bf r}),\end{eqnarray} \tag{ 16 }$$

where ${{\tilde{u}}_{pl}}$ and ${{\tilde{v}}_{pl}}$ are constant amplitudes, and ${{F}_{pl}}({\bf r})$ are LG functions of the form

$$\begin{eqnarray}&&{{F}_{pl}}({\bf r})={{C}_{pl}}L_{p}^{|l|}(X){{X}^{|l|}}{\rm exp} \left( {\rm i}l\theta -\frac{X}{2} \right).\end{eqnarray} \tag{ 17 }$$

Here, C_pl are normalization constants, L_p^l are the associated Laguerre polynomials, and $X={{r}^{2}}/{{w}^{2}}(z)$, where the wave beam waist $w(z)$ slowly varies along the z-axis. It is well known that, choosing appropriate normalization constants, these LG functions satisfy the orthogonality relations

$$\begin{eqnarray}&&\int _{0}^{\infty }r{\rm d}r\int _{0}^{2\pi }{\rm d}\theta {{F}_{pl}}{{({\bf r})}^{*}}{{F}_{p^{\prime} l^{\prime} }}({\bf r})={{\delta }_{pp^{\prime} }}{{\delta }_{ll^{\prime} }}.\end{eqnarray} \tag{ 18 }$$

Such solutions correspond to twisted phonons, carrying a finite amount of angular momentum, as first derived in [13] by using the quantum fluid equations. These modes are re-derived here, as a special example of BdG modes. They correspond to phonons with a beam waist much shorter than the size of the condensate.

Let us now consider the opposite situation, where the wave structure occupies the whole size of the condensate, and the finite boundaries play a significant role. For this purpose, we consider a condensate with a homogeneous cylindrical shape, with radius a. In this case, we can assume solutions of the form

$$\begin{eqnarray}&&u({\bf r})=\tilde{u}(r){\rm exp} ({\rm i}l\theta +{\rm i}\beta z),v({\bf r})=\tilde{v}(r){\rm exp} ({\rm i}l\theta +{\rm i}\beta z).\end{eqnarray} \tag{ 19 }$$

We then obtain differential equations for $\tilde{u}(r)$ and $\tilde{v}(r)$ of the form

$$\begin{eqnarray}&&\left[ r\frac{{\rm d}}{{\rm d}r}\left( r\frac{{\rm d}}{{\rm d}r} \right)+\left( k_{\bot }^{2}{{r}^{2}}-{{l}^{2}} \right) \right]\tilde{u}=0,\end{eqnarray} \tag{ 20 }$$

with $k_{\bot }^{2}={{k}^{2}}-{{\beta }^{2}}$, and a similar equation for $\tilde{v}$. The solutions are Bessel functions of order l, such that

$$\begin{eqnarray}&&\tilde{u}(r)={{\tilde{u}}_{l}}{{J}_{l}}({{k}_{\bot }}r),\quad \tilde{v}(r)={{\tilde{v}}_{l}}{{J}_{l}}({{k}_{\bot }}r),\end{eqnarray} \tag{ 21 }$$

with constant amplitudes ${{\tilde{u}}_{l}}$ and ${{\tilde{v}}_{l}}$. Noting that the wave amplitudes vanish at the boundary r = a, we have to choose the perpendicular wavenumber ${{k}_{\bot }}$ such that ${{J}_{l}}({{k}_{\bot }}a)=0$, and $\quad {{k}_{\bot }}={{\alpha }_{lp}}/a$, where ${{\alpha }_{lp}}$ are the zeros of the Bessel function J_l. It can easily be shown that the solutions for these new BdG modes, valid for a cylindrical configuration, satisfy orthogonality relations similar to equation (18), of the form

$$\begin{eqnarray}&&\int _{0}^{a}r{\rm d}r\int _{0}^{2\pi }{\rm d}\theta {{\mathcal{F}}_{pl}}{{(r,\theta )}^{*}}{{\mathcal{F}}_{p^{\prime} l^{\prime} }}(r,\theta )={{\delta }_{pp^{\prime} }}{{\delta }_{ll^{\prime} }},\end{eqnarray} \tag{ 22 }$$

where we have defined the new orthogonal FB functions

$$\begin{eqnarray}&&{{\mathcal{F}}_{pl}}(r,\theta )={{C}_{pl}}{{J}_{l}}({{\alpha }_{pl}}r/a){\rm exp} \left( l\theta \right).\end{eqnarray} \tag{ 23 }$$

The new normalization constant is ${{C}_{pl}}=a\sqrt{\pi }{{J}_{l+1}}({{\alpha }_{pl}})$. Finally, it should be noticed that these cylindrical BdG solutions can easily be adapted to a toroidal configuration, if the large radius of the torus R, is much larger than the minor radius a. We can then describe the torus as a periodic cylinder, with $z=R\phi$, where ϕ is toroidal angle. For the resulting toroidal BdG modes, the above cylindrical solutions stay valid, but only discrete values of the parallel wavenumber are considered, such that $\beta ={{\beta }_{m}}=m/R$, where the integer m is the toroidal mode number.

Let us now generalize the above description for the case of a non-uniform condensate. We assume that the condensate density is slowly varying (and eventually periodic) along the z-axis. The transverse density profile, ${{n}_{0}}({\bf r})$ is assumed arbitrary and not necessarily cylindrically symmetric. Instead of the above simple assumption, stated in equation (10), we now use a generic mode decomposition of the form

$$\begin{eqnarray}\begin{array}{rcl} \tilde{\psi }({\bf r},t) & = & \mathop{\sum }\limits_{k}\left[ {{u}_{k}}{{F}_{pl}}({{{\bf r}}_{\bot }}){\rm exp} \left( {\rm i}\int {{k}_{z}}(z){\rm d}z-{\rm i}{{\omega }_{k}}t \right) \right. \\ {} & {} & +\left. v_{k}^{*}F_{pl}^{*}({{{\bf r}}_{\bot }}){\rm exp} \left( -{\rm i}\int {{k}_{z}}(z){\rm d}z+{\rm i}{{\omega }_{k}}t \right) \right], \\ \end{array}\end{eqnarray} \tag{ 24 }$$

where the subscript k represents here the set of mode numbers, $(l,p,{{k}_{z}})$, where ${{k}_{z}}=k$ for LG modes, ${{k}_{z}}=\beta$ for cylindrical modes, or ${{k}_{z}}=m/R$ for toroidal modes. Using the modified TF equilibrium of equation (7), such that ${{H}_{r}}=-({{\hbar }^{2}}{{\nabla }^{2}}/2M)+g{{n}_{0}}({\bf r})$, we can get from the linearized wave equation (9) the following pair of BdG equations, for each k mode

$$\begin{eqnarray}\begin{array}{rcl} \left[ \hbar {{\omega }_{k}}+\frac{{{\hbar }^{2}}{{\nabla }^{2}}}{2M}-g{{n}_{0}}({\bf r}) \right]{{u}_{k}}{{F}_{pl}} & = & {{G}_{0}}{{v}_{k}}{{F}_{pl}}, \\ \left[ \hbar {{\omega }_{k}}-\frac{{{\hbar }^{2}}{{\nabla }^{2}}}{2M}+g{{n}_{0}}({\bf r}) \right]{{v}_{k}}{{F}_{pl}} & = & -G_{0}^{*}{{u}_{k}}{{F}_{pl}}. \\ \end{array}\end{eqnarray} \tag{ 25 }$$

At this point we should notice that, each of these k modes satisfy an equation of the form $({{\nabla }^{2}}+{{k}^{2}}){{u}_{k}}{{F}_{pl}}=0$ and $({{\nabla }^{2}}+{{k}^{2}}){{v}_{k}}{{F}_{pl}}=0$. Using the relevant orthogonality conditions (18) or (22), and after integration over ${{{\bf r}}_{\bot }}=(r,\theta )$, we obtain the following result

$$\begin{eqnarray}\begin{array}{rcl} \left( \hbar {{\omega }_{k}}-\frac{{{\hbar }^{2}}{{k}^{2}}}{2M}-{{G}_{k}} \right){{u}_{k}} & = & {{S}_{k}}{{v}_{k}}, \\ \left( \hbar {{\omega }_{k}}+\frac{{{\hbar }^{2}}{{k}^{2}}}{2M}+{{G}_{k}} \right){{v}_{k}} & = & -S_{k}^{*}{{u}_{k}}, \\ \end{array}\end{eqnarray} \tag{ 26 }$$

where we have defined

$$\begin{eqnarray}\begin{array}{rcl} {{G}_{k}} & = & g{{\int }_{\bot }}{{n}_{0}}({\bf r}){{\left| {{F}_{pl}}({{{\bf r}}_{\bot }}) \right|}^{2}}\;{\rm d}{{{\bf r}}_{\bot }}, \\ {{S}_{k}} & = & g{{\int }_{\bot }}\Psi _{0}^{2}({\bf r}){{\left| {{F}_{pl}}({{{\bf r}}_{\bot }}) \right|}^{2}}\;{\rm d}{{{\bf r}}_{\bot }}. \\ \end{array}\end{eqnarray} \tag{ 27 }$$

Here, these quantities can slowly vary along z, and the integration over ${{{\bf r}}_{\bot }}$ extends from r = 0 to r = R, for the cylindrical case, or from r = 0 to infinity for LG modes. These new equations (26) can be re-written in matrix form ${{{\bf \bar{L}}}_{k}}\cdot {{{\bf a}}_{k}}=0$, with

$$\begin{eqnarray}&&{{{\bf a}}_{k}}=\left[ \begin{array}{ccccccccccccccc} {{u}_{k}} \\ {{v}_{k}} \\ \end{array} \right],\quad {{{\bf \bar{L}}}_{k}}=\left[ \begin{array}{ccccccccccccccc} (\hbar {{\omega }_{k}}-{{H}_{k}})\quad \quad -{{S}_{k}} \\ S_{k}^{*}\quad \quad \quad (\hbar {{\omega }_{k}}+{{H}_{k}}) \\ \end{array} \right],\end{eqnarray} \tag{ 28 }$$

where we have defined ${{H}_{k}}=({{\hbar }^{2}}{{k}^{2}}/2M)+{{G}_{k}}$. This leads to the general dispersion relation of BdG in a non-uniform medium, of the form

$$\begin{eqnarray}&&\hbar {{\omega }_{k}}=\sqrt{{{\left( \frac{{{\hbar }^{2}}{{k}^{2}}}{2M}+{{G}_{k}} \right)}^{2}}-{{\left| {{S}_{k}} \right|}^{2}}}.\end{eqnarray} \tag{ 29 }$$

From here, we obtain

$$\begin{eqnarray}&&\omega _{k}^{2}=\omega _{0k}^{2}+c_{sk}^{2}{{k}^{2}}+\frac{{{\hbar }^{2}}{{k}^{4}}}{4{{M}^{2}}},\end{eqnarray} \tag{ 30 }$$

where the cut-off frequency ${{\omega }_{0k}}$, and the modified Bogoliubov sound velocityc_sk, are defined by

$$\begin{eqnarray}&&\omega _{0k}^{2}=G_{k}^{2}-{{\left| {{S}_{k}} \right|}^{2}},\quad c_{sk}^{2}=\frac{{{G}_{k}}}{M}.\end{eqnarray} \tag{ 31 }$$

This is a generalization of the dispersion relation equation (14) to a non-uniform medium. At this point, it should be notice that, in the definition of S_k of equation (27), we can write $\Phi _{0}^{2}({\bf r})={{n}_{0}}({\bf r}){\rm exp} ({\rm i}\varphi )$. This means that S_k only differs from G_k because of the phase factor ${\rm exp} ({\rm i}\varphi )$ inside the integrand. If this phase is assumed constant, we will have $G_{k}^{2}={{\left| {{S}_{k}} \right|}^{2}}$, and the cut-off frequency ${{\omega }_{0k}}$ will vanish. However, this is not necessarily true. It will depend on the solutions of the average wave equation (3).

We conclude that, in general, the BdG modes show a cut-off frequency and a modified sound velocity, which both depend on the shape of the condensate density, and on the shape of the BdG mode. Both quantities will also eventually depend on z, if the condensate density is supposed to vary slowly along the axis of propagation. In this case, the total mode wavenumber k will slowly adapt to the local axial conditions, for a fixed mode frequency. On the other hand, the transverse shape of the equilibrium density can be assumed as arbitrary.

Let us now consider a nonlinear mode analysis. For that purpose, we use a single mode solution, with frequency ${{\omega }_{k}}$, of the form

$$\begin{eqnarray}\begin{array}{rcl} \tilde{\psi }({\bf r},t) & = & {{u}_{k}}(z,t){{F}_{pl}}({\bf r}){\rm exp} (-{\rm i}{{\omega }_{k}}t) \\ {} & {} & +v_{k}^{*}(z,t)F_{pl}^{*}({\bf r}){\rm exp} ({\rm i}{{\omega }_{k}}t). \\ \end{array}\end{eqnarray} \tag{ 32 }$$

Replacing this in the nonlinear wave equation, equation (5), and neglecting the terms describing second and third harmonics evolving as $(\pm 2{\rm i}{{\omega }_{k}}t)$ and $(\pm 3{\rm i}{{\omega }_{k}}t)$, we obtain the following expressions

$$\begin{eqnarray}\begin{array}{rcl} \hbar \left( {{\omega }_{k}}+{\rm i}\frac{\partial }{\partial t} \right){{u}_{k}}{{F}_{pl}} & = & \left( {{H}_{{\rm eff}}}+g|{{v}_{k}}{{|}^{2}}|{{F}_{pl}}{{|}^{2}} \right){{u}_{k}}{{F}_{pl}} \\ {} & {} & +{{G}_{0}}{{v}_{k}}{{F}_{pl}}\; \\ \end{array}\end{eqnarray} \tag{ 33 }$$

and

$$\begin{eqnarray}\begin{array}{rcl} -\hbar \left( {{\omega }_{k}}-{\rm i}\frac{\partial }{\partial t} \right){{v}_{k}}{{F}_{pl}} & = & \left( {{H}_{{\rm eff}}}+g|{{u}_{k}}{{|}^{2}}|{{F}_{pl}}{{|}^{2}} \right){{v}_{k}}{{F}_{pl}} \\ {} & {} & +G_{0}^{*}{{u}_{k}}{{F}_{pl}}, \\ \end{array}\end{eqnarray} \tag{ 34 }$$

where we have used the auxiliary quantity

$$\begin{eqnarray}&&{{H}_{{\rm eff}}}={{H}_{r}}+g\left( |{{u}_{k}}{{|}^{2}}+|{{v}_{k}}{{|}^{2}} \right)|{{F}_{pl}}{{|}^{2}}.\end{eqnarray} \tag{ 35 }$$

We now assume that the mode amplitudes depend on z, according to $({{u}_{k}},{{v}_{k}})=({{\tilde{u}}_{k}},{{\tilde{v}}_{k}}){\rm exp} ({\rm i}{{k}_{z}}z)$, where the new amplitudes also depend on z and t. This means that the nabla-square operator inside the Hamiltonian H_r can be replaced by

$$\begin{eqnarray}&&{{\nabla }^{2}}({{u}_{k}}{{F}_{pl}})={{F}_{pl}}{\rm exp} ({\rm i}{{k}_{z}}z)\left( {{k}^{2}}+2{\rm i}{{k}_{z}}\frac{\partial }{\partial z}+\frac{{{\partial }^{2}}}{\partial {{z}^{2}}} \right){{\tilde{u}}_{k}},\end{eqnarray} \tag{ 36 }$$

with a similar equation for v_k. Notice that, with an appropriate choice of orthogonal functions F_pl, these expressions are valid for LG modes, (if we assume ${{k}_{z}}\equiv k$), as well as for cylindrical Bessel modes with radius a (if we assume ${{k}_{z}}=\beta$ and ${{k}_{\bot }}=({{\alpha }_{lp}}/a)$, with ${{k}^{2}}=k_{\bot }^{2}+k_{z}^{2}$). The toroidal configuration with major radius $R\gg a$, can also be described by the same type of cylindrical Bessel modes, if in addition we use discrete values for the axial wavenumber ${{k}_{z}}=m/R$, as discussed above.

Taking this into account, multiplying equation (33) by $F_{pl}^{*}$, and integrating over the transverse space variables, we obtain the nonlinear BdG relations

$$\begin{eqnarray}\begin{array}{rcl} \left[ \hbar \left( {{\omega }_{k}}+{\rm i}\frac{\partial }{\partial t} \right)-{{H}_{k}} \right]{{\tilde{u}}_{k}} & = & {{S}_{k}}{{\tilde{v}}_{k}}+{{R}_{k}}|{{\tilde{v}}_{k}}{{|}^{2}}{{\tilde{u}}_{k}}, \\ \left[ -\hbar \left( {{\omega }_{k}}-{\rm i}\frac{\partial }{\partial t} \right)-{{H}_{k}} \right]{{\tilde{u}}_{k}} & = & S_{k}^{*}{{\tilde{u}}_{k}}+{{R}_{k}}|{{\tilde{u}}_{k}}{{|}^{2}}{{\tilde{v}}_{k}}. \\ \end{array}\end{eqnarray} \tag{ 37 }$$

Here, we have used the nonlinear Hamiltonian operator

$$\begin{eqnarray}\begin{array}{rcl} {{\mathcal{H}}_{k}} & = & \left( \frac{{{\hbar }^{2}}{{k}^{2}}}{2M}+{{G}_{k}} \right)-\frac{{{\hbar }^{2}}}{2M}\left( 2{\rm i}{{k}_{z}}\frac{\partial }{\partial {{k}_{z}}}+\frac{{{\partial }^{2}}}{\partial {{z}^{2}}} \right) \\ {} & {} & +\left( |{{\tilde{u}}_{k}}{{|}^{2}}+|{{\tilde{v}}_{k}}{{|}^{2}} \right){{R}_{k}}. \\ \end{array}\end{eqnarray} \tag{ 38 }$$

We recognize in the first two terms under brackets the linear Hamiltonian H_k, appearing in equation (28). The quantities G_k and S_k are defined as in the previous section, and we have added a new quantity

$$\begin{eqnarray}&&{{R}_{k}}=g{{\int }_{\bot }}|{{F}_{pl}}({\bf r}){{|}^{4}}\;{\rm d}{{{\bf r}}_{\bot }}.\end{eqnarray} \tag{ 39 }$$

At this point, we should use symmetry arguments, and assume $|{{\tilde{u}}_{k}}{{|}^{2}}=|{{\tilde{v}}_{k}}{{|}^{2}}$, as an additional and necessary condition. Assuming that the linear dispersion relation for the generic BdG modes, as given by equations (29)–(31), is still satisfied, we can then reduce the system of coupled equations (37) to the following nonlinear evolution

$$\begin{eqnarray}&&{\rm i}\hbar \frac{\partial }{\partial t}{\bf u}+\frac{{{\hbar }^{2}}}{2M}\left( 2{\rm i}{{k}_{z}}\frac{\partial }{\partial {{k}_{z}}}+\frac{{{\partial }^{2}}}{\partial {{z}^{2}}} \right){\bf u}-\frac{3}{2}{{R}_{k}}{{\left| {\bf u} \right|}^{2}}{\bf u}=0,\end{eqnarray} \tag{ 40 }$$

where the vector ${\bf u}$ is defined by ${\bf u}=\left[ {{\tilde{u}}_{k}},{{\tilde{v}}_{k}} \right]$. It is now convenient to use a transformation of variables, from $(z,t)$ to the non-dimensional variables $(\xi ,\tau )$, as defined by

$$\begin{eqnarray}&&\xi ={{k}_{z}}\left( z-\frac{\hbar {{k}_{z}}}{M}t \right),\quad \tau ={{\omega }_{k}}t.\end{eqnarray} \tag{ 41 }$$

This transforms equation (40) into the form

$$\begin{eqnarray}&&{\rm i}\frac{\partial {\bf u}}{\partial \tau }+\frac{\hbar k_{z}^{2}}{2\;M{{\omega }_{k}}}\frac{{{\partial }^{2}}{\bf u}}{\partial {{\xi }^{2}}}-\frac{3}{2}\frac{{{R}_{k}}}{\hbar {{\omega }_{k}}}{{\left| {\bf u} \right|}^{2}}{\bf u}=0.\end{eqnarray} \tag{ 42 }$$

This is a defocusing nonlinear Shrödinger equation, which can be rewritten in a more appropriate form by rescaling the space variable ξ and the mode amplitude ${\bf u}$, as

$$\begin{eqnarray}&&\eta ={{\left( \frac{M{{\omega }_{k}}}{\hbar k_{z}^{2}} \right)}^{1/2}}\xi ,\quad {\bf u}={{\left( \frac{2\hbar {{\omega }_{k}}}{3{{R}_{k}}} \right)}^{1/2}}{\bf A}.\end{eqnarray} \tag{ 43 }$$

We are then led to the final standard form

$$\begin{eqnarray}&&{\rm i}\frac{\partial {\bf A}}{\partial \tau }+\frac{1}{2}\frac{{{\partial }^{2}}{\bf A}}{\partial {{\eta }^{2}}}-{{\left| {\bf A} \right|}^{2}}{\bf A}=0.\end{eqnarray} \tag{ 44 }$$

It is well known that such an equation satisfies dark soliton solutions, given by

$$\begin{eqnarray}&&{\bf A}={{{\bf A}}_{0}}{\rm tanh} \left( |{{{\bf A}}_{0}}|\eta \right){\rm exp} \left( {\rm i}|{{{\bf A}}_{0}}{{|}^{2}}\tau \right),\end{eqnarray} \tag{ 45 }$$

where the constant vector ${{{\bf A}}_{0}}$ determines the soliton amplitude. Such solutions describe nonlinear BdG modes in a non-uniform medium, and correspond to twisted phonons propagating in a non-uniform BEC. They generalize similar envelope soliton solutions [13], previously derived from a quantum fluid description. The present approach is, not only more general because it applies for non-uniform condensates, but also more exact in the sense that the GP wave equation is used, instead of the fluid approximation. Furthermore, these new solutions were built over a modified TF equilibrium profile, which is attained in the presence of an arbitrary broadband spectrum of oscillations, assumed with constant amplitudes, as described in section 2. The density perturbations associated with these dark envelope solitons, as given by $\tilde{n}=(1-|A{{|}^{2}})|{{F}_{pl}}{{|}^{2}}$, are illustrated in figure 2.

Zoom In Zoom Out Reset image size

A general approach to turbulence in BEC was proposed in this work. Starting from the GP equation, which adequately describes the condensate in the mean-field approximation, we have derived two different equations describing the steady state and the elementary excitations in a BEC. We have shown that the TF equilibrium, as well as the vortex structure, can be modified due to the existence of a background spectrum of oscillations.

On the other hand, we have derived a general dispersion relation for BdG modes valid in different geometries, which takes into account the modified quasi-static equilibria. We have explicitly discussed the cases of infinite and finite condensates, with or without homogeneous density. We have shown that the twisted phonons considered by us in a recent publication [13] are indeed BdG modes, defined in the paraxial approximation. The cases of cylindrical and toroidal condensates were considered. Finally, we have studied the nonlinear BdG solutions, and shown that envelope dark solitons exist for modes satisfying the phonon dispersion relation. These soliton solutions confirm and generalize our previous results on twisted nonlinear waves.

The possible existence of global oscillations, such as breathing and surface modes, was not considered here. Such modes, are known in both condensed and non-condensed gas [17, 18], and can also play a role in the turbulent medium. They could be treated as a natural extension of the present work. We have also ignored transient phenomena, and have focused our attention on nearly steady state oscillations. An important problem would be the influence of the BdG mode spectrum on stability and growth of vortices [19], as well as the exchange of energy and angular momentum between these structures and the elementary excitations. This will left to a future publication. Another important extension of the present approach concerns coupling with a non-condensed background gas, where modified gap solitons, similar to those found in [20], could eventually be observed. Finally, the present work suggests that a quasi-particle description of the BEC spectrum of oscillations is possible, where the vortex structures and equilibria can be defined on a background of a quasi-particle gas. This can be seen as analogue to the well known two-fluids approach to superfluidity, and could contribute to the understanding of finite temperature effects in a Bose gas [21].

We would like to thank the Brazilian research funding agencies CNPq and FAPESP, for financial support.