Conditions for Bose–Einstein condensation in periodic background

We investigate Bose–Einstein condensation of a noninteracting gas of Bose particles moving in the background of a periodic lattice of delta functions. In the one-dimensional case, where one has no condensation in the free case, this property persists also in the presence of the lattice for all examples which are considered in the present paper and we could only formulate some conditions which are necessary for condensation. We also considered the three-dimensional case and showed that the lattice does not destroy condensation. We calculated, for small coupling, the change in the critical temperature, which is lowered by the lattice. Finally, we took another, more general view on the problem using heat kernel expansion, and discuss BEC for Casimir effect related configurations.


Introduction
The investigation of Bose-Einstein condensation (BEC) in one-dimensional configurations has a long history. In [1] it was shown that long-range order in homogeneous Bose and Fermi systems in one and two dimensions is not possible. A different picture appears in confined one-dimensional systems, with confining potential of finite volume [2], or with interacting bosons (see the more recent [3] and references therein).
In the present paper, we investigate the question in which way a periodic background influences BEC in an otherwise homogeneous background. Such systems are assumed to have an infinite spatial extension such that no boundary effects occur and no confining potential. As background, we take primarily lattices of delta functions, which in one dimension represents the well known Kronig-Penney model. Delta function potentials have the advantage that much more explicit results are possible. As known, in more than one dimension these must be handled with care which was considered in detail, for example, in [4].
Similar systems were recently investigated in connection with the Casimir effect in [5], where the interaction of such lattices in one and two dimensions was calculated, including several limiting cases, for dense or dilute lattices. The methods used that paper will be helpful below.
As for the Casimir effect and entropy, there is a yet not completely settled problem with negative entropy. It consists in the observation that entropy, calculated as negative derivative of the free energy of the electromagnetic field in the presence of polarizable bodies takes negative values. This effect was first observed in the interaction of two dispersive slabs [6], later also in less trivial geometries like a ball in front of a plane (see, for example, table 1 in [7]). In these examples, it was the interaction part of the free energy in the presence of the considered bodies, which was of interest for the Casimir effect at finite temperature. At once the related entropy was calculated and negative values were encountered. However, it was recognized that this way one has only the interaction part and not the whole entropy which still could well be positive.
To clarify this question, past years the free energy and the entropy of single, free-standing objects was calculated, first for a plane carrying a delta function potential [8], later for a similar sphere [9]. Since that work was plagued by some divergences, in [10] for a sphere and in [11] for a plane it was recognized that in the temperature-dependent part of the free energy and in the entropy there are in fact no divergencies and all calculations can be carried out straightforwardly. This way, in the mentioned examples the occurrence of negative entropy was confirmed. It must be mentioned that with these calculations there is still an opened controversy spelled out in [12]. This controversy can be reduced to the question whether the heat kernel coefficient a 1 2 is zero or not (this coefficient does in the considered examples not depend on the background potential). Meanwhile more examples were calculated, a number of onedimensional configurations with background potential, non-periodic in [13] and periodic in [14]; again showing different signs of entropy.
The above-mentioned work served as another motivation to consider different thermodynamics properties of the mentioned systems, namely BEC, having in mind a possible correlation. For that we consider periodic lattices of delta function potentials in both one and three dimensions. The questions to answer are whether in one dimension BEC may appear and how BEC will be modified in three dimensions.
In the present paper, in the next section, we display the basic thermodynamic formulas required to formulate a criterion for the occurrence of BEC. In the third section, we introduce some one-dimensional lattices consisting of delta function, generalized delta functions without and with additional potentials. In section 4 we formulate general conditions the spectral functions must obey for BEC and we investigate whether the considered examples allow for BEC. In section 5 we consider a three-dimensional lattice of delta function and show that such background does not destroy BEC. As example, we calculate the correction to the critical temperature. In section 6 we consider a different approach to BEC using the heat kernel expansion. It was successfully applied to BEC with confining potential, see [15]. We apply it to the above mentioned single, free-standing objects and come to the conclusion that BEC is not relevant in that context in the thermodynamic limit.
Throughout the paper we use units with k B = = c = 1.

