Density profile of a semi-infinite one-dimensional Bose gas and bound states of the impurity

We study the effect of the boundary on a system of weakly interacting bosons in one dimension. It strongly influences the boson density which is completely suppressed at the boundary position. Away from it, the density is depleted over the distances on the order of the healing length at the mean-field level. Quantum fluctuations modify the density profile considerably. The local density approaches the average one as an inverse square of the distance from the boundary. We calculate an analytic expression for the density profile at arbitrary separations from the boundary. We then consider the problem of localization of a foreign quantum particle (impurity) in the potential created by the inhomogeneous boson density. At the mean-field level, we find exact results for the energy spectrum of the bound states, the corresponding wave functions, and the condition for interaction-induced localization. The quantum contribution to the boson density gives rise to small corrections of the bound state energy levels. However, it is fundamentally important for the existence of a long-range Casimir-like interaction between the impurity and the boundary.


I. INTRODUCTION
The boundaries play an important role in one-dimensional quantum liquids affecting, e.g., correlation and response functions, as well as the ground-state energy [1][2][3][4][5][6][7][8]. The boundary has also an impact on the density profile of particles n(x), which is suppressed at the boundary position x = 0. Away from it, the particle density in fermionic systems shows socalled Friedel oscillations [9]. They describe an oscillatory decay of n(x) − n 0 , where n 0 is the average density. The envelop of Friedel oscillations follows 1/x K law for spinless fermions [3,10]. Here K is the Luttinger liquid parameter, which is determined by the interaction between fermions. The periodicity of oscillations is controlled by the average fermion density.
Friedel oscillations in fermionic systems are studied within the harmonic Tomonaga-Luttinger liquid description [3,10]. This is the low-energy theory of both, interacting fermions and bosons in one dimension [11]. It was applied for bosons in Ref. [12] where the same pattern of density oscillations was found for separation longer than the inverse mean density of bosons, 1/n 0 . Applied to the special case of weaklyrepulsive bosons where K 1, the result of Ref. [12] implies very rapid saturation of n(x) − n 0 . Such fast recovery of the density is not expected to occur in weakly interacting superfluids where the density should change on the scale comparable to the healing length ξ ∼ K/n 0 , which is much longer than 1/n 0 . Thus, in order to describe the density profile of bosons one needs to go beyond the harmonic Tomonaga-Luttinger liquid theory.
In this paper we study weakly-interacting bosons in a semiinfinite system. We consider the equation of motion for the field operator. At the mean-field (classical) level it reduces to the Gross-Pitaevskii equation [13]. Its solution reveals that the spatial extent of the depletion of the boson density imposed by the boundary is controlled by the healing length ξ. At distances longer than ξ, the mean-field density exponentially rapidly reaches the constant value n 0 . Accounting for the effect of quantum fluctuations around the mean-field solution within the Bogoliubov-de Gennes formalism, we find that the density reaches n 0 much slower, following 1/x 2 law.
Our approach enables us to obtain an analytic expression for the density profile at all distances.
The nonuniform density profile of the Bose gas determines the potential for a weakly-coupled quantum impurity introduced in the system. For repulsive interaction between the impurity and the bosons, the impurity can be localized. We solve the Schrödinger equation and characterize the impurity by the energy spectrum of the bound states, their wave functions, and the mean position. We find the condition for the appearance of the bound states, which is a threshold for a single dimensionless parameter that involves the masses of the impurity and of the particles of the Bose gas as well as the interaction strengths. We note that a related phenomena of self-trapping of a single impurity in homogeneous Bose-Einstein condensates was studied in Refs. [14][15][16][17][18]. Similarly, the localization of bosonic atoms by fermionic ones in attractive Bose-Fermi mixtures was examined in Ref. [19], while the formation of bounds states of two impurities immersed in one-dimensional liquid has been studied recently in Refs. [20][21][22]. The formation of bound states of two polarons in a Bose-Einstein condensate is considered in Ref. [23].
The paper is organized as follows. In Sec. II we introduce the model of interacting bosons in a semi-infinite geometry. In Sec. III we solve the equation of motion for the single-particle field operator. We find the mean-field solution and the first two quantum corrections. This enables us to evaluate the density of bosons, including its quantum contribution, in Sec. IV. In Sec. V we consider the problem of a mobile quantum impurity interacting with the system. We solve the Schrödinger equation for the impurity and study its properties. Section VI is devoted to discussions, while in Sec. VII we summarize our results. In Appendix A we present a simplified procedure that leads to the density in the regime of large separations from the boundary.

