Quantum gases of dipoles, quadrupoles and octupoles in Gross-Pitaevskii formalism with form factor

Classical and quantum field theory of dipolar, axisymmetric quadrupolar and octupolar Bose gases is considered within a general approach. Dipole, axisymmetric quadrupole and octupole interaction potentials in the momentum representation are calculated. These results clearly demonstrate attraction and repulsion areas in corresponding gases. Then the Gross-Pitaevskii (GP) equation, which plays a key role in the present paper, is derived from the corresponding functional. The zoology of the form factors appearing in GP equation is studied in details. The classes of proper for the description of spatially non-uniform condensates form factors are chosen. In the Thomas-Fermi approximation a general solution of the GP equation with a quasilocal form factor is obtained. This solution has an interesting form in terms of a double rapidly converging series that universally includes all the interactions considered. An important analysis of the condensate stability, in other words the study of condensate excitations, is also performed in this paper. In the Gaussian approximation (from the GP functional), a functional describing the perturbations of the condensate is derived in details. For a probe wave function in the form of a plane wave, a spectrum of Bogoliubov excitations was obtained, from which an equation describing the threshold momentum for the emergence of instability was derived. An important result of this paper is the dependence of the threshold on the momentum of a stationary condensate. For completeness of the presentation, the approximating expression in the form of a rapidly converging series is obtained for the corresponding dependence. Finally, in the conclusion a quantum hydrodynamic theory for dipolar, axisymmetric quadrupolar and octupolar gases is briefly presented, giving a clue to the experimental determination of the form factors.