Free particles and BEC
In this section, we collect basic thermodynamic formulas for free Bose particles in a homogeneous background. This is a topic represented in many textbooks and we restrict ourselves to the simplest criterion for BEC (and follow the book [16], section 12.5). We focus on the particle number and leave out the discussion of long-range order.
Starting point is the grand canonical ensemble of free Bose particles with its partition function (with β = 1/T ) (in notations of [16] we have q = ln Z G , D = Z G ). For instance, the number of particles is with the fugacity z = e βµ , the chemical potential µ and the one particle energies ε i . In order to investigate BEC, one way is to separate the lowest state (ground state) having ε 0 = 0 (possibly, after a redefinition of the chemical potential). In the remaining sum, one increases the number of particles such that the sum can be substituted by an integral, Here, d is the dimension and ε(k) is the one-particle energy; now a function of the momentum k = | k|. The integral in (3) represents the sum over all states. There is no double-counting of the ground state since it is a zero measure set in the integral. If one takes relativistic dynamics, the one particle energy will be and in case of nonrelativistic dynamics one has Doing the angular integrations, the integral in (3) can be reduced to where is the volume of the unit sphere and Γ denotes Euler's gamma function. With the dynamics (4), the temperature dependence can be scaled out, with the frequently used notation and the thermal wavelength Similar formulas appear for the dynamics (5). The first term in (7) is the occupation of the ground state. We divide by the volume, n in the left side being the number of particles per volume, and we can perform the thermodynamic limit, V → ∞, keeping n. In general, there is a solution of equation (10) for any temperature in the range 0 < z 1 and the two terms give the number of particles in the ground state and in the excited states, correspondingly. Condensation takes place if the first term grows with N. This is possible if in the solution z → 1. Accordingly, the contribution from the second terms must be limited. This is the case if its contribution at z = 1 is finite. As a consequence, condensation sets in when any increase in density goes into the first term, i.e. starting from This gives and the critical temperature is T c = 1/λ c . This way, the criterion for the possibility to have BEC is whether the function g d (z) is finite at z → 1. The source of singularity is the zero in the denominator in the integrand in (8) for ε(k) = 0 at k → 0. With the relativistic dynamics (4), obviously one needs d 1 such that the factor k d−1 from the initially d-dimensional integration measure may compensate the zero in the denominator. In that case, from (8) we note with the Riemann zeta function ζ(d), which has a pole for d = 1, and g d (1) is finite for d > 1.
In case of nonrelativistic dynamics, (5), in place of (6) we have and rescaling k → T/2m results in The remaining discussion is identical to the case of relativistic dynamics and we come to the criterion g d/2 (1) must be finite for BEC. With (13) this is the case for d > 2 and the critical temperature T nr c is In the next two section we will use the above formulas in specific systems.

One-dimensional chain of delta functions and some generalizations
In this section we derive the basic formulas for a periodic background composed of delta functions as well as some simple generalizations. We consider massless phonons. Their wave function φ(t, x) obeys, in general, the equation In the case of one spatial dimension, after Fourier transform in time, and with the potential for the periodic background inserted, the corresponding equation reads where φ(x) is the transformed field and we dropped the argument ω. In (18), α is the coupling, a is the spacing of the lattice and the sum goes over all integers. This is a modification of the well known Kronig-Penney [17] model with Dirac delta functions as potentials (Dirac comb). Details and extensions were discussed recently in [18], where also the relevant literature [19], for instance, is cited. We mention that these extensions cover all self adjoint extensions of the wave operator for a potential with point support. Also, there is a close relation to tight binding models in solid state theory.
Having in mind the generalizations to higher dimensions we provide here a short derivation in terms and notation of our earlier paper [4]. A solution of (18) can be searched for in the form where is the free Green's function In (19), the first term is a plane wave solution of the free wave equation and the sum represents the waves scattered from all the delta functions. The coefficients Φ −1 n−n , which form an infinite dimensional matrix, can be found from inserting the ansatz (19) into the equation (18), Because of the translational invariance in equation (18), it is possible to invert the matrix Φ n−n by Fourier transform. Defining we find which may be inserted into the ansatz (21). Now, the eigenfrequencies ω(k) of the field φ(x) follow from the conditioñ (this is equivalent to vanishing determinant of the above mentioned matrix). Inserting (21) into (22) and accounting for (20) we get the representatioñ In this expression, the sum is convergent with ω > 0.
Starting from here, we put the lattice spacing a = 1. The dependence on a can be restored by substitutions ω → aω, k → ak , α → α/a and Φ →Φ/a. With these, we rewrite the sum in the formΦ We apply the Poisson resummation formula to this expression, where we denoted the new summation index by N. It runs also over all integers. In the last line we carried out the integration over t. The sum in the last line converges. The contributions from n = 0 and N = 0 do not pose any problems in accordance with the circumstance that the one-dimensional delta function in a wave equation like (18) is well defined. This representation can be easily generalized to higher dimension as we will show below.
The sum in equation (27) can be carried out resulting iñ Now, from this representation, the condition Φ (ω, k) = 0 can be rewritten in the form which is the well known frequency condition of the Dirac comb. It is known (see, for example [18]) that this frequency condition can be generalized to any potential with support within one lattice cell. Let t(ω) denote the transmission coefficient of the corresponding scattering problem and δ(ω) = 1 2i ln(t(ω)/t * (ω)) the scattering phase shift. Then this generalization reads (equation (34) in [20]) Below we will consider some specific examples and generalizations.

