Exact surface-wave spectrum of a dilute quantum liquid

TL;DR

This paper derives the exact surface-wave spectrum of a dilute Bose-Einstein condensate near a boundary, revealing precise dispersion relations for all wavenumbers and connecting theoretical predictions with experimental observations.

Contribution

It provides the first exact analytical solutions for the surface excitation spectrum of a dilute quantum liquid for all wavevectors.

Findings

Exact dispersion relations for surface excitations at all wavenumbers.

Identification of a maximal binding energy for surface modes.

Close agreement of theoretical binding energy with experimental data.

Abstract

We consider a dilute gas of bosons with repulsive contact interactions, described on the mean-field level by the Gross-Pitaevskii equation, and bounded by an impenetrable "hard" wall (either rigid or flexible). We solve the Bogoliubov-de Gennes equations for excitations on top of the Bose-Einstein condensate analytically, by using matrix-valued hypergeometric functions. This leads to the exact spectrum of gapless Bogoliubov excitations localized near the boundary. The dispersion relation for the surface excitations represents for small wavenumbers $k$ a ripplon mode with fractional power law dispersion for a flexible wall, and a phonon mode (linear dispersion) for a rigid wall. For both types of excitation we provide, for the first time, the exact dispersion relations of the dilute quantum liquid for all $k$ along the surface, extending to $k \rightarrow \infty$. The small wavelength excitations are shown to be bound to the surface with a maximal binding energy $\Delta= \frac18 (\sqrt{17}-3)^2 mc^2 \simeq 0.158\, mc^2$, which both types of excitation asymptotically approach, where $m$ is mass of bosons and $c$ bulk speed of sound. We demonstrate that this binding energy is close to the experimental value obtained for surface excitations of helium II confined in nanopores, reported in Phys. Rev. B 88, 014521 (2013).

Equations

$$i\hslash \frac{\partial \psi }{\partial t}=-\frac{{\hslash }^2}{2m}{\nabla }^2\psi +g{n}_0(|\psi {|}^2-1)\psi .$$

$$-\frac{1}{2}\frac{d^2}{d{x}^2}{\psi }_1+(3{\psi }_0^2-1+{\kappa }^2){\psi }_1$$

$$-\frac{1}{2}\frac{d^2}{d{x}^2}{\psi }_2+({\psi }_0^2-1+{\kappa }^2){\psi }_2$$

$$\displaystyle\hat{A}=\left(\begin{array}[]{cc}-\frac{1}{\sqrt{2}}\frac{d}{dx}-\sqrt{2}\psi_{0}&0\\ 0&-\frac{1}{\sqrt{2}}\frac{d}{dx}+\frac{1}{\sqrt{2}}\frac{1-\psi_{0}^{2}}{\psi_{0}}\end{array}\right)$$

$$\displaystyle\hat{H}_{-}=\hat{A}^{{\dagger}}\hat{A}=\left(\begin{array}[]{cc}-\frac{1}{2}\frac{d^{2}}{dx^{2}}+3\psi_{0}^{2}-1&0\\ 0&-\frac{1}{2}\frac{d^{2}}{dx^{2}}+\psi_{0}^{2}-1\end{array}\right)$$

$$\displaystyle\hat{H}_{+}=\hat{A}\hat{A}^{{\dagger}}=\left(\begin{array}[]{cc}-\frac{1}{2}\frac{d^{2}}{dx^{2}}+\psi_{0}^{2}+1&0\\ 0&-\frac{1}{2}\frac{d^{2}}{dx^{2}}+\frac{1-\psi_{0}^{2}}{\psi_{0}^{2}}\end{array}\right),$$

$$\displaystyle\hat{H}_{\rm SUSY}=\left(\begin{array}[]{cc}\hat{H}_{-}&0\\ 0&\hat{H}_{+}\end{array}\right)$$

$$\displaystyle\hat{Q}=\left(\begin{array}[]{cc}0&0\\ \hat{A}&0\end{array}\right),\,\hat{Q}^{{\dagger}}=\left(\begin{array}[]{cc}0&\hat{A}^{\dagger}\\ 0&0\end{array}\right)$$

$${\hat{H}}_{\mathrm{SUSY}}=\{Q,Q^{†}\};{\hat{Q}}^2=0,({\hat{Q}}^{†}{)}^2=0.$$

$${\psi }_1^{(0)}=A\exp (-{\alpha }_{1}x)\left(\frac{{\alpha }_1^2-1}{3}+{\alpha }_1{\psi }_0+{\psi }_0^2\right),$$

$${\psi }_2^{(0)}=B\exp (-{\alpha }_{2}x)({\psi }_0+{\alpha }_2),$$

$$-\frac{1}{2}\frac{d^2}{d{x}^2}{\psi }_1^{(1)}+\left(3{\psi }_0^2-1\right){\psi }_1^{(1)}$$

$$-\frac{1}{2}\frac{d^2}{d{x}^2}{\psi }_2^{(1)}+\left({\psi }_0^2-1\right){\psi }_2^{(1)}$$

$${\psi }_1^{(1)}=\frac{1}{2}B\left\{{\psi }_0+x(1-{\psi }_0^2)\right\},$$

$${\psi }_2^{(1)}=-A.$$