Discussing the non-relativistic systems of many particles, we are faced with a colossally complex problem of solving the Schrödinger equation for an enormous number of particles. To date, standard methods for solving this problem are the methods of quantum field theory (see, for example, the monographs [1][2][3][4][5][6][7][8][9]): in terms of the functional integral over theorys primary fields, a generating functional of Green functions or n-particle (n ≥ 1) irreducible vertices (the corresponding functional is called an effective action) is formulated. Further calculation of this functional integral is carried out in terms of perturbation theory, which leads to the construction of a diagram technique (Feynman diagrams) for generating functionals, or immediately for the corresponding families of Green functions.
Since the expressions constructed using the diagram technique are asymptotic series, the next step in the obtaining of a meaningful answer is the application of summation methods for asymptotic series, for example, the Borel-Laplace method. However, the latter is more common among the quantum field community ( [6][7][8][9]). In a solid-state physics community, we usually try to sum the subsequences of the diagrams in such a way that the final expressions for the Green functions make sense, in other words, try to find a "re-expansion" with respect to some new effective coupling constant ( [1][2][3][4][5]).
Another method of calculating functional integrals is the saddle-point method. Within this method an effective (for example, already taking into account a series of diagrams of a certain type, or various nonperturbative effects, not "captured" by perturbation theory) equation of motion for a certain Green function (the saddle-point method is well described in [7]). An example of such an equation is the well-known Gross-Pitaevskii (GP) equation, repeatedly applied for description of quantum gases of bosons, fermions and their mixtures ( [10] and [11]). Thus, the GP equation has a fundamental microscopic meaning.
At the same time, we can find ourselves in the following situation: the saddle-point method may not "catch" the required behaviour of the quantum gas. Simultaneously, this behaviour may not be achieved within the diagram technique. These situations are common in literature, and one of the generally accepted methods is the functional (nonperturbative, exact) renormalization group method ( [8,9] and [14]). Alternatively, we can try to develop a semi-phenomenological model. For instance, the simplified (stochastic) quantum field model is widely used for the fundamental turbulence phenomenon description. Taking into account the experience of the semi-phenomenological description, a model based on the GP equation with a phenomenological form factor is proposed in the present paper (the analysis of the form factors in the context of quantum field theory is performed in [12] and [13]. What is known about such a "building block of the theory" from the most general considerations? It should correctly reproduce different distributions of several reference physical quantities in the system and give a qualitative picture. Therefore, there are no univocal rules for choosing this building block (the ambiguity in the choice can be reduced by restrictions due to the renormalization group, the Ward identities or the Schwinger-Dyson equations). In defence of such a vulnerable for criticism state, a number of generally accepted arguments is presented (see [7]).
In microscopic theories, these building blocks must be generated by various microscopic mechanisms, and their characteristics for a particular problem must be computable. However, if there is no microscopic theory of this kind, within an effective description, which is only a simplified semi-phenomenological version of the (hypothetical) accurate theory, a specific choice of the form factor can be justified by general considerations and results.
Further, let us provide a brief literature review reflecting on a current situation. Here it should be noted that any choice is subjective hence, by definition, incomplete. For example, it can be said that the quantum Bose gases science was revitalized because of the relatively recent series of experimental papers ( [15][16][17]) on the direct observation of the Bose-Einstein condensation phenomenon in various atomic gases in traps. A number of theoretical publications in this field has significantly increased since those experimental successes. Only in the last decade studies such as Bose-Einstein condensation in dipole systems ( [18,19]), different modes of electronhole pairing in certain graphene-based structures, in particular in graphene bilayer ( [20][21][22]), and similar pairing in very popular to date topological insulators thin films ( [23]) can be noted. Both theoretical and experimental situations are described in detail in the remarkable review ( [24]). Here we should also mention a very trustworthy theoretical work about exciton-roton excitations in two-dimensional Bose gases of quadrupoles ( [25]). It should be noted, that in this paper the transition from the three-dimensional case to the two-dimensional case is performed strictly and consistently.
Talking about the phenomenon of Bose-Einstein condensation, we must note a brief description of the possibility of obtaining a Bose-Einstein condensate from a more general phenomenon -the BCS-BEC crossover phenomenon. This phenomenon, valid in quantum Fermi gases, consists in changing of the behaviour of the corresponding gas upon the transition from superconducting Bardeen-Cooper-Schrieffer (BCS) pairing to Bose-Einstein condensate of molecules strongly localized in the coordinate space. Due to that we note reviews [26] and [27], and also an interesting review [28] devoted to a wide range of condensed matter physics phenomena using the example of ultracold atomic gases in optical traps. Finally, a very intriguing theoretical approach to describing various condensed matter physics phenomena is the so-called "holographic" approach, originated from fundamental papers [29][30][31]. A lot of publications are focused on this approach, among which we note papers [32,33], focused on holographic picture for d-wave superconductor. This concludes the analysis of the literature and we proceed to the description of the structure of our paper.
This paper consists of five main sections and one appendix. After the introduction there is "Multipole interactions" section, where classical interaction potentials in dipolar, axisymmetric quadrupolar and octupolar gases are calculated in coordinate r representation as well as in momentum k representation, also the expression for the l-pole momentum interaction is presented. We note that various modifications of this interaction are also discussed in our paper. In addition to the classical modification with core (a hard core illustrating a hard repulsion at short distances) an original modification with "split" is proposed. This mechanism implies the splitting of the interaction by a certain feature, for example, by the sign of the initial interaction range.
The next section, "The Gross-Pitaevskii equation", is entral to this paper. It starts with obtaining of the Gross-Pitaevskii equation, which models the interaction action as an effective result of calculation of the functional integral. The GP equation is derived both in the coordinate and the momentum representation. The key feature of this model is that it allows considering almost a general scenario of the spatially non-uniform condensate formation and its excitations (of an arbitrary origin). This semi-phenomenological model allows describing the physics of practically any quasiparticles arising in considered quantum Bose gases.
Then follows a detailed study of various types of form factors. It was considered that for a wide class of physical scenarios the so-called "quasilocal" form factors are sufficient. A general solution of the GP equation is obtained for this type of form factors in Thomas-Fermi approximation. This solution is one of the main results of this paper. This solution is presented in terms of a double rapidly converging series (the last statement can be proved directly by constructing the corresponding majorant). In addition to the above the obtained solution has a universal form in terms of the dimension of the space D (for specific calculations we use D = 3). The obtained solution also universally includes all types of Bose gases considered.
The results are demonstrated in the form of graphs of condensate density functions for the exponentialtrigonometric form factor, because the exponent has good properties for numerical calculations, and the trigonometric argument is perfect for reproducing of a stationary periodic structure in the system. For the sake of completeness the GP equation in the limit of small condensate densities is also considered, together with the case when the Bose gas is in a trap. The potential is chosen in a form of an optical lattice. The limit of small condensate densities does not distinguish between dipolar, quadrupolar and octupolar gases, because the nonlinear term is equal to zero in this limit. The linear limit of the GP equation with a potential of an optical lattice is a single-particle Schrödinger equation in a periodic potential, and its stationary analogue is the Mathieu equation. Analytical results for the energy spectrum of this quantum mechanical system are presented in this paper in the limit of weak coupling as well as strong coupling in terms of the parameters of the trap (with reference to [34] and [35]).
The study of the properties of the stationary solution of the GP equation cannot be considered complete unless the analysis of the stability of the stationary solutions is made. For this reason an important analysis of the condensate stability is performed in the section "The analysis of the condensate excitation". In Gaussian approximation the functional, which is an action for different condensate perturbations, is derived in details from the GP functional. We only consider the case of a quasilocal form factor of a special form for this problem (the amplitude of the form factor depends only on the space coordinate r). The solved problem is alternative to the problem of the Bogoliubov transformations used in the studies of the quantum Bose gases in the operator formalism. At the same time, the considered case is an essential generalization of a standard one. The generalizations are the presence of a structure in the system, as well as the arbitrariness of the probe wave function of the condensate perturbations. To obtain the Bogoliubov spectrum, this wave function was chosen in the form of a plane wave.
As the next step, an equation describing the threshold (critical) perturbation momentum was obtained from the Bogoliubov spectrum, at which instabilities arise. This equation is also one of the main results of the paper. The equation has an original form, since it involves suf-ficiently general form of the theory's form factor. At the same time, various dependencies between the characteristics of the condensate and its excitations can be studied with help of this equation. Unlike the simplest problems focused on the Bogoliubov transformation such an equation includes realistic scenarios of the interactions between two subsystems. In particular, it allows us to determine the dependence of the threshold of the condensate excitation from the specific momentum of the condensate. An approximating formula for this dependence is also derived in this paper. The approximating formula has a form of a rapidly converging series and also a universal form in terms of the dimension of the space D (in specific calculations we again only use threedimensional Bose gases, corresponding to the interactions which we derived in the classical theory) and the types of the gases considered. For the specified dependence of the instability threshold from the specific momentum of the condensate the corresponding graphs are constructed for the exponential-trigonometric form factor.
The section "Conclusion" sums up the results, possible further directions in the study of the quantum Bose gases in a model with a form factor, and also a question about the experimental determination of the form factor is discussed. To answer this question the hydrodynamics of the quantum Bose gases of dipoles, quadrupoles and octupoles is briefly considered (in the spirit of the monograph [10]). It follows from hydrodynamic theory that a form factor can be determined experimentally because it directly determines the corresponding generalized Euler equation. Moreover, the obtained system of equations of hydrodynamics is an alternative approach for determination of Bogoliubov spectrum for an excited system which naturally corresponds with the known analysis of the system's stability. We note that an important consequence of hydrodynamics is the fact that quasilocal form factors turn out to be the most "natural" ones: more complex physical mechanisms must be turned on for the emergence of the non-locality of the generalized Euler equation. For this reason the phenomenology of the model proposed in this paper is minimal. Thus, the proposed description of quantum Bose gases in terms of the Gross-Pitaevskii equation with a quasilocal form factor is interesting both from the theoretical point of view and for practical applications such as experiments with Bose gases of dipoles, axisymmetric quadrupoles and octupoles.
Finally, in the "Apendix" section of this paper focused on the classical part of the problem of describing Bose gases detailed derivations of the interaction potentials of classical gases with different values of multipolarity both in the coordinate r and the momentum k representations.

II. MULTIPOLE INTERACTIONS
We begin with the discussion of multipole interactions. In doing so we will consider systems of identical multipoles with axial symmetry which simplifies the calculations. Expressions for multipolar interactions will be the necessary "building blocks" in the problem of the Gross-Pitaevskii equation with a form factor for ultracold Bose gases of a various multipolarity. Detailed calculations related to this section are given in the "Apendix" section.

A. D-D interaction
In the coordinate representation, the expression for the energy of the dipole-dipole interaction is well-known and has the form of: where ψ is the angle between the axis of symmetry of a dipole and the vector r. Next the Fourier transform has to be done, but before that we need to transform from the angle ψ to the angle α between the vector k and the axis of symmetry of the dipole. Here and further we denote the n-th Legendre polynomial by P n . Without loss of generality and to simplify the calculations the axis of symmetry can be considered as laying on the y = 0 plane, then the equation which relates the angles ψ and α can be written in the form: cos ψ = sin α cos ϕ sin θ + cos α cos θ.

(II.2)
In the momentum representation the expression for the dipole-dipole interaction energy is written in the form [24]: where α is the angle between the direction of the dipole moment d and the vector k and A D = 8πd 2 /3. It should be noted that in the momentum representation the D-D interaction does not depend on the absolute value of the vector k. Now let us consider the Q-Q interaction.

B. Q-Q interaction
Calculating the energy of the Q-Q interaction in the r-representation is a simple problem. Consider the interaction between two quadrupoles with equal quadrupole moment Q µµ . Tensor Q µµ has a form of: Here Q = (2z 2 − x 2 − y 2 )/2 dq. The potential and the Performing the calculations taking into account the connection between the angles θ and ψ we obtain the following equation: Or in a more compact way: Then in the momentum representation we have: where α is the angle between the axis of symmetry of the quadrupole and the vector k. After integrating the equation for the energy will be in the form of: The last expression can be written in a more compact and clear way: It is clear that unlike the D-D interaction a dependence from the squared absolute value of the momentum k occurs.

C. O-O interaction
The octupole moment tensor has a form of: For the octupole moment tensor in case of axial symmetry the following equations occur: In addition since O γγα = 0 the following equations hold: (II.13) Here and further O zzz = O. In what follows, we denote O zzz = O. The calculation of the interaction energy in the coordinate representation is a cumbersome but easy task. The interaction potential is given by: (II.14) After calculations the following equation is obtained: For detailed calculations see the section "Appendix". In momentum representation the interaction is: In case of octupoles the potential energy depends on the forth degree of the absolute value of the momentum k.
The generalization for higher-order multipoles is obvious. Graphical illustrations of the interaction energy are shown in Fig. 1 and 2.
D. Comparative analysis

Interaction energy of multipoles in general terms
Using the properties of Legendre polynomials one can easily obtain the following equation: The interaction energy as a function of the angle in rrepresentation and in k-representation is shown in fig.  1 and 2. The tensor of the 2 l -pole moment in the axially symmetric case has a 2l + 1 non-zero component, 2l of which are equal to −M/2 and (2l + 1) are equal to M . Then the "amplitude" of the interaction, which represents the moment tensor contraction with itself, is simply the sum of the squares of the components, i.e. in the case of the 2 l -pole moment, the interaction amplitude is: For the interaction of two axially symmetric multipoles the following expression can be written: where K(l) is a function of the multipole order. Now let us turn to discussion of the results obtained.

Discussion of interactions in momentum representation
The results of calculations of the interaction energy in the k-representation lead to the conclusion that starting from the Q-Q interaction the potential energy is comparable to the kinetic energy of the particles. In the case of the O-O interaction, the potential energy completely suppresses the kinetic energy at large momenta. It should also be noted that the amplitude of the interaction highly decreases with the increase of the multipole order. For example, for d = Q = O, one gets: (II.20) Here we note that these formulas are useful when one wants to investigate the scaling properties of a quantum Bose gas from general considerations, which is important not only for the Gross-Pitaevskii formalism, but also, for example, for the method of the functional renormalization group (with application to the quantum Bose gas). This concludes our discussion of the classical theory of dipolar, quadrupolar and octupolar gases, and we turn to the central section of this paper focused on the Gross-Pitaevskii equation with a form factor (GP equation is widely discussed in the book [10]).

III. THE GROSS-PITAEVSKII EQUATION
The description of the physics of various excitations of quantum Bose gases with a so-called "form factor" is key to this paper. This idea largely repeats the considerations made when one wants to describe a physical phenomenon using the apparatus of stochastic partial differential equations (SPDEs): we offer a semi-phenomenology in which the form factor (in the SPDEs theory this is the so-called "pump function" which is the Fourier transform of the correlator of a random variable) is chosen based on general considerations and is not derived from the microscopic theory.
What is known about such "building blocks" of theory as a form factor or a pump function from the most general considerations? These blocks must correctly reproduce the various distributions of energy, density and a number of other standard physical quantities in the system, give a qualitative picture (whether different wave structures are formed in the system, or there is stochasticity semi-phenomenologically modeled by random force for example). Thus, there are no univocal rules for choosing this building block (the arbitrariness in the choice can be reduced by restrictions due to the renormalization group). In defence of such a vulnerable for criticism state, a number of generally accepted arguments is presented.
In microscopic theories, these building blocks must be generated by various microscopic mechanisms, and their characteristics for a particular problem must be computable. However, if there is no microscopic theory of this kind, within an effective description, which is only a simplified semi-phenomenological version of the (hypothetical) accurate theory, a specific choice of the form factor or the pump function can be justified by general considerations and results.

A. Equations in r and k representations
The first principle in quantum field theory and statistical physics of quantum gases is the statement of the partition function of the theory in terms of the functional integral ( [6][7][8][9]): (III.1) The action of the model consists of the action of the free theory S 0 and the action of interaction S 1 : The interaction S 1 can be of any form, i.e. be a rather complicated function of fields and their derivatives. In this paper we set a simple goal: refusing from the direct calculation of the integral we consider that the result of the integration can be expressed as a saddle-point equation with a form factor which modulates the action S 1 . Thus, we assume the following: For the condensate density n F : where n (τ, r) =ψ (τ, r) ψ (τ, r). The action of a free theory in this case is: In order to obtain the mean field equation one has to find the functional derivativeψ: After calculations the desired Gross-Pitaevskii equation is obtained (here r 12 = r 1 − r 2 ): We are interested in stationary solutions (III.7) thus it is convenient to introduce: . (III.8) In terms of these denotations the equation in real-space representation has the following form: Also, in momentum-space representation: Here and further we work in terms of: The types of form factors and the solutions of the Gross-Pitaevskii equation with a form factor are discussed below (in the spirit of the books [12,13]), in which a detailed analysis of form factors is carried out within the nonlocal quantum field theory).
B. Form factor