The generalized comb
In this model the potential in the wave equation has in addition a delta'-contribution, (this β is not the inverse temperature). The matching conditions are [21] and the transmission coefficient is The transmission coefficient for a simple delta potential appears with β = 0.

Double delta function
In this model there are two delta functions in each cell, representing for example two species of scattering centers, with The case of a single delta function can be obtained from here by L = 0 and α 1 = α 2 = α/2.

Sturm-Liouville problem framed by generalized delta functions
Consider a potential V(x) with support x ∈ [− L 2 , L 2 ] with L < 1 and the related Sturm-Liouville equation, Besides u(x), this equation has a second independent solution, v(x), with non-vanishing Wronskian w = uv − u v. The scattering setup with two generalized delta potentials at the edges of the interval takes the form, (37) We shifted the periodic lattice to the interval x ∈ − L 2 , L 2 , for convenience. Imposing at x = ± L 2 the boundary conditions corresponding to the generalized delta function, equation (32). After some calculations we come to the transmission coefficient Here we used the notations u ± ≡ u ± L 2 , v ± ≡ v ± L 2 and similar for the derivatives. The primed parameters belong to the delta function in x = −L/2.
With a Pöschl-Teller potential in equation (36) and , such transmission coefficient, which follows from inserting these solutions into (39), was derived in context of the vacuum energy in [22] (formula without number on p 2234). It is also easy to take a step potential For V 0 = 0 we get the transmission coefficient for two generalized delta functions which was considered earlier, for example in [23] and [24].

On the possibility of BEC in one-dimensional combs
In this section we adopt the thermodynamic formulas, represented in section 2 for a homogeneous background, to a periodic background. The general frequency condition is equation (30). For convenience we introduce a new function, for the right side of (30). Now the 'allowed' frequencies follow from |h(ω)| 1 and the spectrum has a band structure where k is the quasi momentum. We denote the solutions of equation (30) by ω n (k), where n = 1, 2, ... counts the allowed bands. The band boundaries are denoted by ω n . The one particle energies (4) are ε = ω n (k) and the expression (6) for the particle number turns into Since now the temperature dependence cannot be scaled out such easy as before, we define, in analogy with (8), a function (45) such that equation (44) after division by the volume turns into (46) Further, the discussion about a possible condensation in the thermodynamic limit goes the same way as in section 2 and reduces to the question whether g (1) is finite. To answer this question we change the integration for ω and come with the Jacobian to the representation where we adjusted the chemical potential such that at the lowest energy we have ω − ω (u) 1 = 0. Now, for z = 1, we have a simple zero in the denominator and have to look for the behavior of the Jacobian at the lowest energy. It should be mentioned that this Jacobian, which is equivalent to the density of states, is the only difference to the case of a free field, mentioned in section 2. A further property of J is that it is real for 'allowed' frequencies only. As mentioned in [14], this property may be used in numerical calculations by first integrating over the whole frequency axis and taking the real part afterward. Now, for investigating g(1), we consider the Jacobian (47) in more detail. Taking the derivative with respect to ω of equation (30) we get i.e. we could express the Jacobian in terms of the function h(ω). Now we can make an expansion at the beginning of the spectrum (assuming no bound states) and come to From here we conclude that for BEC in the one-dimensional case to be possible, the function h(ω), (43), must have vanishing first and also second order (to compensate the zero of the denominator in (48) for z = 1) derivatives in ω = ω

