Exact solutions for the dispersion relation of Bogoliubov modes localized near a topological defect - a hard wall - in Bose-Einstein condensate

We consider a Bose-Einstein condensate of bosons with repulsion, described by the Gross-Pitaevskii equation and restricted by an impenetrable"hard wall"(either rigid or flexible) which is intended to suppress the"snake instability"inherent for dark solitons. We solve analytically the Bogoliubov - de Gennes equations to find the spectra of gapless Bogoliubov excitations localized near the"domain wall"and therefore split from the bulk excitation spectrum of the Bose-Einstein condensate. The"domain wall"may model either the surface of liquid helium or of a strongly trapped Bose-Einstein condensate. The dispersion relations for the surface excitations are found for all wavenumbers $k$ along the surface up to the"free-particle"behavior $k \rightarrow \infty$, the latter was shown to be bound to the"hard wall"with some"universal"energy $\Delta$.


I. INTRODUCTION
Initially the Gross-Pitaevskii equation (GPE) was intended to describe structures and excitations in superfluid helium [1], [2]. Being a non-linear Schrödinger equation it possesses a broad spectrum of applications for various non-linear processes in condensed matter such as bright and dark solitons in Bose-Einstein condensate (BEC) and non-linear optics [2], as well as the waves of finite amplitude on the surface of liquid [3]. BEC disturbances known as Bogoliubov excitations [4], [5] are described by the eigenmodes of the matrix Bogoliubov-de Gennes equation (BdGE) that follows from GPE and has become an archetype in many fields ranging from superconductivity [6] to the gravitational black hole analogy in BEC [7], [8]. The ubiquitous nature of GPE and BdGE demands rigorous analytical solutions though not many of them have been obtained.
Domain wall solutions of GPE such as 2D dark solitons are known to be unstable except the case of a solid wall. Maybe this explains the fact that the BdG equation has not been paid much attention to find the localized solutions near such walls. Yet even the case of the solid wall deserves investigation as far as it is connected with the generic topic of edge excitations in topological phases. The situation may bear some resemblance to two-band models with Majorana bound states that arise as solutions to three dimensional BdG theories. The gapless modes that propagate along a physical boundary, while they are exponentially localized away from the physical boundary are gapless boundary modes or edge states.
The surface excitations in restricted BECs (superfluid helium 4 with BEC confined in pores [9], self-bound BEC at the surface of superfluid helium [10], as well as at the surface of BEC trapped in an external potential [11], or * Electronic address: peter@snu.ac.kr BEC near a solid wall [12]) are of fundamental interest. One way to obtain surface excitations of BEC was considered in the work by Anglin et al. [11] which treats the surface excitation of stable BEC (half analytically, half numerically) in the presence of an external linear trapping potential.
Unlike the "soft" trapping potential [11], in order to make the problem analytically tractable one may consider the extreme boundary condition of the solid wall for the surface of trapped BEC which stability was proven in [12] by showing that the imposition of the boundary condition of zero wavefunction on the wall ensures the stability of the solution near a solid wall in spite of the fact that the "domain wall" itself is essentially unstable. In other words, the "hard wall" condition means that the wavefunction of BEC vanishes at the place where "a steep repulsion with a turn-on length smaller than the healing lengths and penetration depths of the condensates" exists [13]. For example, potentials steeper than harmonic were prepared by using Laguerre-Gauss doughnut-shaped laser beams for a BEC container [14]. An inhomogeneous stationary solution of GPE (the "domain wall") which coincides with the half of the dark soliton at rest (the kink soliton ψ 0 = tanh(x), the distance x is measured in the units of the healing length, see below) [2] may have one of its physical realization as a model for BEC near a solid wall [12] where localized Bogoliubov excitations were proposed to exist [15]. However, an exact analytical solution for correspondent surface-bound excitations (if any) has not been found.
Although the stationary dark soliton in BEC was proven to be unstable with the "snake instability" [12], [16], the hard wall boundary condition may also approximate the sharpness of a self-bound potential at a free surface of liquid helium which was proven to be composed of nearly 100% BEC that satisfies GPE [10]. One may consider the free kink wall of ψ 0 = tanh(x) as a model for a free surface demanding only the topological stability of such a solution in which its nodal surface undergoes weak flexural oscillations. For this case the position of the "hard wall" is flexible (like an impenetrable membrane on the surface of helium II) and it imitates the free surface of the liquid. Then the role of a hard wall container is played by the liquid surface of helium II.
Here we consider the problem of the localized gapless excitation modes by finding analytical solutions of a matrix Schrödinger equation being a slight (but important for what follows) modification of BdGE [5] as such as given in [12], [16], and [15]. The binding energy of localized excitations is of our concern in the present work. Surprisingly, we found that the spectrum of surface excitations can be calculated analytically for any k to be compared with the numerical results and with the analytical results obtained for limiting cases k → 0 and k → ∞. The natural limit of k → ∞ which in the bulk BEC results in the energy spectrum ε = ( k) 2 /2m + µ where m is the mass of the boson and µ = gn 0 is the chemical potential while g and n 0 are the intensity of the repulsion and the BEC particle density, correspondingly, leads to ε = ( k) 2 /2m + µ − ∆. The analytical approach developed below for solving the Bogoliubov-de Gennes equations may turn useful for other applications.