II. THE SEMI-INFINITE BOSE GAS
We are interested in the influence of the boundary at x = 0 on the physical quantities in a one-dimensional system of interacting bosons. We thus study a long system described by Here m is the mass of bosons, while the coupling constant g > 0 describes the repulsive contact interaction between them. The system size is L; however we will eventually consider the thermodynamic limit. The bosonic single particle operatorŝ Ψ andΨ † satisfy the usual equal time commutation relations , while the other commutators vanish. The model (1) is characterized by the dimensionless parameter [24] γ = mg/ 2 n 0 , where n 0 is the mean particle density. The boundary of the system imposes the nullification of the single-particle operator at its position, The condition (2) implies that the boson density also vanishes at x = 0.
Our first goal in this work is to calculate the density profile of the bosons where the average is with respect to the Hamiltonian (1). For simplicity, we introduce the dimensionless coordinates for the position and the time, respectively, defined as Here ξ µ = / √ mµ denotes the healing length, while µ is the chemical potential. Assuming the single particle operator in the formΨ its equation of motion i ∂ tΨ = [Ψ, H] in the dimensionless units becomes We will solve Eq. (6) at γ 1, corresponding to the regime of weak interaction. In this case, one can expand the field operator as [13,25] ψ(X, T ) = ψ 0 (X) + αψ 1 (X, T ) + α 2ψ 2 (X, T ) + . . . , (7) such that [ψ(X, T ),ψ † (X , T )] = α 2 δ(X − X ). Here the small parameter is α = (γgn 0 /µ) 1/4 ≈ γ 1/4 1. The function ψ 0 (X) describes the time-independent wave function of the system in the absence of fluctuations. The field operatorsψ 1 andψ 2 account for its first and the second quantum correction. Note that the Bose-Einstein condensate does not exist in one dimension in the thermodynamic limit due to strong effect of long-wavelength fluctuations. However, in finite-size systems, the inverse system size provides an infrared cutoff. The perturbative expansion (7) is justified as long as [25,26] ln(L/ξ µ ) 1/ √ γ where L is the length of the system. The latter inequality shows that for weakly interaction bosons [i.e., at γ 1] the system size can actually be huge. We point out that the density n(x) [see Eq. (3)] which is to be calculated, is independent of the cutoff and therefore our result applies in the thermodynamic limit as well [27]. The form (7) substituted into Eq. (6) leads to a hierarchy of equations controlled by the small parameter α 1. The equation of motion for ψ 0 (X) is obtained at order α 0 . It reads where we have introduced the operator The expression (8) is known as the Gross-Pitaevskii equation [13]. This second order differential equation should be supplemented by two boundary conditions. One of them follows from Eq. (2) and becomes ψ 0 (0) = 0. The other condition arises from the physical requirement that the density, which is proportional to |ψ 0 (X)| 2 , becomes constant at long separations from the boundary. The real solution [28] of Eq. (8) satisfying such boundary conditions is At weak interaction the chemical potential of the Bose gas is µ = gn 0 . This gives the density (3) at the mean-field level to be n(x) = n 0 tanh 2 (x/ξ), ξ = 1/n 0 √ γ.
This expression shows that the density quickly saturates at distances beyond the healing length ξ (which is ξ µ in the limit γ 1), provided one takes only the leading order term from the expansion (7).