$$\displaystyle\det\left(\begin{array}[]{cc}\alpha_{1}(2-\alpha_{1})(2+\alpha_{1})/3&\varepsilon\\ -\varepsilon&\alpha_{2}\end{array}\right)=0$$

$$\epsilon =\sqrt{\frac{4\sqrt{2}}{3}}{\kappa }^{3/2}.$$

$$-\frac{1}{2}\frac{d^2}{d{x}^2}{\psi }_{\infty }+(3{\psi }_0^2-2+\Delta ){\psi }_{\infty }$$

$$-\frac{1}{2}\frac{d^2}{d{x}^2}{\psi }_{\infty }+({\psi }_0^2-2+\Delta ){\psi }_{\infty }$$

$$\begin{gathered} {\psi }_{\infty }=(1-{\psi }_0^2{)}^{{\alpha }_{\infty }/2} \\ \times {}_2{F}_1\left({\alpha }_{\infty }-s,{\alpha }_{\infty }+s+1,{\alpha }_{\infty }+1,\frac{1-{\psi }_0}{2}\right), \end{gathered}$$

$$\begin{gathered} {}_2{F}_1\left({\alpha }_{\infty }-s,{\alpha }_{\infty }+s+1,{\alpha }_{\infty }+1,\frac{1}{2}\right)= \\ \frac{\Gamma (\frac{1}{2})\Gamma ({\alpha }_{\infty }+1)}{\Gamma (\frac{1}{2}[1+{\alpha }_{\infty }-s])\Gamma (\frac{1}{2}[2+{\alpha }_{\infty }+s])}=0. \end{gathered}$$

$${\psi }_{1,2}=(1-{\psi }_0^2{)}^{\alpha /2}{\varphi }_{1,2},qquadz=\frac{1-{\psi }_0}{2},$$

$$\begin{gathered} z(1-z)\frac{d^2{\varphi }_1}{d{z}^2}+\left\{\alpha +1-2(\alpha +1)z\right\}\frac{d{\varphi }_1}{dz} \\ +\left\{6-\alpha (\alpha +1)\right\}{\varphi }_1+\frac{1}{2z(1-z)}\left(\frac{{\alpha }^2}{2}-2-{\kappa }^2\right){\varphi }_1 \\ =\epsilon \frac{{\varphi }_2}{2z(1-z)}, \\ z(1-z)\frac{d^2{\varphi }_2}{d{z}^2}+\left\{\alpha +1-2(\alpha +1)z\right\}\frac{d{\varphi }_2}{dz} \\ +\left\{2-\alpha (\alpha +1)\right\}{\varphi }_2+\frac{1}{2z(1-z)}\left(\frac{{\alpha }^2}{2}-{\kappa }^2\right){\varphi }_2 \\ =\epsilon \frac{{\varphi }_1}{2z(1-z)}. \end{gathered}$$

$$\displaystyle\hat{\Phi}={\left(\begin{array}[]{c}\phi_{1}\\ \phi_{2}\end{array}\right)},$$

$$\displaystyle\hat{\alpha}^{2}=2\left(\begin{array}[]{cc}2+{\kappa^{2}}&{\varepsilon}\\ {\varepsilon}&{\kappa^{2}}\end{array}\right),$$

$$\hat{c}=\hat{\alpha }+\hat{1},$$

$$\hat{1}+\hat{a}+\hat{b}=2(\hat{\alpha }+\hat{1}),$$

$$\displaystyle-\hat{a}\hat{b}=\left(\begin{array}[]{cc}6&0\\ 0&2\end{array}\right)-\hat{\alpha}^{2}-\hat{\alpha}.$$

$$\displaystyle\hat{\alpha}={2}\left(\begin{array}[]{cc}r&l\\ l&p\end{array}\right),$$

$${\alpha }_{1,2}=\sqrt{2}\sqrt{1+{\kappa }^2\pm \sqrt{{\epsilon }^2+1}}$$

$$z(1-z){\hat{\Phi }}^{\prime \prime }+(\hat{c}-(\hat{1}+\hat{a}+\hat{b})z){\hat{\Phi }}^{\prime }-\hat{a}\hat{b}\hat{\Phi }=0,$$

$$\det \hat{P}=(3r+{\kappa }^2)(3p+{\kappa }^2)-(3l+\epsilon {)}^2=0,$$

$$\begin{gathered} {\left(\hat{1}+\frac{\widehat{U_1}}{2}\right)}_{0,0}{\left(\hat{1}+\frac{\widehat{U_1}}{2}\right)}_{1,1}-{\left(\hat{1}+\frac{\widehat{U_1}}{2}\right)}_{0,1}{\left(\hat{1}+\frac{\widehat{U_1}}{2}\right)}_{1,0} \\ =0 \end{gathered}$$

$$\displaystyle{\left(\begin{array}[]{c}\phi_{1}\\ \phi_{2}\end{array}\right)}={}_{2}F_{1}(a,b,c,z)\,{\left(\begin{array}[]{c}A\\ B\end{array}\right)},$$

$${\varphi }_1$$

Full text

I Introduction