The zoology of form factors
In order to develop some intuition regarding the role of form factors and formulate the conditions for choosing the latter, let us turn to specific examples. As is known, the simplest type of a function of two variables is a separable realization. A separable symmetric core F is given by: (III.12) Let us introduce the following notations: Note thatn is a linear functional of the density n. This fact will be important further in the narration. In terms of notation (III.13) the equation (III.9) takes the form: (III.14) Thus, our first conclusion is that for a separable form factor the equation (III.9) takes the form of an ordinary (linear) Schrödinger equation with an effective potential equal to the sum of the initial and form factor additions, the latter itself contains the integral of the desired solution. If a solution of the equation (III.14) is found, one has to use the second equality in (III.13) and, having obtained the matching condition, find the valuen. Here we note that the role of the trap v can be reduced in favor of the form factor f . The next example is the translation-invariant symmetric core F which is given by: (III. 15) In this case, it is convenient to work in the krepresentation (the momentum representation is the most suitable for translation-invariant problems), since in this representation all the expressions will be simple and clear. The expression for this form factor has the form of If l is considered to be given the equation (III.10) becomes algebraic: This equation is simple to solve and the expression for n is determined by the reverse transformation to the rrepresentation. If l is given (in Thomas-Fermi approximation l is equal to zero) the expression for the density is obtained. For further purposes, this material is the hint to answer the question of the final choice of the form factor. Clearly all of that is not a strict proof, but rather a hint to further action. With this knowledge, more complex classes of form factors can be considered, which is done further.