The simple comb
For a comb of single delta functions, where the function h(ω) is given by the right side of (29), we get We use the direct expansion in powers of ω − ω which demonstrates that the above derived criterion cannot be satisfied. What is more, this Jacobian is singular which makes the condensation even more suppressed as compared to the free case.

Thermodynamic properties of more general combs
In this subsection we consider a more general case of a frequency condition (30). In distinction from the preceding subsection, we use the function h(ω), (43), for the right side. Using equations (49) and (50), as mentioned above, h(ω) must have vanishing first and second order derivatives at ω (u) 1 . At the moment we do not know whether such a condition is in conflict with all other necessary conditions or not. We restrict ourselves to the consideration of some specific examples, mentioned above (which do not fulfill this condition).
First, let us consider the generalized comb, defined in section 2. Inserting the transmission coefficient (33) into (30), we get from (43) In place of investigating this function let us first consider the slightly more general case of where a and b are some constants. We get two conditions, The first equation is the frequency condition and the second is the vanishing of the first derivative. The solution of this system is (56) Comparing (53) with (54), we express the parameters of the generalized comb in terms of a and b, β 2 = a−1 a+1 , α = 4b a+1 . However, as can be seen from (56), a < 1 holds (for all ω (u) 1 ) and no real parameter β exists. This way the generalized Dirac comb will not show BEC.
Second, let us consider the double delta function comb as defined in section 2. Here the function h(ω), following from inserting (35) into (43), reads We consider the condition h(ω) = 1. In this case a solution with ω (u) 1 = 0 exists. It follows from the expansion for ω → 0, (note, we have put the lattice spacing a = 1 such that 0 < L < 1 must hold). From (58) we can have h(0) = 1 if the relation between the coupling parameters holds. We mention that this is not the condition for a threshold state. As a result one of the delta functions in each cell of the lattice is attractive. However, inspection shows that nevertheless the spectrum starts from ω = 0. Now we consider higher terms in the expansion of h(ω). Using (59) and dropping the index on α 2 we get As can be seen, the numerator in front of ω 2 has no real zeros for 0 < L < 1. This way, the Jacobian is not singular in the ground state, but still it does not vanish. This is better than in the case with the generalized comb or with the comb with single delta functions, but still not sufficient to have BEC. Third, we consider a Pöschl-Teller potential framed by two generalized delta functions as defined in section 3.2. The parameters are α 1 , β 1 for the potential in x = −L/2 and α 2 , β 2 for the potential in x = L/2 as well as the width L. We need to check whether the conditions h(0) = 1 and h (0) = 0 can be fulfilled (in all considered examples we have h (0) = 0 for symmetry reasons). Since the resulting system of equations is too complex for a general treatment, we first resolve the condition h(0) = 0 for α 2 , which results in quite simple expressions, and consider h (0) as function of the remaining parameters. In all considered cases we did not find h (0) = 0. As illustration we present figure 1. In the left panel the spectral function h(ω) is shown. For −1 < h(ω) < 1 one has the 'allowed' bands. In the given choice of parameters no bound states appear. In the right panel the second derivative, h (0), is shown as function of the width L. It takes only negative values and generalizes the corresponding expression in equation (60). Variations of the additional parameters do not improve the situation.
Fourth, and finally, we consider a rectangular potential framed by generalized delta functions as also introduced in section 3.3 with transmission coefficient given by (41) and (42). We observe a picture similar to the preceding case. Now we have with the potential's height V 0 an additional parameter and we restrict ourselves again to cases with no bound states. As example we show in figure 2 the dependence of h (0) on V 0 . It is seen that an increase of V 0 makes h (0) even more negative. In figure 3 we demonstrate the dependence on α 1 and see that for both signs the situation does not improve.
This way, in all considered examples we did not find that the conditions for BEC were satisfied. We did not make a complete screening of all possible choices of the parameters but restricted ourselves to a number of examples.