Initially, the Gross-Pitaevskiǐ equation (GPE) was intended as a model to describe structures and excitations in superfluid helium.Gross (1961); Pitaevskii (1961) Being a nonlinear Schrödinger equation, it was however recognized later on that it possesses a variety of applications for various nonlinear processes in condensed matter such as bright and dark solitons in dilute Bose-Einstein condensates (BECs, for which the GPE is accurate on the mean-field level)Pitaevskii and Stringari (2003) and nonlinear optics,Carusotto and Ciuti (2013), as well as finite amplitude waves on the surface of a liquid.Zakharov (1968) Excitations on top of the mean-field ground state representing the BEC, known as Bogoliubov excitations,Bogoliubov (1947) are described by the eigenmodes of the matrix Bogoliubov-de Gennes equations (BdGE). The associated quanta of the perturbation field have become the archetype of quasiparticle excitations in superconductivity Bogoliubov (1958); Valatin (1958); de Gennes and Saint-James (1963) and the theory of dilute quantum gases,Leggett (2001), inter alia also for the formulation of the propagation of quantum fields on effective curved spacetimes.Kurita et al. (2009) The ubiquitous nature of the BdGE makes rigorous analytical solutions highly desirable, but very few, and only in limiting cases, have been obtained.

Domain wall solutions of the GPE such as 2D dark solitons are known to be unstable except for those in the presence of a hard wall. However the case of a hard wall deserves investigation in particular because it is connected with the generic topic of edge excitations in topological phases. Specifically, the corresponding physical situation bears some resemblance to two-band models with Majorana bound states that arise as solutions to a BdG approach. The gapless modes that propagate along a physical boundary, while they are exponentially decaying away from the physical boundary, are gapless boundary modes or edge states.Hasan and Kane (2010)

Examples for the occurrence of surface excitations in bounded BECs comprise, for example, superfluid ${}^{4}\!$He (helium II) confined in pores,Shams et al. (2006) self-bound condensates at the low-density surface of superfluid helium Griffin and Stringari (1996), as well as surface states of a BEC trapped in an external potential Anglin (2001), or surface states of other media with a defocusing nonlinearity.Kuznetsov and Turitsyn (1988) They are of fundamental interest since they reveal the role of quantum effects on the excitation character (i.e., effects which are not existing on the classical level) in restricted geometries.

Considering the boundary condition of a hard wall for the surface of a trapped BEC, the stability of surface bound states was examined in Kuznetsov and Turitsyn (1988), by imposing that the wave function vanishes at the wall. The corresponding surface potentials, much steeper than harmonic, have been prepared by using laser sheets to trap the dilute quantum gas (for example, in Gaunt et al. (2013)). An inhomogeneous stationary solution of the GPE (the “domain wall”) which coincides with the half of the dark soliton (kink) at rest,Pitaevskii and Stringari (2003) may have as one of its physical realizations a hard wall Kuznetsov and Turitsyn (1988) where localized Bogoliubov excitations were proposed to exist.Pikhitsa (1992) However, the full analytical solution for the corresponding surface-bound excitations has not been found before.At large wavelengths, one class of these excitations represents a surface phonon and the other a ripplon. Our approach is inherently quantum, as it operates near the node plane of the domain wall-soliton, and is hence based on an inherently nonclassical (vector-valued) wavefunction, and is not restricted to large wavelengths, where the (essentially quantum) kinetic terms are small. We note that the existence of a short-wavelength surface excitation (a “surface roton”) was previously conjectured,Reut and Fisher (1971) but its possible connection to capillary waves was then stated as being doubtful. We will see below that for both classes of excitations, starting either from surface phonon or ripplon at large wavelengths, small-wavelength surface excitations exist, with a binding energy approached by both types of excitation at large momenta.

The hard wall boundary condition approximates the steepness of the effective potential at the free surface of liquid helium, which was proven to be composed of a nearly pure condensate of dilute bosonic gas that satisfies the GPE.Griffin and Stringari (1996) The wave function of the BEC is a quantum order parameter that approximately describes the condensate in real liquid helium below the superfluid transition. The helium background (including a well-defined surface) fixes the natural boundary conditions for the BEC. Therefore, the BEC concept accomodates both liquid helium II and a dilute superfluid Bose gas bounded by an external wall.

One may consider the free kink wall with profile $\psi_{0}=\tanh(x)$ extending into the bulk of the liquid $(x\geq 0)$ to model the free surface, demanding only the topological stability of such a solution for which its nodal surface undergoes weak flexural oscillations. Then the position of the hard wall is flexible (like an impenetrable membrane on the surface of helium II) and imitates the free surface of the liquid. The liquid surface of helium II is under these provisos equivalent to a hard wall container.

Here we consider the problem of localized gapless excitation modes by finding analytical solutions of a matrix Schrödinger equation which we show to be equivalent to the BdGE.Chen et al. (1998); Kuznetsov and Turitsyn (1988); Muryshev et al. (1999); Pikhitsa (1992) While recently, Ref. Takahashi et al. (2015) obtained such an analytical solution in the presence of a domain wall, it is restricted to large wavelengths, and furthermore faces the difficulty of extrapolation to the case of an infinite-size surface. We stress that even the classical ripplon (fractional power law) spectrum at small wavenumbers is not trivially obtained from the BdGE, where no classical (phenomenological) surface tension is assumed a priori. In a BEC, the surface tension itself is expressed using Planck’s constant and thus is of an inherently quantum nature.