Example of translational invariance violation
Outside the translation-invariant case, many beautiful ways of violation of this invariance can be proposed. We begin by considering the core of the form factor in the krepresentation of the type of (III. 16), but now the Dirac delta function will choose a non-zero value for the total momentum k 1 + k 2 : In this case the left-hand side of the equation (III.10) says: This expression is overloaded for the primary analysis, and the simplest configuration of the parameters has to chosen in (III.18). Let us choose to momenta (N = 2) in a symmetric way: p 1 = p and p 2 = −p, and consider the amplitude function f to be constant, f a (k 1 , k 2 ) = f 0 . The left-hand side of the Gross-Pitaevskii equation in the k-representation will take a form of: The expression (III.20) is a complex functional relation, since it contains the density n taken at three different points. Fortunately, such complexity is illusory. To see this we return to the r-representation with 2f 0 = 1. The form factor F takes the form of: The Gross-Pitaevskii equation (III.9) now contains an integration with respect to only one radius vector r : Considering the function l to be given, this equation can be solved if we divide all by the cosine appearing on the left-hand side of the equation, and then do the Fourier transform. From the obtained algebraic equation the Fourier image of the product of the desired density and the remaining cosine is determined, and then the inverse transformation is done, which gives: An important conceptual moment of our paper follows from the equation (III.23): up to some details, the form factor of type (III.21) is most preferable, since it contains all the necessary features for modeling of roton physics. As it will be shown further, the expression (III.23) remains valid on a qualitative level for a more general class of quasilocal form factors. Here we note that the kdependent form factors do not generate similar answers.
To conclude the subsection, we note the following: in order to avoid the zero in the denominator (inserted there by us), we can select (using momenta p a ) smoother behavior of the integrand, for example: This is the behavior that should be targeted because it follows from the most general considerations. But before proving the last statement we derive the Gross-Pitaevskii equation in the class of quasilocal form factors and solve it.

C. General solution of the GP equation with a quasilocal form factor
Below a general solution of the equation (III.9) for a quasilocal form factor (symmetric in both arguments) will be obtained. The latter is an operator function F from r and ∂ r acting on the Dirac delta function: For this particular form factor the stationary Gross-Pitaevskii equation in r-representation contains an integration with respect to only one radius vector r : (III.26) The equation (III.26) can be rewritten in the form: The obtained equation can be solved with respect to the modulated density n F by doing the Fourier transform.
In the k-representation the equation (III.27) transforms into the algebraic analog the solution of which is found automatically: Now one should go back to the r-representation and then express the density n in terms of n F . As a result n is given by: . (III.29) Finally, the final equation for the condensate density function n is given by: . (III.30) The expression (III.30) has a remarkably simple analytical form, from which we can see all particular cases obtained in the process of obtaining the general solution.
Here we should also mention the important property of the solution (III.30): if the operator function f contains a dependence on r, as in the case of a separable form factor, the trap v does not play a primary role. The obtained general solution (III.30) still contains a large arbitrariness in the form of a function F . To eliminate this arbitrariness, we will now prove some general statements about form factors, and, having determined the final form of the form factor, calculate the integral for the density n.