B. First quantum correction
We now consider the effects of quantum fluctuations and determineψ 1 . Its equation of motion is We seek a solution forψ 1 (X, T ) as a superposition of excitations of energy k using the ansatz based on Bogoliubov transformation [13] Since Eq. (12) is linear and homogeneous, we must account for the normalization factor N k in Eq. (13), which should be determined in such a way to satisfy the proper commutation relations betweenψ 1 fields. This is discussed further below. The bosonic operatorsb k andb † k satisfy the standard commu- The boundary condition (2) at order α is given in terms of u k and v k by Substitution of Eq. (13) into Eq. (12) leads to two coupled equations for u k and v k . They are known as the Bogoliubovde Gennes equations and are given by In order to simplify them, we introduce the functions This enables us to decouple the equations and lead to the fourth order differential equation Note that D is given in terms of S as Four independent solutions of the forth-order equations (16) are [29] S n (X) = (−ik n + 2 tanh X)e iknX , (17a) where n ∈ {1, 2, 3, 4}, while k = k 2 + k 4 /4 is the energy dispersion. The four roots entering S n are k 1,2 = ±k and Note that S 4 and D 4 do not appear since they diverge at large X. The two unknown coefficient A and B are determined using the boundary condition (14) that is terms of S and D becomes S(k, 0) = D(k, 0) = 0. We find A = 1 and B = 0. Therefore the general solution of Eqs. (15) that satisfy the boundary condition (14) are the real functions The normalization in Eq. (13) is obtained by requiring [13] This leads to N k = (ξ µ /4L k ) 1/2 at L ξ µ . One can then verify the equal time commutation relation [ψ 1 (X),ψ 1 (Y )] = 0. The evaluation of the other commutation relation is more involved: since δ(X + Y ) always equals zero for X, Y > 0. The second equality in Eq. (21) is obtained after performing the integral , which turns out to be zero. The remaining part after integration over k then gives the delta function in Eq. (21). There we use k>0 (· · · ) = (L/πξ µ ) ∞ 0 dk(· · · ).

C. Second quantum correction
The first quantum contribution to the density (3) is proportional to α 2 and thus is determined by the first two corrections of the field operator in Eq. (7). We thus now consider the second quantum correction to the field operatorψ denoted byψ 2 . Its equation of motion is obtained from Eq. (6) at order α 2 . Since ψ 0 (X) is real [see Eq. (10)], it is sufficient to consider the real part of the expectation value ψ 2 which enters into n(x). For this purpose we introduce the notation ψ 2 = Re ψ 2 . The equation of motion for ψ 2 is where At zero temperature one has Note that the source term f (X) in Eq. (22) is timeindependent and thus ψ 2 is only a function of X. We note that evaluation of f (X) requires a small-k cutoff (λ ∼ ξ µ /L) since it is divergent at k → 0.
The solution of the linear equation (22) can be expressed as where G is the Green's function of the operator L 3 (X). It satisfies Moreover, the Green's function is symmetric G(X, Y ) = G(Y, X) and satisfies to account for the δ(X −Y ) on the right hand side of Eq. (26).
Here Y ± = lim δ→0 + Y ± δ. Also, the Green's function obeys G(0, Y ) = 0 due to the boundary condition imposed by the end of the system, i.e., ψ 2 (0) = 0. For the solution of Eq. (26) we find Therefore, the solution of Eq. (22) for ψ 2 (X) is obtained from Eq. (25), where one should substitute the Green's function and the source function (23). Since the latter depends on the infrared cutoff, ψ 2 has also this feature. In Appendix A we derive ψ 2 (X) in the regime X 1.