The binding energy of localized excitations is a primary quantity of interest. Recent experiments that prove the common physical origin of the Landau description of a superfluid and the BEC description Diallo et al. (2014) support the view that the binding energy is relevant. Furthermore, neutron scattering experiments in helium II Prisk et al. (2013) reveal a surface excitation that directly gives the binding energy. Remarkably, we show that the spectrum of surface excitations can be calculated analytically for any wavevector $k$, reproducing the numerical results and with the analytical results obtained for the limiting cases $k\rightarrow 0$ and $k\rightarrow\infty$. We have solved the BdGE for the case of the domain wall (see Eqs. (4.16-4.19) in Takahashi et al. (2015)). The limit of $k\rightarrow\infty$, which in the bulk BEC results in the energy spectrum $\varepsilon={\hbar^{2}k^{2}}/{2m}+\mu$ where $m$ is the mass of the boson and $\mu=gn_{0}$ is the chemical potential while $g$ and $n_{0}$ are the coupling constant and the BEC particle density, respectively, then leads to $\varepsilon={\hbar^{2}k^{2}}/{2m}+\mu-\Delta$.

II Bogoliubov-de Gennes equations

II.1 Basic setup

The GPE of a scalar quantum gas can be written as:Ginzburg and Pitaevskii (1958)

[TABLE]

