aa r X i v : . [ qu a n t - ph ] O c t Drude model and Lifshitz formula
M. Bordag ∗ Universit¨at Leipzig, Institute for Theoretical PhysicsBox 100 920, 04009 Leipzig, Germany
Abstract
Since nearly 10 years, it is known that inserting the permittivity ofthe Drude model into the Lifshitz formula for free energy causes a viola-tion of the third law of thermodynamics. In this paper we show that thestandard Matsubara formulation for free energy contains a contributionthat is non-perturbative in the relaxation parameter. We argue that thecorrect formula must have a perturbative expansion and conclude that thestandard Matsubara formulation with the permittivity of the Drude modelinserted is not correct. We trace the non-perturbative contribution in thecomplex frequency plane, where it shows up as a self-intersection or a bi-furcation of the integration path.PACS 03.70.+k,11.10.Wx, 65.40.gd
Since nearly 10 years, we observe a problem with thermodynamics in the theoryof the Casimir effect. This problem appears when calculating the free energy F of the electromagnetic field using the Lifshitz formula in the presence of metallicbodies described by the Drude model. The same happens for dielectric bodieswith nonzero dc conductivity. According to the third law of thermodynamics(Nernst heat theorem), one must expect the entropy S = − ∂ F ∂T (1)to vanish for T →
0. In fact, with the Drude permittivity ǫ D = 1 − Ω ω ( ω + iγ ) (2) ∗ [email protected] For γ = 0 it turns into the permittivity of the plasma model. γ ( T ), vanishing for T → T → S = S , (3)which constitutes the violation of the third law of thermodynamics. This wasfound for plane parallel interfaces in [1, 2]. A similar violation was found fordielectrics with finite dc conductivity in [3]. Recently, a similar finding was madefor a sphere in front of a plane [4].It must be mentioned that in a number of papers, [5, 6, 7, 8], nonzero S isnot considered as a violation of thermodynamics or it is attributed to a perfectlysymmetric crystal lattice with no impurities, which is not realized in nature [9].However, these discussions are not the topic of the present paper.This violation constitutes not only a serious theoretical problem. Calculat-ing the free energy as described above, anomalous large thermal contributions atshort separation appear first predicted in [10]. In a very careful short-separationexperiment with metal test bodies [11, 12] these have been shown to be excludedon a high confidence level. Similar large thermal corrections arising from account-ing dc conductivity of dielectrics were also experimentally excluded [13, 14]. In alarge-separation experiment with centimeter-sized metal bodies [15] the predic-tions following from the use of the Drude permittivity were not supported. Inline with this, it must be mentioned that in another large-separation experiment[16] the opposite result, supporting thermal corrections predicted by the Drudepermittivity, was claimed .In the present paper we reconsider the foundation of the theory of the Casimireffect at finite temperature and we show in which place the standard approachwith the Drude permittivity fails. We start from reminding the basic formulasof this approach. In application to two dielectric half-spaces with permittivity ǫ separated by an empty gap with plane parallel surfaces the free energy can beexpressed by the Lifshitz formula, F = T ∞ X l =0 ′ Z d k || (2 π ) X i =TE , TM ln (cid:0) − r i e − ηL (cid:1) , (4)where we have put ~ = c = k B = 1. The prime on the summation sign meansthat the term with l = 0 has to be multiplied by 1/2. In (4), the l -summationgoes over the Matsubara frequencies, ξ l = 2 πlT, (5)and the sum over i accounts for the two polarizations of the electromagnetic field.The reflection coefficients are r TE = η − p ( ǫ − ξ l + η η + p ( ǫ − ξ l + η , The method used in [16] was questioned in [17]. TM = ǫη − p ( ǫ − ξ l + η ǫη + p ( ǫ − ξ l + η , (6)and the notation η = q ξ l + k || (7)is used, where k || is the momentum in the translational invariant directions inparallel to the interfaces. For ǫ one has to insert an expression according tothe model chosen. If ǫ depends on the frequency ω like (2) does, one needs toinsert ω = iξ l . This follows since the Matsubara representation (4) assumes asummation over the imaginary frequencies.It should be mentioned that the free energy (4) with the Drude model per-mittivity, ǫ D , eq. (2), inserted is real. This follows because on the imaginaryfrequency axis, the permittivity ǫ D is real unlike as on the real frequency axis.Equation (4) is one of a number of various equivalent representations of the Lif-shitz formula. These may differ by a change of variables or by integration byparts. This formula was initially derived more than 50 years ago [18]. At finitetemperature, it uses the common representation of a quantum field theory atfinite temperature after Matsubara. Looking at the history and on the many ap-plications, there seems to be no reason for any doubt in the validity of (4). Theremight be only a question about the permittivity that enters through the reflectioncoefficients. The basic idea is that (4) is the free energy of the electromagneticfield interacting with the dielectric bodies. This interaction is accounted for bythe permittivities in the sense of macroscopic electrodynamics. Now, the permit-tivities can be taken from some theoretical model or from measurements, as well.Especially, the latter are well founded by a huge amount of data. In this way,inserting the Drude model permittivity should give a reasonable result.However, there is the much discussed question on what about the dissipationof energy inherent to the Drude model. This dissipation is due to physical pro-cesses like scattering and transfer of energy from the electromagnetic field intoheat. In general, dissipation makes the considered model incomplete as long asthe return of this heat to the electromagnetic field is not accounted for or a heatbath is used like in the Huttner-Barnett type of approach [19]. In physics, thisreturn happens in a natural way since the heat is radiated back.In general, in the quantum field theory, it is common that the energy or thefree energy acquires an imaginary part if the system is not closed and loosesenergy, by creation of particles for example. So, if not including the returnof energy, one should expect the free energy to have an imaginary part. TheLifshitz formula (4) with permittivity (2), interpreted as or derived from thevacuum energy, accounts only for the degrees of freedom of the electromagneticfield in the presence of the lossy dielectrics. By this setup, there is no return ofenergy included. As a consequence, one should expect the free energy to have animaginary part. It is a merit of the formula (4) to have none.The question on how to include dissipation was discussed intensively in [20, 21](see also the explanation in [22], section 7.3). The main idea is to include an3uxiliary field, which acts like a fluctuating source for the electromagnetic fieldand returns, in this way, the dissipated energy. With some restrictions on theparameters, eq. (4) was confirmed also for the case of dissipation [20, 21]. Thisdiscussion is generally considered as giving the correct result.In our discussion, we go a similar way by dividing the system into two. Thefirst is the electromagnetic field in the presence of the dielectric described by theDrude model. The second is some mechanism that performs the return of thedissipated energy. Then, as an approximation, we consider the first part only.Of course, this procedure is quite unsatisfactory insofar as we do not have anyinformation on how good this approximation is. Nevertheless we will be able togain some information on representation (4). In dividing the system, we will havetwo contributions to the free energy. Now, both will have imaginary parts thatmust compensate in the sum. In this sense, it is natural to have free energy withan imaginary part even in equilibrium if considering the first part alone.In line with this setup, we make a point that the free energy must be pertur-bative in the relaxation parameter γ , at least in the sense that for γ →
0, the freeenergy turns into that of the plasma model. This follows from physical reasons interms of the expectation that, for sufficiently small dissipation, the system shouldbe described by the plasma model permittivity. We argue that in the underlyingphysics, there is no process that could create a nonvanishing contribution forvanishing dissipation. In the following, we will use the notion ’perturbative’ ina generalized meaning by allowing for expansion terms like γ ln γ . The point isthat all contributions beyond the plasma model are assumed to vanish for γ → γ .It is interesting to mention that the representation (4) is nonanalytic in γ . Thiscan be inferred already from the known observation that for T → S of the entropy does not depend on the details of how fast γ goes tozero. In terms of the representation γ ( T ) = γ T α + . . . for T → γ drops out for α > ǫ . As observed in [23], for theTE mode, the limits T → ∞ and ǫ → ∞ (ideal conductor limit) do not commute.Considering the ideal conductor as a zeroth order approximation, the dependenceon 1 /ǫ is non–perturbative. However, this is not a good example. One can arguethat this model is not physical since a dielectric with ( ǫ −
1) not decreasing forhigh frequencies cannot be realized in nature. Indeed, as soon as ( ǫ −
1) decreaseslike, for example, in the plasma model, this problem disappears.In the following section, we focus on the TE mode and describe how to formu-late the free energy in the Drude model in terms of physical frequencies and howto set up a perturbative expansion. In the third section, we discuss the relationto a representation in terms of imaginary frequencies. The conclusions are drawnin the last section.Throughout this paper we put ~ = c = k B = 1.4 Free energy in terms of physical frequencies
We start our discussion of free energy from first principles, not resorting to theMatsubara formalism. In general, the free energy F is given by the Gibbs sum, F = − T ln Tr e − βH , (8)where β = 1 /T is the inverse temperature, H is the Hamilton operator, and thetrace goes over the corresponding Hilbert space. For bosonic excitations withone-particle energies resp. frequencies ω J , where J is a generic index numberingthe excitations, this representation turns into F = T X J ln (cid:18) (cid:18) βω J (cid:19)(cid:19) . (9)Also, it can be rewritten in the form F = 12 X J ω − sJ + T X J ln (cid:0) − e − βω J (cid:1) , ≡ E + ∆ T F , (10)where we separated the vacuum energy, E , which remains for vanishing tempera-ture. In (10), a regularization parameter s is used. We do not touch the questionabout the temperature dependence that comes into E through the permittivity.We are interested only in the second part, ∆ T F , which is the temperature–dependent part of the free energy. Of course, representation (8) is not new, seeeq. (5.29) in [24] for example. This formula is very general and it holds for anybosonic system. In application to our system, we restrict our consideration tothe excitations ω J of the electromagnetic field. For these, we can write ∆ T F ina more specific form,∆ T F = T Z d k || (2 π ) (cid:20) X j ln (cid:0) − e − βω j (cid:1) + Z ∞ dk π ln (cid:0) − e − βω k (cid:1) ∂∂k δ ( k ) (cid:21) , (11)where, like in (4), k || is the momentum parallel to the planes and in passing from(10) to (11), ∆ T F became the energy density per unit area. The sum over j goesover the discrete (for a given k || ) frequencies ω j corresponding to the waveguideand evanescent modes. These are all modes whose wave functions decrease in thedirections z → ±∞ . The integral over k accounts for the photonic modes, themomenta are related by k = q k || + k (12)and δ ( k ) = 12 i ln t ( k ) t ( − k ) (13)5s the scattering phase shift expressed in terms of the transmission coefficient, t ( k ) = 4 k q ( k + q ) e − iqL − ( k − q ) e iqL , (14)of the corresponding one-dimensional scattering problem for the TE polarization.Here, L is the width of the gap. The scattering coefficient for TM polarizationdiffers by q in the parentheses to be multiplied by ǫ . In this setup, k is themomentum of the asymptotic states and q = s (1 − ǫ ) k || + k ǫ (15)is the momentum of the photons in the gap. These modes, having a real momen-tum k , correspond to photons propagating in the whole space and we call themphotonic modes. We denote their frequency by ω k .The discrete frequencies ω j appear for imaginary momenta k = i κ j corre-sponding to the poles of the transmission coefficient. Here, we have k = q k || − κ j (16)and the κ j are solutions of the known equations, κ q = tan qL , (symmetric solution) κ q = − cot qL , (anti-symmetric solution) (17)with q given by q = s (1 − ǫ ) k || − κ ǫ . (18)For real q the solutions are waveguide modes, and for imaginary q , the solutionsare the evanescent waves. The latter appear in the plasma model and for TM po-larization only. The eqs. (17) are for TE polarization. That for TM polarizationfollow by multiplying the q in the denominators by ǫ .Now, we turn to the frequencies of these modes. From the Maxwell equations,the frequencies are solutions of the equation ǫ ω = k , (19)where k , eq.(12), is the absolute value of the three–dimensional spatial momen-tum. It is essential to understand that, for a frequency–dependent permittivity, ǫ ( ω ), eq. (19) takes the form ǫ ( ω ) ω = k and that just this equation determinesthe physical frequencies. 6he frequencies appearing in eq. (11) are just solutions of this equation,whereby ω j are the frequencies of the evanescent and waveguide modes withmomentum k = q k || − κ j and ω k is the frequency of the photonic modes with k given by eq. (12).For the Drude model with permittivity (2), we have to insert ǫ → ǫ ( ω ) =1 − Ω / ( ω ( ω + iγ )) in eq. (19), and this equation can be rewritten in the form ofa third–order polynomial, (cid:0) ω − k (cid:1) ( ω + iγ ) − ω Ω = 0 . (20)It has three roots, ω a ( γ, k ) ( a = 1 , , γ . ω ( γ, k )) < ω ( γ, k )) >
0. For γ = 0, it turns into the positivefrequency of the plasma model, ω (0 , k ) | γ =0 = √ Ω + k . (21)The second root has Im( ω ( γ, k )) < ω ( γ, k )) <
0, and it turns for γ = 0 into the negative frequency of the plasma model. These two roots describepropagating modes that are damped according to dissipation. Their startingpoints are ω , ( γ,
0) = 12 (cid:16) − iγ ± p − γ (cid:17) . The third root ω ( γ, k ) is purely imaginary. It starts in ω ( γ,
0) = 0, and it endsin ω ( γ, ∞ ) = − iγ . This root represents the over-damped mode in the Drudemodel. It has no analogy in the plasma model. We mention that these rootshave analytic expressions, which are, however, not really helpful, and we refrainourselves from displaying them. We mention only the relation X i =1 ω i ( γ, k ) = iγ, (22)which holds for these roots. In application to free energy (11), we consider ω ( γ, k )as the frequency of the physical modes.Now, we return to the point of the perturbative expansion for small γ . Fromeq. (20), it follows that the solutions ω i ( γ, k ) all have power series expansions in γ with a finite radius of convergence. The first few orders of these expansions are ω ( γ, k ) = √ Ω + k − iγ Ω + k + γ + 4 k )Ω (Ω + k ) / + O ( γ ) ,ω ( γ, k ) = −√ Ω + k − iγ Ω + k − γ + 4 k )Ω (Ω + k ) / + O ( γ ) ,ω ( γ, k ) = − iγ k Ω + k + O ( γ ) . (23)These are Taylor series expansions with a finite radius of convergence, i.e., theseseries do converge inside a circle, | γ | ≤ γ , in the complex plane of γ . Radius7 Ω Ω Γ - - Ω i Figure 1: Functions ω i ( γe iα , k ) as functions of γ for α = 0, α = 0 . π , α = π and k = 0 .
3. The solid lines are the real parts, and the dashed lines are the imaginaryparts. γ can be found from the explicit formulas for ω i ( γ, k ). Again, these formulasare quite inconvenient. It is easier and more instructive to consider the real andimaginary parts of ω i ( γe iα , k ) as functions of γ for several values of α (see Fig. 1).For γ .
1, the curves do not intersect. For γ &
2, we observe bifurcations. So, weconclude that singularities in the complex γ -plane are all on a finite separationfrom the origin. Hence, the expansion of each ω i ( γ, k ) in powers of γ has a finiteradius of convergence.From the above formulas, it is seen that the perturbative expansion requiresΩ + k = 0. However, this is just the frequency in the plasma model (21), whichis known to be never zero. We mention that this holds not only for the photonicmodes but also for the waveguide and evanescent modes too. At once, one cansee from eqs. (23) that the coefficients decrease for k → ∞ . This allows theconclusion that the expansion can be made under the sign of the integrals in(11).We continue with mentioning that all functions of ω entering the free energy,including the solutions of eq. (17) and the phase shift (13), have power seriesexpansions in ω . This follows from their explicit form for eq.(13), or from theform of the defining equations (17). Now, if we insert into these the expansion of ω i ( γ, k ), we again obtain a power series expansion with a finite radius of conver-gence. Finally, we insert the expansion obtained into the free energy (11). Sincethe integrations converge for each order of the expansion, we conclude that ∆ T F γ with a finite radius of convergence.The zeroth order of this expansion is the free energy of the plasma model.The first order is proportional to iγ and purely imaginary. The second order isproportional to ( iγ ) , and it is real. Higher orders behave accordingly. So, ourconclusion in this section is that the free energy of the first of the consideredsystems goes to zero for the vanishing relaxation parameter γ . Since we expectthe complete system to have the same property, we can conclude that the secondpart, which we do not calculate, has this property too. In this section, we discuss the question about turning of the momentum inte-gration over k in (11) to the imaginary axis. The aim is to investigate how topass from representation (11) in terms of physical frequencies to the standardMatsubara formulation (4) which is in terms of imaginary frequencies. We startfrom reminding this procedure in the plasma model, where it is possible withoutdifficulty. So, we put γ = 0 for the moment and look on the structure of therepresentation (11). We follow the procedure applied in [25]. The integral over k can be split into two integrals by writing the logarithm in (13) as a differenceof two, δ ( k ) = 12 i (ln t ( k ) − ln t ( − k )) . (24)Now, we turn the integration path in the first contribution upwards, k → i κ ,and in the second contribution downwards, k → − i κ . We mention that thereare no contributions from the large circles and that there are no singularitiesthe paths might cross. The complex k -plane is shown in figure 2. The initialintegration path goes along the real axis. It is shown by the dashed lines for thetwo contributions in eq.(24). On the positive imaginary axis, we have the polesof t ( k ), on the negative imaginary axis, poles from t ( − k ). Further, after thepoles, there come the cuts starting from ± i q κ − Ω − k || . Now, the integralsfrom below the cuts have contributions from the poles that just cancel the con-tributions from the discrete modes, i.e., the first part in the square bracket in(11). After that, the integrals along the cuts remain. Here, the two logarithmsin (24) become equal – however, the other logarithms are now different since thefrequency becomes imaginary, ω = iξ in the contribution turned upwards and ω = − iξ in the other, with ξ = q κ − Ω − k || . (25)For these logarithms, we noteln (cid:0) − e − βω (cid:1) = − i βξ (cid:12)(cid:12)(cid:12)(cid:12) I sin βξ (cid:12)(cid:12)(cid:12)(cid:12) + f s ( κ ) (26)9 e( ω )Im( ω ) Figure 2: Complex k -plane for contour rotation in the beginning of section 3.for ω = i κ and ln (cid:0) − e − βω (cid:1) = − i βξ − iπ + ln (cid:12)(cid:12)(cid:12)(cid:12) i sin βξ (cid:12)(cid:12)(cid:12)(cid:12) − f s ( κ ) (27)for ω = − i κ . Function f s ( κ ) = iπ ∞ X l =1 Θ( κ − κ l ) (28)with κ l = q ξ l + Ω + k || (29)is a stair function appearing from passing the points βξ = πl ( l integer), wherethe sinus changes sign. Obviously, ξ l is the Matsubara frequency (5). Next, weobserve that term ( − iβξ ) in the parenthesis just cancels the vacuum energy E ,which was separated in (10). So, we continue with the complete free energy F .Finally, we integrate by parts and come to F = Z dk || (2 π ) (cid:26) −
12 ln t ( i κ ) − Z ∞ q Ω + k || d κ iπ ∂f s ( κ ) ∂ κ ln t ( i κ ) (cid:27) . (30)The first contribution in the curly brackets is the surface term from the ( − iπ )-contribution. Now, the derivative of the stair function gives the sum of delta10unctions, which allow to carry out the κ -integration and, as a result, the Mat-subara sum occurs, F = T Z dk || (2 π ) ∞ X l =0 ′ ln t ( i κ l ) . (31)This expression is known to be a version of the Matsubara sum representation offree energy in the plasma model. For example, it coincides up to notations witheq. (5.5) in [26]. Also, it coincides with (4) up to a contribution that does notdepend on separation L and up to the sum over polarizations that we did notindicate in this section.In this way, we have seen for the plasma model how to pass from the rep-resentation of the free energy as a sum/integral over real frequencies, i.e, overthe physical excitations, to the Matsubara representation. Now, we consider thequestion on whether this is possible for the Drude model too. We take represen-tation (11) with the frequency ω ( γ, k ) as that of the physical excitations of theelectromagnetic field, which need to be accounted for in the free energy.So, we start from considering solution ω ( γ, k ) defined in the preceding section.When turning k → i κ , we have to turn the complete spatial momentum k = q k || + k → iξ ≡ i q k − κ || (32) ω depends on. Therefore, we first investigate ω ( γ, k ) for k = ξe iα under defor-mation α = 0 . . . π .This task is a bit cumbersome, however, it does not pose any principal diffi-culties. For the understanding of deformation, it is necessary to consider all threesolutions, ω a ( γ, k ) ( a = 1 , , ω -plane as functions of ξ for several values of α . The relaxation parameter is taken to be γ = 0 .
1, and theplasma frequency is put Ω = 1. First, we consider α ≤ . π , which is shownin figure 3. Solutions ω , start from the points ω , ( γ,
0) = 12 (cid:16) − iγ ± p + γ (cid:17) , (33)whereas ω starts from ω ( γ,
0) = 0. These points stay fixed under the deforma-tion of α . For ξ → ∞ , solutions ω , behave like the corresponding frequenciesin the plasma model, i.e., ω , ( γ, k → ∞ ) = ± k + . . . . The third solution takesimaginary values for all ξ while α = 0. It goes on the imaginary axis from 0 to − iγ .Now we increase α . The curves deform and are denoted by numbers 2 to 5for increasing α . It is seen that the first two solutions move towards the corre-sponding parts of the imaginary axis, while the third solution deforms to the rightacquiring a positive real part. At some critical value, α ∗ ∼ . ω and ω touch. Beyond, for α > α ∗ , they detach, but in achanged order. That means that the correspondence between starting points (for11 Ω Ω Ω
12 3 4 51 2345 - - - H Ω L - - H Ω L Figure 3: Parametric plots of functions ω a ( ξe − iα ) in the complex ω -plane as afunction of ξ for deforming the path. Curves 1 to 5 correspond to α/ ( π/
2) =0 , . , . , . , . ω = 0 . − . i represent the first solutions, ω , that joining ω = 0 with ω = − i represent thethird solution, ω , and that in the third quadrant represent the second solution, ω . The value of the relaxation parameter is γ = 1. ξ = 0) and end points (for ξ → ∞ ) changed. Hence, a turning of the integrationpath beyond α ∗ is impossible. So, we have to conclude that within the Drudemodel it is not allowed to turn from the physical frequencies to the imaginaryones. Equivalently, it is not possible to come to the Matsubara representation.It is interesting to see what could be done about this problem. Indeed, it ispossible to turn the contour beyond α ∗ . For this, one needs to include from thebeginning the contribution from the over–damped mode, ω ( γ, k ), in eq. (11).In that case, for α > α ∗ , we have a part of the integration path composed ofthe former ω and ω going from ω ( γ,
0) to ω ( γ, ∞ e iα ). These are the curvesin the first quadrant in figure 4. For α = π the final curve of deformation goesjust along the positive imaginary axis. This is obviously the curve used in thestandard Matsubara representation for the Drude model (4). In addition, wehave another contribution from the path joining the starting point of ω withthe end point of ω . These are the curves in the fourth quadrant. These arecomplex also for α = π and, if accounted for, deliver a complex contribution tothe free energy. In the standard Matsubara representation, this contribution isnot present. In this way, the Matsubara representation accounts for a part of thefrequencies originating from ω ( γ, k ) and from ω ( γ, k ) and disregards anotherpart of these.In addition, we would like to mention that solution ω ( γ, k ) (eddy currents)being inserted into free energy (11) causes a problem with the convergence ofthe integrals for large momenta. First, we mention that the convergence in (11)results from the decrease in factor ln(1 − e − βω ) for ω → ∞ . Now, frequency12
234 1 23 41 2 3 4 - - H Ω L - - H Ω L Figure 4: Parametric plots of functions ω a ( ξe − iα ) in the complex ω -plane as afunction of ξ for deforming the path. Curves 1 to 4 correspond to α/ ( π/
2) =0 . , . , . , . γ = 1. ω ( γ, k ) is purely imaginary. Hence, the logarithm does not decrease. One is leftwith the hope that it oscillates sufficiently fast to ensure convergence. However,this is also not the case since ω ( γ, k ) is bounded, 0 ≤ − iω ( γ, k ) ≤
1, for real k . Finally, one could argue that the convergence in (11) could come from thedecrease in the phase shift. However, the integration over k || remains divergent.In [27], the modes corresponding to solution ω ( γ, ω ) (these are the ’eddy’currents), were included in addition to ω ( γ, ω ) in the sum over modes in (10) or(11). However, this alone is also not sufficient to get a real free energy, resp. thestandard Matsubara representation, what is the aim of [27]. Even for α = π/ ω ( γ, ω ) a part of the path is inthe complex plane (curve 4 in Fig. 3). In addition, in [27], also the modes corre-sponding to the solutions ω ( γ, ω ) are included. In that case, all the imaginarycontributions compensate each other and one comes to a real free energy (alreadyon the level of eq.(10) or (11)). However, we cannot include mode ω ( γ, ω ) be-cause in that case, when setting γ = 0, we would get the plasma model twice.In addition, there is another problem with the modes corresponding to solutions ω ( γ, ω ). Their negative real part would cause the sum over the eigenvalues J in eq. (10) to diverge. Of course, this divergence does not enter the force (itdoes not depend on separation L , as can be seen in eq. (9)); however, it wouldproduce an unusual temperature–dependent divergence.13n the remaining part of this section, we will bring an argument showing thatthe Drude model inserted into the standard Matsubara formula has a contributionthat is non perturbative in γ . For that, we consider a contour rotation like above,but in the opposite direction. We start from the standard Matsubara representa-tion (4) of the free energy in the Drude model. It can be represented by eq.(31),which is the free energy for the plasma model in Matsubara representation, whereone has to insert (2) at imaginary frequency, ǫ = 1 + Ω ξ ( ξ + γ ) . (34)The corresponding frequency condition, which gives the connection with theimaginary momentum perpendicular to the plane, reads − ξ = Ω γ/ξ + k || − κ . (35)In opposite to the third–order equation (20) for ω , this is a second–order equationfor κ , which can be solved explicitly, κ ( ξ ) = s ξ + Ω γ/ξ + k || . (36)Now, we go the way back from representation (31), or from eq. (4), to a repre-sentation in terms of real frequencies. This can be done by going from eq. (31)backwards or by applying the Abel-Plana formula to the sum over l in (31). Thelatter procedure was used, for example, in [4], eqs. (11) and (13). In any case, wehave to turn the integration path for ξ towards its imaginary axis, ξ → ± iω . Letus follow what happens to κ ( ξ ), eq.(36). We show in figure 5 the deformation offunction κ ( ωe − iα ) in the complex κ -plane for α = 0 . . . π . Again, we take values γ = 0 . k || = 1 and Ω = 1. For α = 0, the initial curve is along the real axisfrom 1 to ∞ . For increasing α , the curve goes down into the fourth quadrant.These are curves 1 and 2. Further increasing α , at some critical value of α , thecurve develops a cusp and, further, a self-intersection, see curve 3. Increasing α further (curve 4) until the final value of π/ γ ) to the realand imaginary axes. It is clear that, in this way, one comes to an integration overreal frequencies ω , however, one does not come to the frequencies of the Drudemodel which are the first root of eq. (20). Instead, one comes to a representationwith integration over real ω but complex k = i κ ( − iω ). This reflects the factthat we do not consider the complete free energy.In fact, this representation is not really interesting except for the appearanceof the self-intersection and of the loop since the loop is a contribution which isnon-perturbative in γ . We demonstrate this in figure 6, which displays the finalpath (for α = π/
2) for decreasing values of the relaxation parameter γ . The path14 H Κ L - - - - - H Κ L Figure 5: Parametric plot of function κ ( ξe − iα ), eq. (36), in the complex κ -plane as a function of ξ for deforming the path. Curves 1 to 5 correspond to α/ ( π/
2) = 0 . , . , . , . ,
1, respectively. H Κ L - - - - - - H Κ L Figure 6: Parametric plot of function κ ( − iξ ), eq. (36), in the complex κ -planeas a function of ξ for decreasing values of γ . Curve 1 corresponds to γ = 0 .
1, thecurves 2,3,4,5 to γ = 0 . , . , . , . γ = 0, just the free energy ofthe plasma model. Thereby, the contributions from the real axis dropped out.The contribution from the loop does not disappear for any finite γ , no matterhow small, and it constitutes just the non-perturbative contribution.This non-perturbative contribution can also be calculated analytically usingthe methods of [4]. In fact, in [4], limit T → γ → T fixed) can be treated by the same method.We delegate this calculation to the Appendix. The result is∆ T F Drude ( γ → − ∆ T F plasma = T πa f D (0) , (37)where function f D (0) given by eq. (55) in [4] and by eq. (47). Equation (37) is,of course, just the contribution violating the third law of thermodynamics. In the foregoing sections, we considered free energy using the permittivity of theDrude model. We restricted our consideration to that part of the free energy,which results from the excitations of the electromagnetic field in the presence ofthe dielectric bodies and, in this way, we ignored the return of the dissipatedenergy. In this approach, the free energy, or more exactly, a part of it, is givenby eq. (11) and it has an imaginary part. We have described the physical modesof the electromagnetic field and identified that which enter (11). Next, we haveshown that this free energy has a perturbative expansion for small relaxationparameter γ . Combining with the discussion in the Introduction, where we arguedthat the complete free energy must have such an expansion, we can conclude thatthe second part of the free energy, which we did not consider in this paper, alsomust have a perturbative expansion in γ . This is the main conclusion from section2. In section 3, we discussed the relation between representation (11) of the freeenergy in terms of physical frequencies and the standard representation (Lifshitzformula) (4), which is in terms of imaginary frequencies. We started from remind-ing how this transition goes in the plasma model, where it is possible withoutany difficulties. Then we tried to do the same in the Drude model. Here, weobserved a principal difficulty from the behavior of the integration path underrotation in the complex plane. At some angle of rotation, a kind of bifurcationappears and the path looses its uniqueness. From here, we conclude that it isimpossible to pass to the Matsubara representation, at least by contour rotation.Further, we considered the question what happens if ignoring this problem. Itturns out that in doing so one comes to the standard Matsubara representationif integrating along a part of an integration path belonging to the over-damped16ode in the Drude model and, if, at once, dropping a contribution from somepart of the initial integration path belonging to the physical frequencies.We continued in section 3 with a consideration of the above procedure theother way round. We started from the standard Matsubara representation (4)and tried to rotate the integration path back to physical frequencies. Needlessto say that this works fine in the plasma model. However, in the Drude model,under rotation, the integration path develops a self-intersection, see figure 5. Asa result, the integration path contains a loop in the complex plane, which givesa finite contribution no matter how small the relaxation parameter γ is made.From here, we have two conclusions. First, and not unexpected from the abovediscussion, we conclude that one cannot pass from the Matsubara representationto a representation in terms of physical frequencies. Second, we conclude that theMatsubara representation contains a non-perturbative contribution. Also, this isnot unexpected after the discussion at the end of the Introduction.We mention that the procedure, described in the above paragraph, results inthe know representation of the Lifshitz formula in terms of real frequencies (see,e.g., eq. (12.37) in [24], or eq.(40) in the appendix) for the Drude model too.This is, however, different from the representation in terms of physical frequenciesdiscussed in section 2.Finally, we consider the question whether in the perturbative approach thereis a way to go to the Matsubara representation. In section 2, we have seen thatthe free energy of the first part of our system has a perturbative expansion forsmall γ . Thereby, we have seen that factor Ω + k appears in the denominator.This is seen explicitly in eq. (23). From this, the natural conclusion follows thatperturbative expansion is not possible in a representation where this parameteris small. Fortunately, in a representation in terms of physical frequencies, thisparameter never becomes small. Next, we mention that the plasma model is thezeroth order approximation for the Drude model. In the plasma model, we haveseen in section 3 that, in passing from physical to imaginary frequencies, i.e.,when turning the integration path, this path goes finally along the imaginaryaxis starting from ξ = 0. However, in this point, the perturbative expansionis not possible. From here, we draw the conclusion that in the perturbativeapproach to the Drude model, it is impossible to pass to the standard Matsubaraformulation. This backs the above conclusion that (4) contains a non-perturbativecontribution.It is also possible to localize the non-perturbative contribution to (4) by ob-serving that one can pass from (4) to (11) if dropping the contribution resultingfrom integration over the mentioned loop (see Fig. 6). This contribution is com-plex and it has a finite, nonzero limit for γ →
0. At the end of section 3 we haveshown the corresponding analytic expression. This is just the contribution fromthe loop in the integration path, which causes the trouble with thermodynamics.From the above, our conclusion is that for the Drude model one does not havea Matsubara representation. The commonly used one (4) cannot be correct since17t contains a non-perturbative contribution. With eq. (11), one has a part of thecorrect representation. However, one needs to add the other, still unknown, partthat is responsible for the return of the dissipated energy to the electromagneticfield.We would like to mention two papers, [28] and [29], where it was found thatthe complete system is perturbative in γ in the sense we discussed above. In thesepapers, dissipation is included by coupling of the polarization fields (oscillatorsconstituting the medium) to a heat bath like in the Huttner-Barnett model. Theequations of motion for the quantized heat bath field are solved and inserted intothe Maxwell equations. As a result, dissipation and noise polarization appear.Then, the energy is divided into two parts, W and W , in the notations of [29],corresponding to the electromagnetic field with dissipation ( W ) and the noisepart ( W ) which provides the return of energy. That part of the free energy,which we discuss in this paper, just corresponds to W . Now, in that papers,it is shown that in the final answer one can tend γ to zero and arrives just onthe original system, without dissipation. In this sense, the models considered inthese papers are perturbative in γ . Regrettably, these papers are of restrictedapplicability. In [28] a (1+1)-dimensional model is considered, and in [29] themodel is considered for a homogeneous dielectric medium filling the whole space.In the present paper, a number of questions are left unanswered. The mostimportant one is about the derivation of the Lifshitz formula (4) in case of dissi-pation. In light of the forgoing discussions, we have to expect an inconsistency atsome place in the standard derivation. It is clear this is a quite strong statementthat needs further investigation.Another opened question is about a dielectric with dc conductivity, which isknown to have similar problems with the third law of thermodynamics. Here,we expect a similar problem with the Matsubara representation. The point isthat for dc conductivity, representation (11) is perturbative in conductivity σ as can be checked easily. On the other hand, in the Matsubara representation,the dependence on σ is clearly non perturbative. So, the prediction is that thetransformation from physical to imaginary frequencies has a loop hole somewherelike that we found in this paper for the Drude model.Finally, we return to the discussion in the Introduction of considering the firstpart of the free energy as an approximation to the complete one. We have seenthat this allowed us to get qualitative conclusions about the standard formulation.However, it is not possible to draw quantitative conclusions on how good theapproximation is. This means we can calculate the free energy (11) numerically.We can do that also perturbatively. The zeroth order is the plasma model,the first order in γ is purely imaginary and it describes dissipation. It will becompensated if accounting for the second part of the system describing the returnof energy. The second order in γ is real, and it gives a part of the first correctionbeyond the plasma model to the complete free energy. However, the second partis still unknown. Moreover, there are no a’priori reasons to expect its contribution18o be small as compared to the first one. Hence, the only conclusion we have sofar, is that the correction is second order in γ . Since in the experimental situations γ is small (the smallest of all dimensional parameters), this information might behelpful. Acknowledgement
The author benefited from exchange of ideas by the ESF Research NetworkCASIMIR, especially on the highly inspiring workshop
Observability and theo-retical grounding of thermal Casimir forces , Trondheim, Jan 26-27, 2011.The author acknowledges very fruitful discussions with G. Klimchitskaya, V.Mostepanenko and, especially, with G. Barton.
Appendix
In this appendix we calculate the limit γ → F Drude in the standard Matsubara representation with permittivity ǫ D , (2), of theDrude model inserted. Technically, we follow the calculations in [4], eqs. (4)-(16).Applying the Abel-Plana formula we split the free energy, F Drude = E Drude0 + ∆ T F Drude , (38)into the vacuum energy (this is the Matsubara sum substituted by an integral)and the temperature–dependent part,∆ T F Drude = 14 π Z ∞ dx e βx − i ( ϕ ( ix ) − ϕ ( − ix )) . (39)Function ϕ ( ξ ) is given by ϕ ( ξ ) = Z ∞ dk || k || ln (cid:0) − r e − ηL (cid:1) (40)for the real argument and ϕ ( ± ix ) is its analytic continuation.As compared to [4] we changed the notations slightly. So k in [4] is now k || and q in [4] is now η . The reflection coefficient is given by r = q ω p γ/ξ + η − η q ω p γ/ξ + η + η (41)with η = q ξ + k || . (42)19he important step for the following is to divide the integration region in theanalytic continuation ϕ ( ix ) into two parts. The first is k || ∈ [0 , x ] and the secondis k || ∈ [ x, ∞ ).In the first part, we observe that r , eq. (41), can be expanded in powers of γ . This expansion is in fact in powers (cid:0) γix (cid:1) n . However, since k || ≤ x , there is nosingularity in the integrations at least up to the first order ( n =1). From the zerothorder ( n = 0), we get just the plasma model contribution (note η = i q x − k || here). This can be seen using a procedure similar to that in the beginning ofsection 3 and the difference is only in the variables used.In the second contribution,∆ T F Drude2 = 14 π Z ∞ dx e βx − i ( ϕ ( ix ) − ϕ ( − ix )) . (43)we made in the part of the function ϕ ( ix ), (40), resulting from k || ∈ [ x, ∞ ), achange of variable for η = q k || − x (which is real) and denote the result by ϕ ( ix ). It can be written in the form ϕ ( ix ) = Z ∞ dη η ln (cid:0) − r e − ηL (cid:1) , (44)where now r = q ω p γ/ix + η − η q ω p γ/ix + η + η . (45)The integrals in (43) with (44) inserted do converge. If we put γ = 0 in (44), weget formally zero since ϕ ( ± ix ) | γ =0 is real. In this way, ∆ T F Drude2 , eq. (43), doesnot appear in the plasma model at all.The picture is different if we consider (43) for non-zero γ and perform limit γ →
0. This limit cannot be performed under the sign of the integrals directly.The expansion of the intergrand would start with γ/ix and the x -integrationwould linearly diverge for x →
0. However, following the equation that goes in [4]without number between eqs. (51) and (52), one can substitute the x -integrationfor a ζ -integration by x = γζ . One comes to the representation∆ T F Drude2 = 14 π Z ∞ dζ γe βγζ − i ( ϕ ( iγζ ) − ϕ ( − iγζ )) . (46)From (44) and (45), it is seen that ϕ ( iγζ ) is in fact independent on γ . Now thelimit γ → T F Drude2 ( γ →
0) = T π Z ∞ dζζ i ( ϕ ( iγζ ) − ϕ ( − iγζ )) , ≡ T π L f D (0) , (47)which justifies (37). We mention that the rhs of eq.(47) is, up to notations, justeq. (55) in [4]. 20 eferences [1] V. B. Bezerra, G. L. Klimchitskaya, V. M. Mostepanenko. Phys. Rev. A ,062112 (2002)[2] V. B. Bezerra, G. L. Klimchitskaya, V. M. Mostepanenko, C. Romero. Phys.Rev. A , 022119 (2004)[3] B. Geyer, G. L. Klimchitskaya, V. M. Mostepanenko. Phys. Rev. D ,085009 (2005)[4] M. Bordag, I. G. Pirozhenko. Phys. Rev. D , 125016 (2010)[5] M. Bostrom, B. Sernelius. Physica A , 53 (2004)[6] J. S. Hoye, I. H. Brevik, S. A. Ellingsen, J. B. Aarseth. Phys. Rev. E ,051127 (2007)[7] D. A. R. Dalvit, S. K. Lamoreaux. Phys. Rev. Lett. , 163203 (2008)[8] F. Intravaia, C. Henkel. J. Phys. A: Math. Gen. , 164018 (2008). Specialissue conference QFEXT (Leipzig Sep 2007)[9] K. A. Milton. Amer. J. Phys. , 697 (2011)[10] M. Bostrom, B. Sernelius. Phys. Rev. Lett. , 4757 (2000)[11] R. S. Decca, D. L`opez, E. Fischbach, G. L. Klimchitskaya, D. E. Krause,V. M. Mostepanenko. Phys. Rev. D , 077101 (2007)[12] R. S. Decca, D. Lopez, E. Fischbach, G. L. Klimchitskaya, D. E. Krause,V. M. Mostepanenko. Eur. Phys. J. C , 963 (2007)[13] F. Chen, G. L. Klimchitskaya, V. M. Mostepanenko, U. Mohideen. Opt.Express , 4823 (2007)[14] F. Chen, G. L. Klimchitskaya, V. M. Mostepanenko, U. Mohideen. Phys.Rev. B , 035338 (2007)[15] M. Masuda, M. Sasaki. Phys. Rev. Lett. , 171101 (2009)[16] A. O. Sushkov, W. J. Kim, D. A. R. Dalvit, S. K. Lamoreaux. Nature Phys. , 230 (2011)[17] V. B. Bezerra, G. L. Klimchitskaya, U. Mohideen, V. M. Mostepanenko,C. Romero. Phys. Rev. B , 075417 (2011)[18] E. M. Lifshitz. Zh.Eksp.Teor.Fiz. , 94 (1956). [Sov.Phys.JETP. , 73(1956)] 2119] B. Huttner, S. Barnett. Phys. Rev. A , 4306 (1992)[20] Y. S. Barash. Radiophysics and Quantum Electronics , 836 (1973). [Ra-diofizika, (1973) 1086][21] Y. S. Barash, V. L. Ginzburg. Sov.Phys.Usp. , 305 (1975). [Usp.Fiz.Nauk , 5 (1975)][22] P. Milonni. The Quantum Vacuum. An Introduction to Quantum Electrody-namics. (Academic Press, San Diego, 1994)[23] J. Schwinger, L. DeRaad, Jr., K. Milton. Ann. Phys. , 1 (1978)[24] M. Bordag, G. L. Klimchitskaya, U. Mohideen, V. M. Mostepanenko.
Ad-vances in the Casimir Effect (Oxford University Press, 2009)[25] M. Bordag. J. Phys. A: Math. Gen. , 755 (1995)[26] M. Bordag, U. Mohideen, V. M. Mostepanenko. Phys. Rept. , 1 (2001)[27] F. Intravaia, S. A. Ellingsen, C. Henkel. Phys. Rev. A , 032504 (2010)[28] D. Kupiszewska. Phys. Rev. A , 2286 (1992)[29] F. S. S. Rosa, D. A. R. Dalvit, P. W. Milonni. Phys. Rev. A81