IV. LOCAL DENSITY
Having solved the equation of motion for the terms ψ 0 ,ψ 1 , andψ 2 of the single-particle operator (7), we are now in position to evaluate the spatial density profile of the semi-infinite Bose gas. It is given by the expansion Here the mean-field contribution is while the quantum contribution to the density has the form or more explicitly In Eq. (32) one should use Eqs. (19) and (28), while we recall k = k 2 + k 4 /4. Unlike ψ 2 that requires an infrared cutoff, the density contribution (32) does not. Each of the two terms in Eq. (31) is divergent, however their sum is finite. This should be the case since the density fluctuations are expected to be finite unlike the fluctuations of the phase of the single particle operator. Note that a similar situation also occurs in an infinite system [27]. The evaluation of Eq. (32) is rather tedious. We first perform the integration over Y using various trigonometric identities and the integration by parts. As a result, we obtain a cumbersome expression involving several hypergeometric functions. Such expression is an integrable function of k.
To perform the integration we split the integrand into several summands. We should stress that, unlike the whole integrand, the summands can be nonintegrable due to singularities and one should regularize them by adding and subtracting some functions. As a result of such procedure, we obtain the local density that can be expressed in the form where ξ µ = / √ mµ and Here 2 F 3 is the hypergeometric function while I ν (z) and L ν (z) are the modified Bessel function of the first kind and the modified Struve function, respectively.
We now analyze the behavior of the local density (35). At short separations, x ξ, the quantum fluctuation contributions to the density can be neglected and the leading contribution is given by This classical result originates from the mean-field contribution (30) to the density. At longer separations, x ξ, the classical result (36) quickly saturates to a constant n 0 . However, the approach of the local density toward n 0 is actually much slower than one obtains from the mean-field result (36). At x ξ the quantum contribution (31) is the dominant term in the density deviation n(x)−n 0 since it decays algebraically with the distance. This follows from the asymptotic expansion Therefore, at x ξ the leading correction to the local density is given by The crossover distance x c between the two regimes can be estimated by equating the two expressions (36) and (38). As a result we get Note that the crossover distance depends only logarithmically on γ and practically is on the order of few ξ. This is illustrated in Fig. 1 where the density profile of the Bose gas is shown as a function of the distance from the end of the system at several values of γ. The density in the crossover regime between the classical and the quantum one is given by Eq. (35). A simplified derivation of n (1) (X) at X 1 is given in Appendix A. We finally notice that the ground-state energy calculation of Ref. [30] leads to the limiting forms (36) and (38).

V. INTERACTION-INDUCED LOCALIZATION OF A SINGLE IMPURITY
In this section we study the problem of a quantum impurity in the semi-infinite Bose gas. We consider the impurity that is locally coupled to the boson density (35). In the case of repulsion, the density profile of the Bose gas forms an attractive potential for the impurity, which can lead to the bound states in the spectrum. In that case the particle becomes localized by the surrounding medium.
The impurity wave function is a function of x ant t and is governed by the Schrödinger equation. After separation of variables it reduces to the eigenvalue problem on the positive semi-axis x > 0 Here M is the impurity mass, G > 0 denotes the coupling constant, while the density of the Bose gas n(x) is given by Eq. (35). An alternative formulation of the eigenvalue problem (40) is to consider unrestricted x with the symmetric potential Gn(|x|), and to account only for odd eigenfunctions. They have a node at x = 0 and thus automatically satisfy the boundary condition (40b). In order to simplify the notations, we rewrite Eq. (40) in the form where we have introduced In the following we will find the bound state spectrum for the eigenvalue problem (41).