We introduce dimensionless quantities by measuring distances in units of the healing length $\xi=\hbar/mc$ and energies in units of the “rest mass energy” $gn_{0}=mc^{2}$ where $c=\sqrt{gn_{0}/m}$ is the sound velocity. The stationary version of Eq. (1) for a kink with node at the position $x=0$ gives the wavefunction $\psi_{0}=\tanh(x)$ of the soliton. We will impose perturbations on this solution to investigate its Bogoliubov excitations by representing $\psi$ of Eq. (1) as a sum of plane waves:Pitaevskii (1961) $\psi=\psi_{0}(x)+\vartheta(\vec{r},t)$ with $\vartheta(\vec{r},t)=a_{\omega,\vec{k}}(x)\exp(i\vec{k}\cdot\vec{\varrho}-i\omega t)+b_{\omega,\vec{k}}^{*}(x)\exp(-i\vec{k}\cdot\vec{\varrho}+i\omega t)$, where $\vec{r}=(x,\vec{\varrho})$, $\vec{\varrho}$ lies in the plane orthogonal to the $x$ direction (we consider the situation that all functions decay exponentially with increasingly larger positive $x$), $\vec{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_{\omega,\vec{k}}(x)$ and $b_{\omega,\vec{k}}(x)$. Introducing the functions $\psi_{1}=a+b$ and $\psi_{2}=a-b$, after linearizing Eq. (1) we get a pair of coupled Schrödinger equations:Pikhitsa (1992)

[TABLE]

where $\kappa=|\vec{k}|\xi/\sqrt{2}=k\xi/\sqrt{2}$ and $\varepsilon=\omega$. This pair of equations is identical to the corresponding Bogoliubov-de Gennes equations (see Dziarmaga (2004); Chen et al. (1998)) if one rewrites them for the functions $a$ and $b$. To the best of our knowledge, Eqs. (2) and (3) have never been solved exactly before for arbitrary nonzero $\kappa$ and $\varepsilon$. We find a formal general solution for these equations and illustrate its viability by obtaining a rigorous expression for 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 $\varepsilon_{\rm b}=\kappa\sqrt{2+\kappa^{2}}.$ For $\kappa\to 0$ this gives the bulk phonon dispersion $\varepsilon_{\rm b}\simeq\sqrt{2}\kappa+\kappa^{3}/2\sqrt{2}$ and for $\kappa\to\infty$ it reads $\varepsilon_{\rm b}\simeq\kappa^{2}+1$, which represents a free boson plus chemical potential. The localized excitations to be derived, by definition, have an energy spectrum lying lower than the bulk one.

II.2 Supersymmetry at an exceptional point

We first remark that at the exceptional point of symmetry $\varepsilon=0$ and $\kappa=0$, Eqs. (2) and (3) are the parts of a supersymmetric Hamiltonian with zero ground state energy. Indeed, on introducing the matrix operator

[TABLE]

so that the left-hand side of Eqs. (2) and (3) takes the form of a matrix Hamiltonian

[TABLE]

with its partner Hamiltonian

[TABLE]

we produce a supersymmetric (SUSY) Hamiltonian

[TABLE]

that may canonically be expressed through the supercharges

[TABLE]

as an anticommutator

[TABLE]

The supersymmetry is explicitly broken when either $\varepsilon$ or $\kappa$ (or both) are not zero which, as we will discuss in detail below, leads to a splitting of the SUSY-degenerate ground state into two gapless excitations (a “light” one with $\varepsilon\propto\kappa$ and a “heavy” one with $\varepsilon\propto\kappa^{3/2}$), both bound to the wall.Pikhitsa (1992)

II.3 Boundary conditions

The boundary conditions for $\psi_{1}$ and $\psi_{2}$ in (2) and (3) form two distinct classes. At the node of the kink $\psi=0$, that is both ${\rm Re}~{}\psi=0$ and ${\rm Im}~{}\psi=0$, and therefore also $\psi_{1}=0$ and $\psi_{2}=0$.

However, an additional possibility exists: For $\varepsilon=0$ and $\kappa=0$, Eqs. (2) and (3) have the solutions $\psi_{1}^{(0)}=1-\psi_{0}^{2}$ and $\psi_{2}^{(0)}=\psi_{0}$, the first of which is the so-called “zero mode”,Pikhitsa (1992); Dziarmaga (2004) which leads to Goldstone gapless modes (ripplons and phonons) when the SUSY is broken. This corresponds to a translation of the kink $\psi_{0}$ as a whole along $x$, resulting in the displaced kink $\psi_{0}$ to read as follows: $\psi_{0}(x+\delta x)\simeq\psi_{0}(x)+\psi_{1}^{(0)}\delta x$. Thus the condition ${\rm Re}~{}\psi=0$ turns into $\psi_{0}^{\prime}\delta x(\vec{\varrho},t)+{\rm Re}~{}\vartheta(\vec{r},t)=0$ which determines the shape of the loci of nodes $\delta x(\vec{\varrho},t)$ (the shape of the surface). The derivative of such a mode with respect to $x$ is zero at $x=0$. The mode with the mixed boundary conditions $\frac{d}{dx}\psi_{1}\mid_{x=0}=0$ and $\psi_{2}\mid_{x=0}=0$ allows the “rippling” of the soliton and is thus called ripplon mode.Pikhitsa (1992) As we shall see below, its energy spectrum at low $\kappa$ coincides with the one for a classical capillary wave. The mode with “zero” boundary conditions $\psi_{1}\mid_{x=0}=0$ and $\psi_{2}\mid_{x=0}=0$, which correspond to a flat hard wall will be called surface phonon mode (with a spectrum starting linear).Pikhitsa (1992) Finally, the flat hard wall excludes the possible solution $x\psi_{0}-1$ Muryshev et al. (1999) of Eq. (3) at $\kappa=0$, $\varepsilon=0$ which could lead to the so-called snake instability,111The snake instability amounts to a moving wall (a nodal plane) with its transverse parts moving at different velocities, which is hence acting to destroy the wall cf. Refs. Kuznetsov and Turitsyn (1988); Muryshev et al. (1999). This latter solution does not satisfy zero boundary conditions.

III Asymptotic solutions

We first derive the large and small wavelength solutions of the BdGE, noting that solely the large wavelength case has been considered before.Pikhitsa (1992); Takahashi et al. (2015)

III.1 Large wavelengths

First consider the case of $\kappa\to 0$. For the ripplon spectrum we make an ansatz for $\psi_{1,2}$ in the form of a series in $\varepsilon$: $\psi_{1}\simeq\psi_{1}^{(0)}+\varepsilon\psi_{1}^{(1)}+\mathcal{O}(\varepsilon^{2})$ and $\psi_{2}\simeq\psi_{2}^{(0)}+\varepsilon\psi_{2}^{(1)}+\mathcal{O}(\varepsilon^{2})$. A zeroth-order approximation is the solution of the homogeneous equations Eqs. (2) and (3) with $\varepsilon=0$. This solution can be found for any $\kappa$ (which is verified by direct substitution):

[TABLE]

where $\alpha_{1}=\sqrt{2}\sqrt{2+\kappa^{2}}$ and $\alpha_{2}=\sqrt{2}\kappa$. To determine $\psi_{1}^{(1)}$ and $\psi_{2}^{(1)}$, we have to solve the inhomogeneous equations that follow from Eqs. (2), (3) when $\kappa=0$:

[TABLE]

With the help of the Green functions of the homogeneous equations the inhomogeneous solutions are found as:

[TABLE]

Finally, the derivative with respect to $x$ of $\psi_{1}$ at $x=0$ is found from Eqs. (23), (27) to be $\psi_{1}^{\prime}=A\alpha_{1}(2-\alpha_{1})(2+\alpha_{1})/3+B\varepsilon$, which according to the mixed boundary conditions should be zero together with $\psi_{2}=-A\varepsilon+B\alpha_{2}$, according to Eqs. (24) and (28). A vanishing determinant of the $A,\,B$ linear equations matrix

[TABLE]

gives the ripplon spectrum. Taking into account that $\alpha_{1}\simeq 2+\kappa^{2}/2$ for $\kappa\to 0$ and retaining only the lowest power of $\kappa$, we obtain the fractional dispersion

[TABLE]

The spectrum (32) is shown in Fig. 1. Note that the localization of the ripplon at low $\kappa$ is governed by $\alpha_{2}=\sqrt{2}\kappa$. The spectrum (32) coincides with the well-known expression for the frequency of capillary waves (in the deep water limit), which reads in dimensionful form $\varepsilon=\hbar\sqrt{\sigma/mn_{0}}~{}k^{3/2}$ where $\sigma=\frac{2}{3}\hbar cn_{0}$ is the surface energy density of the stationary soliton $\psi_{0}$. Ginzburg and Pitaevskii (1958) We note that $\sigma$ is exactly half of the energy of the dark soliton at rest [see Eq. (5.59) in Pitaevskii and Stringari (2003)].

Zero boundary conditions lead to surface phonons, for which we obtain the whole spectrum analytically in Sec. IV.1 below. We here only mention in connection to the above discussion that $\alpha_{2}$ for phonons at low $\kappa$ is proportional to $\kappa^{2}$, indicating a much weaker localization as compared to the ripplons.

III.2 Small wavelengths

When $\kappa\rightarrow\infty$, we introduce the function $\chi$ and constant $\Delta$ so that $\psi_{1}=\psi_{2}+\chi/k^{2}$, $\psi_{2}=\psi_{\infty}$ and $\epsilon=\kappa\sqrt{\kappa^{2}+2}-\Delta\asymp\kappa^{2}+1-\Delta$. Then Eqs. (2) and (3) turn into

[TABLE]

with $\chi=-\psi_{0}^{2}\psi_{\infty}$, which after adding and subtracting both equations leads to

[TABLE]

where the hypergeometric function contains $\alpha_{\infty}=\sqrt{2\Delta}$ and $s=(\sqrt{17}-1)/2$ is one of the solutions of the equation $s(s+1)=4$ (see Landau and Lifshitz (1977)). The second solution leads to the same result. The boundary condition $\psi_{2}=0$ at $x=0$ imposes the following identity:

[TABLE]

which demands (for fixed $s$) $1+\alpha_{\infty}-s=0$ in order to have the infinity in the denominator from the corresponfing Gamma function, and therefore $\alpha_{\infty}=(\sqrt{17}-3)/{2}\simeq 0.562$ while $\Delta=\alpha_{\infty}^{2}/2\simeq 0.158$. Finally, the hypergeometric function in (35) reduces to $\psi_{0}$ so that $\psi_{\infty}=\psi_{0}(1-\psi_{0}^{2})^{\alpha_{\infty}/2}=\tanh(x)/\cosh(x)^{\alpha_{\infty}}$ (see Fig.1 (b)) and $\chi=-\psi_{0}^{3}(1-\psi_{0}^{2})^{\alpha_{\infty}/2}$. Therefore, both $\psi_{2}$ and $\psi_{1}$ satisfy the zero boundary conditions. Analogously, one can show222By making use of the known solutions of the homogeneous equations (23) and (24) to satisfy the boundary conditions that the function $\psi_{\infty}$ is also the limiting function for large $\kappa$ in the case of mixed boundary conditions, so that the difference between the functions appears only in close proximity to the boundary $x=0$, at a typical distance $1/\kappa$ [see Fig.1 (c); $\psi_{1}$ deviates from $\psi_{2}$ and hits the $\psi$ axis with zero derivative]. Thus the dimensionful binding energy of the excitation localized near the surface depends only on the bulk parameter $mc^{2}$.

IV Exact solution of the full BdGE

It is well established that many exact solutions of Schrödinger equations with various types of potentials can be directly related to solutions of hypergeometric equations (see, e.g., Ref.Ishkhanyan (2015) for a list); hence factorizations used in quantum mechanics can be obtained from factorizations employing hypergeometric operators. Cotfas and Cotfas (2011) Here, using hypergeometric matrices (which we discuss in detail in the Appendix), we derive below an exact solution of the BdGE.

We aim at finding the exact solution of Eqs. (2) and (3) at arbitrary nondimensionalized momentum $\kappa$. To do so, let us transform these equations into a single matrix hypergeometric equation, where we employ the fact that matrix generalizations of both hypergeometric function and Gamma function were previously shown to be mathematically viable tools.Tirao (2003); Jódar and Cortés (1998) We introduce a wavefunction ansatz by analogy with Eq. (35):

[TABLE]

with a formal parameter $\alpha$. Below this single scalar parameter will be replaced with a matrix, which constitutes the key starting point of finding our exact solution. We now rewrite Eqs. (2) and (3) as

[TABLE]

To turn (38) into a matrix hypergeometric equation, we introduce the vector-function $\hat{\Phi}$, the identity matrix $\hat{1}$, the matrix $\hat{\alpha}$, and matrices $\hat{a},\hat{b},\hat{c}$ derived from it, as follows

[TABLE]

Taking the square root of the matrix $\hat{\alpha}^{2}$ gives, choosing the positive sign,

[TABLE]

where $r=\sqrt{1+\kappa^{2}/2-l^{2}}$, $p=\sqrt{\kappa^{2}/2-l^{2}}$, and $l={\varepsilon}\sqrt{\kappa^{2}+1-(\kappa^{2}(\kappa^{2}+2)-\varepsilon^{2})^{1/2}}/2\sqrt{\varepsilon^{2}+1}$. The two positive eigenvalues of the matrix $\hat{\alpha}$ are

[TABLE]

and $\exp(-\alpha_{2}x)$ determines the asymptotic decay of $\psi_{1}$ and $\psi_{2}$ as $x\to\infty$ (corresponding to the lower sign above). After introducing the matrices, Eq. (38) becomes the canonical Gauss hypergeometric equation in matrix form

[TABLE]

where primes mean differentiation with respect to $z$.

IV.1 Surface phonons

Eq. (54) is formally solved by Eq. (68) contained in the Appendix. We can then obtain the spectrum of the surface phonon localized near a flat hard wall as follows. The boundary condition at $x=0$ (that is at $z=1/2$) will be fulfilled when $\hat{\Phi}=0$, which demands that the determinant of matrix function (70) be zero. The spectrum is then given by the equation

[TABLE]

where the matrix $\hat{P}$ is derived in the Appendix, see Eq. (74). Eq. (55) reproduces the spectrum calculated before for the two limiting cases $\kappa\to 0$ and $\kappa\to\infty$ in section III. When $\kappa\to 0$ the spectrum is $\varepsilon=\sqrt{2}\kappa+\mathcal{O}(\kappa^{5})$, so that the $\kappa^{3}$ term is missing, while the bulk phonon starts with higher energy as $\varepsilon_{\rm b}=\sqrt{2}\kappa+\kappa^{3}/2\sqrt{2}+\mathcal{O}(\kappa^{5})$. Let us define the binding energy as $\delta\varepsilon=\varepsilon_{\rm b}-\varepsilon$. Then the latter starts as $\kappa^{3}/2\sqrt{2}$ [see Fig. 2 (c)]. Now let us consider the other limit $\kappa\to\infty$. It is easy to see that seeking the solution in the form $\varepsilon\asymp\kappa^{2}+1-\Delta$ leads to $r=p\asymp\kappa/2+\sqrt{2\Delta}/4$ and $l\asymp\kappa/2-\sqrt{2\Delta}/4$ which after substitution into Eq. (55) give $2\Delta+3\sqrt{2\Delta}-2=0.$ This has the same root as found before from Eq. (36), $\sqrt{2\Delta}=(\sqrt{17}-3)/2=\alpha_{\infty}\simeq 0.562$ and therefore $\delta\varepsilon_{\infty}=\Delta=\alpha_{\infty}^{2}/2\simeq 0.158$.

The coincidence with the exact asymptotic results obtained in section III confirms the correctness of Eq. (55). Note that the slower decay exponent $\alpha_{2}$ can be approximated by a simple expression $\alpha_{2}=\kappa^{2}/(1+\kappa^{2}/\alpha_{\infty})$ that fits the exact expression of Eq. (53) with $\varepsilon$ from the exact solution of Eq. (55) within $0.2\%$. We plot the decay exponential in Fig. 2(d) in a broad range of wavenumbers.

Finally, it is interesting to note that the first order approximation in powers of $z=1/2$ in the case of the flat hard wall boundary conditions for surface phonons can be represented by the equation

[TABLE]

which, distinct from the case of ripplons discussed below, accidentally gives the exact spectrum of Eq. (55), where $\hat{U}_{1}$ is defined below in Eq. (67).

IV.2 Ripplons

Now consider the case of mixed boundary conditions corresponding to ripplons. Let us first rephrase the general form of the solution (68) provided in the Appendix in the form of a vector function. With the help of (70), we get

[TABLE]

where $A$ and $B$ are arbitrary constants.

Writing the surface phonon boundary conditions explicitly,

[TABLE]

the spectrum for surface phonons is obtained after equating the determinant of the above equation for $A,\,B$ to zero at $z=1/2$. We then proceed analogously as for this case of the flat wall for the rippled wall, except that we differentiate with respect to $z$ the first equation for $\phi_{1}$.

Expanding the hypergeometric function of Eq. (68) and its $z$ derivative to first order in $z=1/2$, one gets from Eq. (LABEL:eq:s19)

[TABLE]

where $\hat{U_{1}}$ and $\hat{U_{2}}$ are calculated according to Eq. (69) for $m=0,1$ with $\hat{\alpha}$ and $\hat{a}\hat{b}$ taken from Eqs. (52) and (49), respectively:

[TABLE]

IV.3 Comparison to numerics

The energy spectrum and binding energy–decay parameter for ripplons in a broad range of wavenumbers are shown in Figs. 3 (a),(b), respectively, in comparison to their values obtained with numerical solutions of the differential equations Eqs. (2) and (3), represented by the symbols. For the numerics, we used PTC’s MathCad 11, applying proper boundary conditions at the surface, imposing an exponentially fast decay at infinity (the latter leads to underestimate the binding energy values, see for a discussion below). One can see that even to lowest nontrivial order in the series on $z=1/2$ [see Eq. (64)], the results shown by solid lines in Fig. 3 (a),(b) are rather close to the numerical solutions.

On the other hand, for the surface phonon, the numerical results can be rendered closer to the exact spectrum from (55), as displayed in Figs. 2 (a),(b), although a slight systematic deviation is still noticeable. These deviations stem from the fact that the numerical solution of the differential equations relies on the criterium of localization: the solution should decay into the bulk, implying that another boundary condition is that the wavefunction should approach zero at infinity. Numerical calculations are imposing boundary conditions at a finite distance, however large. The numerics therefore slightly exaggerates the decay; hence the numerical energy is slightly lower than the exact energy at a given wavenumber.

V Conclusion

In summary, starting from the matrix hypergeometric equation (54), we obtained its formally exact solution at the boundary (70). Many exactly solvable Schrödinger equations with various potentials have solutions of the hypergeometric variety. Often there are also supersymmetric partners in the Hamiltonian operator, as in the case of a hydrogen atom with its Coulomb potential. When a continuous wavenumber is present, creating a bandgap structure, gapless states that stem from (or are accompanied by) Goldstone zero energy modes may exist.Hasan and Kane (2010) In the case of the BdGE that we considered here, we were not only able to obtain an exact solution, but also to express the dispersion relation of the Bogoliubov surface excitations for surface phonons in the closed form of the algebraic equation (55). We have furthermore shown that for ripplons, even a lowest nontrivial order truncation of the hypergeometric series produces results close to numerical solutions of the BdGE.

We now discuss the relation of the analytically obtained binding energy of surface phonons to the experimental finding of Ref.Prisk et al. (2013) for helium II confined by the hard walls of cylindrical pores. Even though the present BEC model with contact interactions does not reproduce the roton minimum in the bulk dispersion curve, it provides a correct estimate for the binding energy in the low-density surface region. Indeed, the binding energy (the difference between bulk and surface excitation energies) has been measured to be $0.15$ meV at the roton-region wavevector $k=2/{\textup{\AA}}$.Prisk et al. (2013) The latter corresponds to $\kappa\simeq 1$, using the estimate $\xi=\hbar/mc\simeq 0.7\textup{\AA}$, with the “bulklike” speed of sound $c=228\,$m/s at full pore.Prisk et al. (2013) We can read off Fig. 2(a) a dimensionless binding energy of approximately 0.07 at $\kappa\simeq 1$, which agrees to good accuracy with the experimental value (using that $mc^{2}\simeq 2.2$ meV). This agreement was obtained at small wavelengths, to which previous approaches did not apply, and which in the bulk correspond to the roton minimum. Therefore, while the latter bulk dispersion feature is not accurately described by our mean-field model, we conclude that the quantum mechanism of trapping excitations close to a surface gives the correct magnitude of the binding energy. We note that the quantum mechanism of binding surface excitations occurs in the solid-state physics of electrons as well, where the bound states are called Tamm and Shockley states.Tamm (1932); Shockley (1939)

The present method for exactly solving the Bogoliubov-de Gennes equations is potentially also useful in more sophisticated cases than the presently considered one. Further extensions of the present approach are for example conceivable by incorporating effectively nonlocal interactions modelling rotons, which occur in dilute quantum gases dominated by dipole-dipole interactions.Fischer (2006) Furthermore, it would be of interest to investigate to which extent the present matrix hypergeometric equation approach can be applied to other physical systems of current widespread interest. For instance, to topological insulators, superconductors, and even to exotic topological mechanical materials.Kane and Lubensky (2013)

Acknowledgements.

The work of PVP was supported by the Global Frontier Center for Multiscale Energy Systems funded by the National Research Foundation of Korea (NRF) Grant No. 2012M3A6A7054855. URF has been supported by the NRF under Grant No. 2017R1A2A2A05001422. *

Appendix A The matrix-valued hypergeometric function

Equation (54), taking the canonical form of a hypergeometric equation, has a formal solution as a matrix-valued hypergeometric function of Gauss Tirao (2003)

[TABLE]

where $(\hat{a},\hat{b},\hat{c})_{0}=\hat{1}$, and higher matrix coefficients are

[TABLE]

In the above equation, the matrices are ordered in a specific way, taking into account their generally noncommutative nature. The related intricate mathematical questions have been discussed in detail when introducing the matrix hypergeometric function in Tirao (2003). Yet this noncommutative nature is not important for our purpose of obtaining the dispersion relations, inasmuch as we deal to this end with the determinant of the matrix hypergeometric function. The latter determinant is expressed below through products of the matrix Euler Gamma function and its inverse in Eq. (70). The Euler Gamma function itself is in turn a product of matrices and their inverse according to Eq. (A4). The determinant of the matrix hypergeometric function is hence independent of the order of the matrices occurring in it.

The hypergeometric function at $z=1/2$ can be expressed through the matrix Gamma functionJódar and Cortés (1998) [because $\hat{c}=(\hat{1}+\hat{a}+\hat{b})/2$, see Eq.(46)], so that

[TABLE]

Then the condition of $\hat{\Phi}=0$ imposed for surface phonons in Sec. IV.1 implies that the matrix (70) has an eigenvalue zero and that therefore its determinant vanishes. We then use that matrix Gamma functions can be represented as follows Jódar and Cortés (1998)

[TABLE]

The role of the matrix $\hat{M}$ is played by either $\hat{1}+\hat{a}$ or $\hat{1}+\hat{b}$. Because the determinant of (70) is required to vanish, the determinant of a Gamma function in the denominator should be infinite. By (71), this is only possible if either the determinant of $\hat{1}+\hat{a}$ or that of $\hat{1}+\hat{b}$ is zero.

One can prove that the determinant of either $\hat{1}+\hat{a}$ or $\hat{1}+\hat{b}$ being zero gives the same spectrum. However an analytical solution for the matrix equations (46), (49) for $\hat{a}$ and $\hat{b}$ is difficult. To obtain analytical results we instead utilize the product $\hat{P}=((\hat{1}+\hat{a})(\hat{1}+\hat{b}))/2=(\hat{1}+\hat{a}+\hat{b}+\hat{a}\hat{b})/2$. We readily get the matrix $\hat{P}$ from Eqs. (46),(49) and (52):

[TABLE]

and taking $\det\hat{P}=0$ yields Eq. (55) of the main text.