D. The choice of the form factor
Let us return to the very beginning of our path -the expression for the interaction action S 1 and formulate a more general case than the one considered above. Let the form factor modulate not the density n, but the fields themselvesψ. In this case, the expression for the S 1 action is given by: At the same time N F (τ, r) is given by: Since a description in terms of a functional integral allows for (functional) change of variables under the integral sign, the definition of primary fields is not unique. We turn to the description of our theory in terms of modulated fields (this corresponds to the simplest change of variables). In other words, now let ψ F andψ F be new independent fields (denote them as ϕ andφ). In this case, the Gaussian action is rewritten as follows: The expression (III.33) contains a new Gaussian propagator defined by the following expression: The last expression is extremely important conceptually: it shows that for the reverse form factor the assumption of analyticity is natural (otherwise we get a propagator that does not have the limit of free theory). Now we use the obtained experience in our original problem. Moreover, we will refrain from the kdependence of the form factor entirely. The role of this dependence was demonstrated above, and consisted of running constants in an optical trap. Now we need to use a different degree of freedom, which is the dependence on r. In the light of what has been said, let us consider the following class of form factors, which is: Based on the experience of the equation (III.20), without loss of generality we can assume that the function f depends on one external momentum p in terms of the scalar product pr. This dependence can be organized in several ways but it is technically easier to assume that this dependence is provided by the exponent. It is also convenient to single out the dependence on two exponents which differ in the sign of argument. Finally, the function f must be expanded into a (double) series in terms of its arguments (they should be considered independent when expanding): be rewritten in the form: Then substituting (III.36) into (III.37), we get a simple integration of the infinite sum of the Dirac delta functions, after which we get the desired answer for n. A remarkable property of the latter is that, even for the zero trap, we obtain a natural model of the roton condensate, which is given by: Like the expression (III.38), the equality (III.40) should be considered one of the main results of this paper. As an example, let us consider an exponential form factor. The exponent has good properties for numerical calculations and graphic illustrations. Qualitative conclusions for the density n, made using the example of a concrete form factor, will also be valid in the general case. The exponential form factor itself is given by: f (r) = e ξ(e ipr +e −ipr ) = e 2ξ cos(pr) , f n,m = ξ n+m n!m! , f n = ξ n n! .

(III.41)
Now the graphs of n for various values of the parameters appearing in (III.41) can be constructed. However, before doing this, it is necessary to go back to the denominator which is the most important building block of the expressions (III.38), (III.40) and the exponential density realization. To study this denominator, we turn to the next subsection.

Core and split
We offer in our opinion a more elegant solution, modifying the idea of the core: Using the idea of a "split" in the expressions for the density of a condensate, we get a series of graphical dependencies, which is given below.

E. Condensate density
Here we discuss the behavior of the condensate density as a function of the order of the interaction and the fixed momentum p, using the exponential form factor with ξ = 1. We also study the the rate of convergence of the series in the expression for the density, using the example of graphs. The graphs are shown in Fig. 3-14. Here θ D , θ Q and θ O are the angles at which the corresponding interaction vanishes.
This concludes the discussion of the solution of the Gross-Pitaevskii equation with the form factor in the Thomas-Fermi approximation and we turn to the discussion of quadrupoles in the optical lattice.

F. Quadrupoles in the optical lattice
We consider a stationary condensate thus the transition ψ → ψ exp(−iµ/ ) is valid, and also in the Gross-Pitaevskii equation due to the smallness of the density the last term can be neglected.
The chemical potential µ can vary continuously, but in the course of the analysis it is shown that the formation of condensate in the low-density limit is possible only in the case of discrete values of µ, similar to the spectrum of bound states when considering the singleparticle Schrödinger equation. where k 1 , k 2 , k 3 are the vectors along the x, y, z axis respectively. The justification for the chosen trap is as follows: we consider such a trap, since in the limit of a small k i r → 0 the given potential transforms into the ordinary (anisotropic) harmonic potential, which is: Besides the potential can be represented as a sum of a form: where each term depends only on one coordinate. Due to this, we can search for the wave function in the form of ψ(r) = ψ 1 (x)ψ 2 (y)ψ 3 (z), in other words, use the method of separation of variables. Then, writing out the Laplace operator and using µ = µ 1 + µ 2 + µ 3 , three equations are obtained: These equations can be reduced to the following form: Thus, the Gross-Pitaevskii equation splits into three onedimensional Schrödinger equations in a periodic potential.

Schrödinger equation in a periodic potential
We have the Schrödinger equation in a periodic potential. After the transformations, it can be rewritten in the following form (for certainty, consider the x-coordinate and omit the indices): However, this form is deceptive: the obtained equation is the Mathieu equation, which is very difficult for analytical analysis. Now let us consider the change z → z + π. Since the cosine value will not change, we have two proportional solutions ψ(z) and ψ(z + π). This can be written in the form:T π ψ(z) = f ψ(z + π), f = const. (III.52) If we consider z = 0 and z = π an important expression f 2 = ψ(2π)/ψ(0) = 1 is found. Now the case of small values of the parameter h can be considered.