II. BASIC EQUATIONS
GPE can be written as [1]: We introduce dimensionless quantities by measuring distances in the units of the healing length ξ = /mc and the energy in the units of gn 0 = mc 2 where c = gn 0 /m is the sound velocity. The stationary equation (1) for the kink with the node at the position x = 0 gives ψ 0 = tanh(x) of the soliton . We will disturb this solution to investigate its Bogoliubov excitations by presenting ψ of Eq. (1) as a sum of plane waves [17]: , ̺ lies in the plane orthogonal to x direction (we consider x ≥ 0 and all the functions decaying exponentially in this direction), k is the wave vector along this plane and * denotes complex conjugation. We will suppress the indices and simplify the notation by using a and b instead of a ω, k (x) and b ω, k (x). Introducing the functions ψ 1 = a + b and ψ 2 = a − b, after linearizing Eq. (1) we get the pair of coupled Schrödinger equations [15]: where κ = ξ| k|/ √ 2 and ε = ω/(mc 2 ). This pair of equations is identical to the corresponding Bogoliubovde Gennes equations (see [4] and [5]) if one rewrites them for the functions a and b, and also sets κ = 0 and ε = 0. As far as we know, Eqs. (2),(3) have never been solved before for arbitrary non-zero κ and ε.
We find a formal general solution for these equations and illustrate its viability by obtaining the rigorous solution of the spectrum of localized phonons. The spectrum of bulk excitations can be easily found from (2) and (3) when neglecting the derivative terms far from the boundary x = 0 to obtain the well-known Bogoliubov spectrum and for κ → ∞ ε B ≈ κ 2 + 1 which is a free boson plus chemical potential. Any localized excitations should have the energy spectrum lying lower than the bulk one.

III. SUPERSYMMETRY OF BDGE
It is interesting to note that Eqs. (2), (3) with ε = 0 and κ = 0 are parts of a supersymmetric Hamiltonian with zero ground state energy. Indeed, on introducing the matrix operator so that the l.h.s. of Eqs.

) with its partner Hamiltonian
to produce a supersymmetric (SUSY) Hamiltonian that may canonically be expressed through the super-chargesQ as anticommutator The supersymmetry is explicitly broken at ε > 0 and κ > 0 which eventually leads to splitting the degenerate zero ground state into two gapless excitations (a "light" one with ε ∝ κ and a "heavy" one with ε ∝ κ 3/2 ) [15], both bound to the wall.

IV. BOUNDARY CONDITIONS
The boundary conditions for ψ 1 and ψ 2 can be of two kinds [15]. At the node of the kink ψ = 0 that is both Im ψ = 0 and Re ψ = 0 therefore ψ 2 = 0 and ψ 1 = 0. However, for ψ 1 the additional possibility exists. Indeed, for κ = 0 and ε = 0 Eqs. (2),(3) have the solutions ψ 0 1 = 1 − ψ 2 0 and ψ 0 2 = ψ 0 the first of which is the so-called "zero mode" [15], [4] or Goldstone gapless mode corresponding to a translation of the kink ψ 0 as a whole along x, being the derivative of the kink ψ 0 : ψ 0 (x + δx) ≈ ψ 0 (x) + ψ 0 1 δx. Thus the condition Re ψ = 0 turns into ψ ′ 0 δx( ̺, t) + Re ϑ( r, t) = 0 which determines the shape of the loci of nodes δx( ̺, t) (the shape of the surface). The derivative of such a mode with respect to x is zero at x = 0. Thus the mode with the mixed boundary conditions d dx ψ 1 | x=0 = 0 and ψ 2 | x=0 = 0 would allow the rippling of the soliton and will be called the "ripplon" mode (predicted in [15]). As we shall see below its energy spectrum at low κ coincides with the one for the capillary wave. The mode with the zero boundary conditions ψ 1 | x=0 = 0 and ψ 2 | x=0 = 0 which correspond to a flat boundary of the solid wall will be called the "surface phonon" mode (as far as its spectrum starts linear) and was predicted in [15]. Finally, the very condition of the solid wall excludes the possible solution xψ 0 − 1 [16] of Eq. (3) at κ = 0, ε = 0 which could be responsible for the "snake" instability of the "domain wall" and which does not satisfy the zero boundary conditions.