A. Bound states of Pöschl-Teller potential
In the limit where the quantum correction to the density h(z) is neglected, the potential of Eq. (41a) describes a hole of modified Pöschl-Teller type which admits an exact solution [31,32]. We now derive the bound states for this special eigenproblem defined by From the definition (43) follows λ > 1 and thus the potential −λ(λ − 1)/ cosh 2 z in Eq. (45a) is negative. We also notice that for unrestricted z, the eigenstates in such symmetric potential can be even or odd functions of z. For the problem (45) only odd solutions have a node at z = 0. They are of the form [32] f (z) = cosh λ z sinh z where 2 F 1 is the Gauss's hypergeometric function.
In order to discuss the appearance of bound states, we consider the asymptotic expansion of Eq. (46) in the limit z → ∞, which is given by We notice that there are two exponential terms ∝ e ±κz at z → ∞. Taking into account the definition of κ given by Eq. (44), we conclude that a bound state can be realized only if κ 2 > 0, i.e., if Gn 0 > E. In the opposite case, Gn 0 < E, the solution at long separations from the boundary is an oscillating trigonometric function since κ is purely imaginary. In the following we focus on the localized (bound) states and derive their energy spectrum. We assume without loss of generality that κ > 0. Thus, the coefficient in front of the term e κz must vanish in order to obtain a nondivergent (normalizable) wave function at large z. This occurs when the arguments of the Gamma functions in the denominators are 0 or negative integers. Since λ > 1 and by assumption κ > 0, only the argument of Γ 2−λ+κ 2 matters. This yields where η = 0, 1, 2, . . .. From the condition κ > 0 we obtain that integer η satisfies Equation (49) determines the condition λ > 2 necessary for the appearance of the first bound state [33]. The number η denotes the number of nodes of the odd bound-state eigenfunction in the region of interest z > 0. Thus, the localized solutions of the eigenvalue problem Eq. (45) are given by Eq. (46) where κ satisfies Eq. (48), while η is a nonnegative integer satisfying the condition (49). The latter enables us to reexpress the eigenfunction (46) in the form a, b; c; z). The preceding expressions enable us to find the spectrum of bound states of the impurity defined by Eq. (40) provided the system has bound states. The condition for that is sufficiently strong coupling between the impurity and the Bose gas, which should satisfy The latter follows from Eq. (49). The set of allowed values for η (non-negative integers) in the spectrum (51) are determined by the inequality (52). In particular, the first bound state (corresponding to η = 0) occurs for GM /gm > 1. The second one (corresponding to η = 1) exists for GM /gm > 6, etc.
In Fig. 2 we show the bound state energies for the first five levels. We notice that the specific bound state energy expressed in units of Gn 0 , which is the value of the potential at infinity, is a function of the single parameter GM/gm. In Fig. 3 we show the normalized wave functions for the bound state levels for the specific value GM/gm = 30. We notice that the wave functions for energies near the top of the potential, Gn 0 , are weakly localized as their spatial extension increases. The level (quantum) number η corresponds to the number of nodes of the wave function. For a given wave function, its spatial extension decreases with increasing GM/gm, since in that case the potential becomes deeper. This is reflected in the mean distance of the particle from the origin, We plot this dependence in Fig. 4. Close to the threshold for the appearance of the bound state, the mean distance diverges according to the law The denominator of Eq. (54) denotes the distance from the threshold for the appearance of the bound state measured from  Notice that the wave function corresponding to η = 3 is spatially less localized since the value of GM/gm is close to the threshold for its appearance.
the "localized" side where the particle is in the bound state.
Indeed, x η determines the localization length. This easily follows from the asymptotic expansion of the wave function [see Eq. (47)], where the localization length is 1/κ close to the "transition". The localization length diverges with the first power of the distance from the "transition", which physically denotes the disappearance of the bound state as GM/gm is decreased.
In the preceding part we considered a localized impurity, i.e., its bound states which exist when the eigenenergies are smaller than Gn 0 (corresponding to the value of the potential at infinity). We should stress that the impurity has a continuum of scattering states with energies higher than

B. Quantum correction to the potential
In the previous subsection we have solved the eigenproblem for the impurity using the mean-field contribution to the potential. We now account for the quantum correction h(z), see Eq. (41). We were not able to solve analytically the whole eigenproblem, since the form of h(z) is very complicated. However we can take advantage of the small parameter √ γ that controls it and find the corrections of the bound state energy (51) using perturbation theory. The energy correction is given by where h is defined in Eq. (34), while f η by Eq. (50). The quantum correction √ γ Gn 0 h(x/ξ) to the effective potential [see, e.g., Eq. (40a)] becomes practically important at distances x ξ (see Fig. 1) with respect to the classical one −Gn 0 / cosh 2 (x/ξ). Contrary to the function h(x/ξ) that has a long-range tail ∝ −1/x 2 , the wave function is localized and exponentially small at x ξ. As a result the numerator in Eq. (55) is small thus the quantum correction to the bound state energy is negligible. Notice that very near the threshold for the bound state appearance (see, e.g., the curve for η = 3 in Fig. 3) the wave function may have some overlap with the tail of h(x/ξ), and the correction δE η can be somewhat important with respect to E η − Gn 0 . Nevertheless, δE η is negative since the quantum correction to the potential broadens the well. The quantum corrections to the different energy levels of the bound states are numerically evaluated using Eq. (55) and shown in Fig. 5 [34].