Weak coupling mode
If the parameter h is small the Mathieu equation is given by:ψ + Eψ = 0.
(III.53) It is obvious that the general solution of this equation is the function ψ = C 1 sin( √ Ez) + C 2 cos( √ Ez). Let us denote ψ + = a cos( √ Ez) and ψ − = b sin( √ Ez). In this case the operatorT π acts as follows: (III.54) Using the solutions for ψ + and ψ − , and also the action of the translation operator on them the Wronskian of these solutions can be constructed, and W [ψ + , ψ − ] = ab √ E. Having written the Wronskian in the explicit form, a quadratic equation for f is obtained. Its solutions are given by: Besides, f + = 1/f − . Saying that f + = exp(iπν), we get that: f + + f − = 2 cos πν = 2ψ + (π) a , (III.56) i.e. ν = √ E. Then we find that: This shows that E = s 2 , where s is an integer. As a result a restriction on µ is obtained: Thus, in the weak-coupling mode (for small values of h) the condensate wave functions match with the wave functions of the particle in the potential well. The last unknown constant is calculated from the normalization of the wave function: We investigated the case of small values of h and now turn to the case of strong coupling mode, i.e. large values of h.

Strong coupling mode
The potential 2h 2 cos(2z) has its minimum in z = ±π/2 and we will use its expansion. Thus it is convenient to rewrite the Mathieu equation as follows: ψ + E + 2h 2 cos (2z ± π) ψ = 0. (III.60) Using the expansion and introducing a new variable the Mathieu equation takes the form of This is the Weber equation, the solution of which is a parabolic cylinder function. But provided that E + 2h 2 2h = 2s + 1, (III.64) parabolic cylinder functions D s become exp(−x 2 /4)He s , where He s is a modified Hermite polynomial. Therefore the energy spectrum has the form of a harmonic oscillator spectrum and µ must satisfy the following condition: The constant is calculated from the normalization of the wave function.

Comparative analysis of weak and strong coupling modes
In the course of the study of the potential v(r) = − i A i sin 2 (k i r) we obtained different spectra of the where the index w denotes the weak coupling mode, and the index s denotes the strong coupling mode. Schematically, the difference between the strong and weak coupling can be represented by the Fig. 15.
This figure shows the different modes of quantizing of the chemical potential from linear to quadratic. It can also be shown that the spectrum of the chemical potential has a fine structure, and this splitting of the levels makes it possible to connect the strong and weak coupling modes. In this case, the splitting of the levels is proportional to exp(−h). This study is given in the work [34], and it is classical, for this reason it is also given in the monograph [35] devoted to the course of quantum mechanics. This completes the section "The Gross-Pitaevskii Equation" and we proceed to the analysis of the condensate excitations.

IV. ANALYSIS OF CONDENSATE EXCITATIONS
In this section we consider condensate excitations in our model and look for their spectrum using the Gross-Pitaevskii functional. Further in the section, we study the dependence of the critical momentum, destroying the stability of solutions of the Gross-Pitaevskii equation, on the condensate momentum. The analysis begins with a revision of the main aspects of the already described material.

A. Intermediate results
The action of the free theory S 0 is given by: The Hamiltonian has a standard form of the sum of kinetic energy and the potential of the trap: The equation for the interaction action S 1 is given by: In the last equation the effective (modulated by the form factor) interaction Γ and the density n are introduced: For the field ψ we choose an ansatz in which it is the sum of a stationary solution of the Gross-Pitaevskii equation in the Thomas-Fermi approximation and a small perturbation δψ: ψ (t, r) = ψ 0 (t, r) + δψ (t, r) , ψ 0 (t, r) = e − iµt ψ 0 (r) . (IV.5) At the same time we look for the perturbation δψ in the following form: We also need the solution of the Gross-Pitaevskii equation itself which in this section is given by: In the next subsection, the Gross-Pitaevskii functional for a specific form factor will be obtained, which will later be used for analysis.
B. Derivation of the GP functional for a form factor of a special form Substituting the above ansatz for ψ and δψ into expressions for the action of the free theory S 0 , the difference ∆S 0 of the action on the perturbed and stationary field configurations: This seemingly cumbersome expression is greatly simplified when taking into account periodic boundary conditions, which give: Such boundary conditions lead to the fact that all linear terms in δψ do not contribute to the difference of actions, from which it follows that they do not contribute to the Gross-Pitaevskii functional also. Substituting the same ansatz into the expression for the interaction action S 1 and again calculating the difference on the perturbed and stationary field configurations in the Gaussian approximation, we get the expression for the difference ∆S 1 : where for brevity χ (t, r) =ψ 0 (t, r) δψ (t, r) + + ψ 0 (t, r) δψ (t, r) . (IV.10) The difference of the actions ∆S 1 is simplified due to periodic boundary conditions (linear in δψ do not contribute to ∆S 1 ), and also due to the Gross-Pitaevskii equation, which gives: When substituting the phase exponents, this term cancels out with the same in the difference of actions ∆S 0 . Thus the Gross-Pitaevskii functional is given by: Implementing the substitutions left for ψ and δψ and also rewriting χ (t, r) = e −iωt C (r) + e iωtC (r) , C (r) =ψ 0 (r) A (r) + ψ 0 (r)B (r) , (IV.12) by integrating over time we get the desired functional: Here we also used the fact that the field ψ 0 is real. It is the Gross-Pitaevskii functional that we will study further in our paper. With its help the equation for the critical momentum, creating an instability in the solution of the corresponding GP equation, will be derived further.