V. LONG-WAVELENGTH APPROXIMATION SOLUTION
First consider the case of κ → 0. In case of the ripplon spectrum we find ψ 1,2 as a series in ε: ψ 1 ≈ ψ 0 1 + εψ 1 1 + O(ε 2 ) and ψ 2 ≈ ψ 0 2 + εψ 1 2 + O(ε 2 ). A zero approximation is the solution of homogeneous equations Eqs. (2), (3) with ε = 0. The solutions can be found for any κ (this can be verified by the direct substitution): where α 1 = √ 2 √ 2 + κ 2 and α 2 = √ 2κ. To find ψ 1 1 and ψ 1 2 we have to solve the inhomogeneous equations that follow from Eqs. (2), (3) where κ = 0: With the help of the Green functions of the homogeneous equations the inhomogeneous solutions are found as: Finally, the derivative with respect to x of ψ 1 at x = 0 is found from Eqs. (10), (14) to be ψ ′ 1 = Aα 1 (2 − α 1 )(2 + α 1 )/3 + εB, which according to the mixed boundary conditions should be zero together with ψ 2 = −Aε + Bα 2 , as it follows from Eqs. (11), (15). The zero determinant with respect to A and B det gives the ripplon spectrum taking into account that α 1 ≈ 2 + κ 2 /2 for κ → 0 and retaining only the lowest power of κ ε = 4 Spectrum (17) is shown in Fig. 1. Note that the localization of the ripplon at low κ is governed by α 2 = √ 2κ which is compared below with the numerical solution for α 2 . As it was already noted in [15] the spectrum (17) exactly coincides with the well-known expression for the frequency of the capillary waves multiplied by the Planck constant when written in a dimensional form: ε = σ/mn 0 k 3/2 where σ = 2/3 cn 0 is the surface energy of the stationary soliton ψ 0 [1]. In fact, σ is exactly the half of the energy of the dark soliton at rest (see Eq. (5.59) in [2]).
The zero boundary conditions lead to surface phonons for κ → 0 [15] and below we obtain the whole spectrum analytically. Here we only mention that α 2 for phonons at low κ is proportional to κ 2 which indicates much weaker localization, in contrast to the ripplons.

VII. EXACT SOLUTION OF BDGE
Let us now find the exact solution of Eqs. (2), (3) at arbitrary κ. To do this let us transform these equations into a single matrix hypergeometric equation. Matrix generalizations of both hypergeometric function and gamma function were shown to be mathematically correct (see Refs. ( [19]) and ( [20])). Introducing z = (1 − ψ 0 )/2 and .
After introducing the matrices, Eq. (22) becomes where primes mean differentiation with respect to z. Equation (30) has a formal solution as the matrixvalued hypergeometric function [19] where (â,b,ĉ) 0 = 1; We use this formal solution to obtain the spectrum of the phonon localized near the soliton. The boundary condition at x = 0 (that is at z = 1/2) will be fulfilled whenΦ = 0. As far as the hypergeometric function at z = 1/2 can be expressed through the matrix gamma function ( [20]) (becauseĉ = (1 +â +b)/2, see Eq.(26)) so that then the condition of zero Φ means that the matrix (33) has a zero eigenvalue, that is the determinant of (33) should be zero. In turn, matrix gamma functions could be presented as the products of matrices [20] Γ To obtain the analytical results we will instead utilize the productP = ((1 +â)(1 +b))/2 = (1 +â +b +âb)/2 that uses already calculated matrices. It is easy to obtainP from (26),(27) and (28) Finally, the spectrum of the surface phonons is determined by the equation One can make sure that Eq. (36) reproduces the spectrum calculated before at κ → 0 and κ → ∞. Indeed, at κ → 0 the spectrum is ε = √ 2κ + O(κ 5 ) so that the κ 3 term is missing as it was predicted in [15], while the bulk phonon starts with higher energy as ε B = √ 2κ + κ 3 /(2 √ 2) + O(κ 5 ). Let us define the absolute value of the binding energy as ∆ε = ε B − ε. Then it starts as κ 3 /(2 √ 2) (see Fig. 2(c)). Now let us consider the other limit κ → ∞. It is easy to see that seeking the solution in the form ε ≍ κ 2 + 1 − ∆ leads to r = p ≍ κ/2 + √ 2∆/4 and l ≍ κ/2 − √ 2∆/4 which after substitution into Eq. (36) give 2∆ + 3 √ 2∆ − 2 = 0. It has the same root as we found before from Eq. (21) √ 2∆ = ( √ 17 − 3)/2 = α ∞ ≈ 0.562 and therefore ∆ε ∞ = ∆ = α 2 ∞ /2 ≈ 0.158. Such a coincidence with exact asymptotic results found before brings confidence to Eq. (36) which is enhanced by a good correspondence with the results obtained by direct numerical solution of the coupled Shrödinger equations (Fig. 2). Note that the slower decay exponent α 2 can be approximated by a simple expression α 2 = κ 2 /(1 + κ 2 /α ∞ ) that fits the exact expression of Eq. (29) with ε from the exact solution of Eq. (36) within 0.2%. We plot the decay exponential in Fig. 2(d) together with the results of the numerical solution. Unfortunately, the case of mixed boundary conditions for ripplons can not be treated in the same way.

VIII. CONCLUSION
In conclusion, we have found the analytic spectra and wavefunctions of localized Bogoliubov elementary excitations existing near the inhomogeneous stationary solution of the Gross-Pitaevskii equation which may represent a physical model for the surface of liquid helium. In this status our solutions predict surface modes -ripplons and surface phonons -that may contribute to the thermodynamics of the surface at low temperatures [15]. We believe that our results and the method for exact solving the Bogoliubov-de Gennes equations could be useful in more sophisticated cases involving solitons.