Three-dimensional lattice of delta functions
We turn now to the case of a three-dimensional lattice of delta functions and we have, in place of (18), the wave equation where n is a vector of integers and the delta function is three-dimensional. In this section bold letters denote three-dimensional vectors, x ∈ 3 , n ∈ 3 . As well know, a delta function as potential in the wave equation is ill defined in three dimensions. There are several possibilities to give them a precise meaning. These possibilities were recently discussed and compared in [4]. In general, these are all equivalent. We use in the present paper the regularization method  described in sections II.D and III in [4]. According to that method one makes the same ansatz (19) as in the one-dimensional case, which reads now where G (ε) ω (x) is the regularized free Green's function which for ε = 0 turns into In (62) we made use of the translational invariance of the lattice which allows the unknown coefficients Φ −1 n−n to depend on the difference n − n only. The Green's function (63) is singular at x = 0 in opposite to (20). Because of this one cannot put x = am in (62) in order to determine the matrix Φ −1 n−n . This singularity is a consequence of the ill defined delta function in the wave equation (61). The mentioned method assumes the use of a regularized G (ε) ω (x) in place of (63).
We insert the ansatz (62) into equation (61) and arrive at the equation Defining a Fourier transform bỹ we arrive at the expressioñ where we separated the contribution from n = 0 and the prime on the sum denotes to drop the contribution from n = 0. Now the mentioned method implies a renormalization, with the renormalized coupling α r , after which one needs to remove the regularization by putting ε = 0, which gives after the renormalization a finite result. It must be mentioned that the initial coupling α lost its meaning completely and with α r a new parameter came in (for more details see [4]). We arrive at the representation (this is in fact equation (109) in [4]) with the notation and we dropped the subscript 'r' at the coupling. The function J s (ω, k) is a typical lattice sum.
With equation (68) we return to our task to determine the spectrum. It is given by the solutions of the equatioñ Starting from here we put, again, the lattice spacing a = 1. The representation of the lattice sum given by (69) is not very convenient. We rewrite it the following way. First, we use an integral representation for the free Green's function (63), entering (69) for s = 1, We add and subtract the term with n = 0. For that we need to introduce temporarily a factor t s with s → 0 in the end and get (72) The second integral in the last line can be done explicitly, in the first we apply the Poisson resummation formula, (73) The last integral over t could be done explicitly. We denoted the summation index appearing from the Poisson resummation formula by N. Now the analytic continuation in s to s = 0 can be done. It delivers no pole (because of odd dimension). The first term is explicit, in the second the continuation can be done as follows. First, we separate the contribution from N = 0, (74) Then we rewrite the remaining sum, (75) The first sum is an Epstein zeta function Z(2(1 + s)) (with p = 3) as defined in equation (1.15) in [25]. For s = 1 it is a number which we denote by z E . The second sum in the right side in (75) is convergent for s = 1 (if accounting also for N → −N). We denote it by S * (ω, k), and get from (73), Finally, inserting into (68) we arrive at which is a sufficiently convenient representation for the investigation of the spectrum.
In the equation (70), we have the frequency ω and the quasi momentum k. The latter has components restricted to the interval k i ∈ [0, π] as can be seen from (69). The solutions, ω n (k), exhibit a band structure (the bands being numbered by n). However, no explicit formula like (29) or (30) in the one-dimensional case, is available here. Moreover, rotational invariance is broken by the lattice and in place of representation like (7) we are left with which is a generalization of (45). As shown in section 2, to determine the onset of BEC we have to calculate g (1). It is finite in three dimensions like in the free case, i.e. without lattice. However, calculation of thermodynamic quantities like the critical temperature T c must be performed numerically since no analytic expressions are available. The only what can be done analytically is a perturbative approach for small coupling α, which we will do now.
For the solution of equation (70) we make the ansatz where N 0 is the number of the considered band. We insert this ansatz into the equation (70) and use (78).
First we consider the case N 0 = 0, Expanding for small α gives from which we get For N 0 = 0 we need to separate the corresponding contribution in the sum S * , (76), Now, inserting (81) into (70) gives and we identify which completes the perturbative expansion of the spectrum up to second order.
Inserting the perturbative expansion (81) into the function g(z), (80), we get in leading order for the integration/summation measure n 0 ki<π restoring the case of empty space. This gives, up to obvious factors, the function g 3 (z), (8). In the next order we have Since there is no other dependence on k as through (k + 2πN 0 ) 2 in this order, the same argument as in the leading order and (89) applies. The expansion (90) cannot be inserted directly, rather we have to insert this expansion into (80). Thus we have to expand to first order in α with the heat kernel coefficients a k , and insert it into (95). In each term the sum over k can be carried out resulting in polylogarithms, and we arrive at We use this representation to determine the critical temperature as before from z = 1. Using Li s (1) = ζ(s), the Riemann zeta function, we see the same criterion as in section 2; the dimension d must be larger than two (K(t), (96), implies nonrelativistic dynamics) since we otherwise hit the pole of the zeta function at s = 1.
Using the heat kernel expansion we get for the critical temperature T c the equation This equation can be solved by iteration, resulting in the expansion The leading order gives the known expression since the heat kernel coefficient a 0 is proportional to the volume, a 0 ∼ V . The next order is determined by the next non-zero coefficient. If this is a 1 2 , which is proportional to the surface S, we have a smallness of the next order proportional to S/V. If a 1 2 = 0, the correction is even smaller since a 1 is proportional to the perimeter.
This way, in the thermodynamic limit, the critical temperature is determined by the heat kernel coefficient a 0 , which represents the contribution from empty space. All structures, like the mentioned plasma plane or sphere, may give only corrections to the critical temperature. As a consequence, there is no way to relate BEC to negative entropy observed in some Casimir effect related configurations after throwing away the empty space contribution. However, it is possible interesting to mention that the heat kernel coefficient a 1 2 , around which there is the controversy mentioned in the Introduction, is important for the first correction to the critical temperature.