C. Derivation of the equation for the critical momentum
The functional derivative of the perturbation is given by: . (IV.14) In this case the Gross-Pitaevskii functional is given by where for brevity Thus, the problem of analysis of the stability of the Gross-Pitaevskii equation comes down to the problem of finding the Bogoliubov spectrum. The extremum conditions are as follows: from which is clear thatū = u andv = v, and the equations connecting these to variables are From the condition of the existence of a non-trivial solution for this system it follows that the following equation is valid (the determinant vanishes): The obtained equation is the equation for the perturbation spectrum. Now let us derive the equation for the critical momentum, preliminarily transforming the expression for I as follows: Here we use the fact that ψ 0 is given by: The equation for the critical momentum can be written in a compact form: But in the explicit form this equation is rather cumbersome: 2µ .

(IV.22)
This equation determines, for example, the dependence of the absolute value of k 0 on the direction n 0 . Though we study it in another way. Let us consider the limiting cases of this equation, which allow us to discuss the dependence of k 0 on the momentum of the stationary condensate p.

D. Analysis of condensate stability
The equation (IV.22) is quite complicated, so to solve it we use a trick, which is well known from the theory of superconductivity. Let the direction n 0 be so that the main contribution is given by attraction (in the case of the attraction the equation (IV.21) has a solution). In this case, we make the following substitution: (IV.23) Then the equation (IV.22) takes a simple for analysis purposes form: The solution of the obtained equation is given by: Then we use the earlier chosen form factor: f (r) = e ξ(e ipr +e −ipr ) = e 2ξ cos(pr) .
(IV. 26) We emphasize here that the qualitative conclusions for k 0 made using the example of a specific form factor are also valid in the general case. To calculate the integral, we take into account the following equality: To make the transition to the Kronecker delta the following correspondence is used: Now the integral can be calculated: Because of the rapid convergence of the series, we can confine ourselves to a few first terms, which is confirmed numerically. After calculating the integral, the dependence of the critical momentum k 0 on the condensate momentum p should be considered graphically. Before this, it is important to note that in the case of D = 3 for D-D interaction k 0 does not depend on p because the interaction potential U D (k) does not depend on the absolute value of the momentum k.
The above graphs show that the dependence has an asymptote, which was expected, because the Q-Q interaction ∝ k 2 and the O-O interaction ∝ k 4 . In the case of quadrupoles for small p, we obtain that k 0 p. For octupoles, for a small p, the critical momentum is practically constant, but at a certain value of p sharply decreases, reaching a plateau smoothly. Also, the indicated dependence is valid for angles α which are not equal to those angles at which the interaction vanishes (α Q and α O , respectively). Graphs of the dependence of the critical momentum k 0 on the momentum of the stationary condensate p are shown in Fig. 16 and 17. In this paper the properties of classical and quantum Bose gases of dipoles, axisymmetric quadrupoles and octupoles with different multipole order are studied. Classical interaction potentials in corresponding gases in the coordinate and momentum representations were calculated. In this case, various modifications of this interac-tion are discussed in details, which are the classical modification using core and the original modification using "split" -splitting the interaction by a certain feature. In this paper it is the sign of the corresponding interaction range (this representation is very convenient if we want, for example, to investigate the threshold of instability in the considered gas).
Next, the quantum theory of Bose gases of dipoles, axisymmetric quadrupoles and octupoles in the Gross-Pitaevskii formalism is discussed. The starting point of this consideration is the Gross-Pitaevskii functional, from which the same name equation is derived in both coordinate and momentum representations. The original point of this paper is that we studied not the "ordinary phonon" condensate and its excitations, but considered a more general scenario for the appearance of a spatially inhomogeneous condensate and its excitations (of an arbitrary nature). To achieve this, we introduced the socalled "form factor". This semi-phenomenological model of the Gross-Pitaevskii equation with a form factor allows to describe the physics of practically any quasiparticles arising in the considered quantum Bose gases.
The zoology of the form factors appearing in the GP equation is studied in detail. It is concluded that for a wide class of physical scenarios, the so-called quasilocal form factors are sufficient. We note that in the general case this form factor leads to a violation of translational invariance in the system, which is physically transparent: the distribution of the condensate is spatially nonuniform. Further, in the Thomas-Fermi approximation, a general solution of the GP equation with a quasilocal form factor is obtained. This solution has an interesting form in terms of a double rapidly convergent series, which can easily be shown by direct construction of the majorant of the corresponding series. Let us note some more beautiful properties of the obtained solution. First of all, it has a universal form with respect to the dimension of the space D (for specific calculations we limited ourselves to three-dimensional Bose gases, the corresponding interactions in which were obtained in the classical theory). Also, the solution obtained universally includes all considered types of Bose gases.
To illustrate the results obtained, graphs of condensate density functions for the exponential-trigonometric form factor are constructed. The exponent has good properties for numerical calculations and graphical constructions, and the periodic argument simulates the presence of a condensate density wave in the system. For the sake of completeness the GP equation with the optical lattice potential in the limit of small condensate concentrations is also considered in this paper. This limit does not distinguish between dipolar, quadrupolar and octupolar gases, and the equation itself is the well-known Schrödinger equation in the periodic potential (its stationary case is the Mathieu equation). The paper gives a brief discussion of the latter.
Then an important analysis of the stability of condensate was performed, in other words, a study of condensate excitations. In the Gaussian approximation, a functional describing the perturbations of the condensate is derived in detail from the Gross-Pitaevskii functional. We only consider the case of a quasilocal form factor of a special type. We note that this problem is a generalized analog of the Bogoliubov transformation used in the study of quantum Bose gases in operator formalism. In addition to the presence of a structure in the system, another generalization is that the probe wave function of the condensate perturbation does not have to be a plane wave, which, however, was chosen in this paper in order to obtain the spectrum of Bogoliubov excitations.
From the Bogoliubov spectrum, an equation describing the threshold perturbation momentum for the onset of the instability is obtained. This equation has an original form, since it includes the form factor of the theory in a complicated way. Another important result of our paper is that this equation makes it possible to establish the dependence of the threshold on the parameters of a stationary condensate. The latter is demonstrated by the example of the dependence of the threshold on the characteristic momentum of the condensate. For the sake of completeness, an approximating expression for the corresponding dependence is obtained in the paper. The approximating equation has the form of a certain rapidly converging series. The last statement is again easily proved by direct construction of the majorant. An interesting property of the equation obtained is that it has a universal form with respect to the dimension of the space D (for specific calculations we again limited ourselves to three-dimensional Bose gases, the corresponding interactions in which were obtained in the classical theory). The graphs of the corresponding series for the exponential-trigonometric form factor are constructed.
In conclusion of the paper, let us note the question of the experimental determination of the form factor of the theory. What approach can give an appropriate hint? To answer this question, let us consider the hydrodynamics of quantum Bose gases of dipoles, quadrupoles, and octupoles. To this end, we represent the wave function of a Bose gas in the form: The hydrodynamic velocity is given by: Then the continuity equation for condensate density is derived from the Gross-Pitaevskii equation: This equation is one of the two equations describing the hydrodynamic of ultracold gases.The second equation is the generalized Euler equation, which is given by: n (t, r) , Let us consider the last term of the Euler equation and do the substitution: Then this term is given by: So we have obtained the equalities that answer the question. It follows from these equalities that the form factor can be determined experimentally by conducting experiments on the hydrodynamics of ultracold gases, since they explicitly determine the corresponding generalized Euler equation. Moreover, these equations are an alternative approach for determining, in particular, the Bogoliubov spectrum for system excitations. The latter is obtained by a well-known analysis of the stability of the system of equations of hydrodynamics. Hydrodynamic equalities also show that quasilocal form factors are the most "natural": for the appearance of non-locality of the generalized Euler equation, more complex physical mechanisms may be required. For this reason, the phenomenological nature of the model is minimal. Moreover, even this model is not fully described in the paper: some calculations are performed for a quasilocal form factor of a special type. Calculations in the general case can be the subject of a separate publication. Also, a deeper study of the hydrodynamics and thermodynamics of quantum Bose gases in a model with a form factor deserves special attention, in particular, the consideration of temperature. Another interesting problem arises here: the modification of the form factor in the case of finite temperatures. Thus, the proposed description of quantum Bose gases in terms of the Gross-Pitaevskii equation with a form factor is interesting both from the theoretical point of view and for practical applications such as experiments with Bose gases of dipoles, axisymmetric quadrupoles and octupoles.