VI. DISCUSSION
We have calculated the bound states of the impurity in a semi-infinite Bose gas. We assumed that the density of the Bose gas it not affected by the presence of the impurity [cf. Eq. (40)]. This approach is valid in cases of weak coupling G between the impurity and the particles of the Bose gas. A more general model would be the Hamiltonian (1) supplemented by the impurity part From the solution of the whole problem one could obtain the spectrum of the bound states at any G. However this is in practice difficult to do analytically. The parameter region where our approach applies can be obtained by calculating the correction to the Bose gas density (35) due to the coupling with the impurity. For a heavy impurity such correction is small at G g/ √ γ [30]. In this regime is therefore justified to study the simplified problem (45). Notice that for weakly interacting bosons we have γ 1 and thus our work covers a big range of values for G.
The quantum contribution to the boson density that universally behaves as 1/x 2 at x ξ [see, e.g., Eq. (38)] leads to a small energy correction to the bound state levels of the localized particle. However, the knowledge of such spatial dependence of the density of particles far from the boundary is fundamentally important for quantities where the density gradient matters. An example is the Casimir-like interaction [27,[35][36][37][38] between the impurity and the boundary of the system mediated by density fluctuations of the medium. In the regime x ξ its form was found in Ref. [30], while our result (35) determines the interaction law at all distances. It takes the form It is expected that under the influence of the long-range potential (57), a heavy impurity immersed in the system far from the boundary slowly drifts towards it and eventually get localized. During this process, the impurity energy that equals Gn 0 at long distances is dissipated by exciting the Bose gas. The boundary in our system can exist naturally as the system's end. Alternatively, it can be created by a heavy impurity strongly coupled to the Bose gas that creates impenetrable potential and leads to the complete depletion of the boson density at its position. Thus the effective interaction between the boundary and the impurity can also be interpreted as the Casimir-like interaction between two very different impurities. We notice two related very recent experimenal works [39,40], which have demonstrated the existence of the induced interaction between quantum particles mediated by the surrounding quantum gas.
The Hamiltonian (1) can be studied within the dual model of attractive fermions [41,42]. The density of bosons in the original model is equal to the density of fermions in the dual one. Such mapping is particularly useful in the regime of strong interaction between bosons, corresponding to γ 1 and the Luttinger liquid parameter K 1 + 4/γ. In this case the fermions are weakly-attractive, which enables one to study their density and obtain the characteristic form with Friedel oscillations. It would be then interesting to study the fate of Fridel oscillations as the value of K increases and find a connection with our result (35). A numerical study of finite number of bosons [43] also suggests the above picture. We finally note that the model (1) is integrable in the box geometry [5] and thus it is in principle possible to use the exact Bethe ansatz solution to find the exact density profile in the thermodynamic limit. We are not aware of such study in the literature.

VII. SUMMARY
In this paper we first calculated the density profile of the weakly-interacting semi-infinite one-dimensional Bose gas. We solved the equation of motion for the field operator. At the mean-field level, it leads to the density profile (11). Taking into account the effect of quantum fluctuations around the mean-field solution we calculated the quantum contribution to the density, see Eqs. (34) and (35). It shows the universal 1/x 2 behavior at long distances.
We then studied the spectrum of the bound states of a quantum impurity which experiences the potential created by the surrounding Bose gas. We were able to exactly solve this problem accounting for the mean-field boson density. The discrete spectrum is given by Eq. (51), while the corresponding wave function are written in Eq. (50). Bound states exist for sufficiently strong repulsion between the impurity and particles of the Bose gas, see Eq. (52). We also found how the mean particle distance from the boundary behaves, and in particular its divergence (54) as one approaches the threshold for the appearance of the bound state levels. We showed that the quantum contribution to the density gives rise to the negligible corrections of the bound state levels. However, the quantum contribution to the density is important since it leads to the long-range Casimir-like interaction between the impurity and the boundary.