Conclusions
In the forging sections, we investigated the possibility to have BEC in the background of lattices of delta functions. In the one-dimensional case, where there is no BEC in the free case (empty space, without lattice), we found no condensation for the considered examples. These were a lattice of 'simple' delta functions, a lattice of generalized delta functions (i.e. including derivative of delta function) and a lattice with two different delta functions in each primitive cell of the lattice and potentials framed by generalized delta functions.
Also, we formulated a more general condition for BEC in the one-dimensional case in terms of a function h(ω), (43), which is defined in terms of the scattering coefficient t(ω) and which determines the Jacobian (49), or, equivalently, the density of states. The condition states that h(ω) must have vanishing first and second order derivatives at the beginning of the spectrum. Whether any such system exists remains an opened question. We found only systems with no BEC.
In the three-dimensional case one has BEC in the empty space. We showed that the lattice of delta function does not destroy this feature but changes its thermodynamic characteristics. As example, we considered the critical temperature. Since no general analytic results are possible we used a perturbative approach and considered small coupling. We calculated to first order in α the correction to the critical temperature T c , equation (93), which showed a lowering of T c .
Finally, we considered with the heat kernel expansion a different approach to BEC. We showed that for single, free-standing background, the condensation is dominated by the heat kernel coefficient a 0 which is proportional to the volume and which does not depend on the details of the background. These enter at best in the first correction, resulting from a 1 2 , and are proportional to the ratio S/V where S is the surface of the volume V or, for example, the surface of the plasma sphere in [8] and [10]. Taking all considered cases together, we could not spot any relation to the negative entropy observed in Casimir effect related configurations.
We conclude with a speculation on BEC in one dimension. As we have seen in Sect. 4, the simple comb and the comb with two delta functions in each cell may have different suppression of the condensation depending on whether the Jacobian is singular at the origin or finite. It would be interesting to investigate this suppression in more realistic models for bosons on a lattice like extended Hubbard or Bose-Hubbard models [27]. The first integral is a Riemann zeta function and in the second integral we can directly expand, In the first term, P, we substitute y → εy and integrate by parts, Here, the first term can be expanded directly. In the second term we can expand in the logarithm under the sign of the integral and arrive at