A. Calculation of the k-representation of the D-D interaction
In this appendix, detailed derivations of the interaction potentials of classical gases with different multipole values both in the coordinate r and in the momentum k representations are given. We begin with calculating the D-D interaction in the k-representation.
The expression for the D-D interaction energy in rrepresentation is given by: In a general case the angle ψ between the vector of the dipole moment d and the vector r is related to the angle between the vector d and k as follows: cos ψ = sin α sin θ cos(ϕ − γ) + cos α cos θ. (VI.2) Indeed, it follows directly from the Fig. 18. The projection on some unit vector a is given by: ar r = sin α cos γ sin θ cos ϕ+ + sin α sin γ sin θ sin ϕ + cos α cos θ. (VI.3) From this expression after the transformations (taking into account the spherical coordinate system) the connection between ψ and α is obtained. From this formula it is clear that one may put γ = 0, because in the course of calculation of the k-representation, one needs to integrate the expression with cos(ϕ − γ) over the period 2π. The integral that has to be calculated is given by: (VI.4) After the integration over ϕ and the substitution x = cos θ with further integration over x in corresponding limits, we see that the integral has a divergence at zero when integrating with respect to r. Thus we integrate beginning from a certain r 0 . This leads to the expression: − 2π(1 + 3 cos 2α)(r 0 k cos r 0 k − sin r 0 k) r 3 0 k 3 , (VI.5) which in the limit of r 0 → 0 gives U D (k, α) = 4πd 2 3 (3 cos 2 α − 1). (VI. 6) It is clear that the expression does not depend on k.

B. Calculation of the k-representation of the Q-Q interaction
In a general case the angle ψ between the quadrupole axis and the vector r is related to the angle α between the vector d and k by the equation (VI.2). α is the angle between the quadrupole axis and the vector k.
The integral that has to be calculated is given by: This is the expression for the Q-Q interaction energy in the k-representation.

C. Calculation of the r-representation of the O-O interaction
The calculation of the interaction has no conceptual complexity, but it is very cumbersome, even if taking into account the symmetry of the considered case. We will describe the general idea of calculations, and the details can be implemented by using any of computer algebra systems and programming languages.
Calculations begin with an expression for the potential created by the quadrupole: ϕ(r) = O λνµ r λ r ν r µ r 4 . (VI.11) Next, let us remind that the expression for the energy of the O-O interaction in the r-representation is given by: (VI.12) Thus, the following equation needs to be calculated first: Then, we should use the fact that the octupole moment tensor has only 7 independent components and symmetry properties which were mentioned earlier. This lets us separate the resulting terms into groups, and then individually simplify each group. Finally, it is necessary to substitute the spherical coordinates into the resulting expression and simplify it.

D. Calculation of the k-representation of the O-O interaction
To calculate the O-O interaction in momentum space, one must substitute the expression for the angles ψ and α. After this, the problem simplifies to calculating the integral: Integration over φ is obvious, integration over θ is simplified by the substitution x = cos θ. When integrating over r, there are no problems with divergence. The whole process is easy to do in any suitable system of computer algebra, there is no point in writing the intermediate calculations because of their cumbersomeness.
The final result is given by: After simplifications one gets the expression given in the corresponding section.