Plasma model and Drude model permittivities in Lifshitz formula

At the physical level of rigour it is shown that there are no substantial theoretical arguments in favour of using either plasma mode permittivity or Drude model permittivity in the Lifshitz formula. The decision in this question rests with the comparison of theoretical calculations with the experiment. In the course of the study the derivation of the fluctuation-dissipation theorem is proposed where it is displayed clear at which reasoning stage and in what way the dissipation is taken into account. In particular it is shown how this theorem works in the case of the system with reversible dynamics, that is when dissipation is absent. Thereby it is proved that explicit assertion according to which this theorem is inapplicable to systems without dissipation is erroneous. The research is based on making use of the rigorous formalism of equilibrium two-time Green functions in statistical physics at finite temperature.


I. INTRODUCTION
The problem of theoretical description of the Casimir forces [1] remains, as before, actual.
To a certain extent that is promoted by accumulation of the experimental data in this field.
At the beginning (Casimir, 1948) these forces were associated with the zero point oscillations of electromagnetic field in vacuum [2]. Later on (Lifshitz, 1955) a more broad physical base of these forces was revealed, namely, its close connection with electromagnetic fluctuations in material media was ascertained. This puts in the order of the day correct account of the material characteristics of the continuous media (first of all, its permittivity and permeability) in calculating the Casimir forces and employment here consistent mathematical methods (quantum field theory in media, fluctuation-dissipation theorem and others). The internal structure of the medium proves also to be important [3].
The substantial difficulties arise in these studies when trying to take into account the dissipative properties of the material. A few years ago it was found out that the experimental data on the Casimir forces are described enough good by the Lifshitz formula with the plasma model permittivity. An attempt to take into account the dissipation by substituting into the Lifshitz formula the Drude model permittivity makes worse the description of the experimental data (see, e.g., [4]). These topics were discussed in detail in many publications, we cite here only some of them [5][6][7][8][9]. There arose also certain difficulties in calculation of the Casimir entropy by making use of the Drude model permittivity [5,6,10,11].
Such a situation turned out to be a surprise because the Drude model, taking account of dissipation, is more realistic in comparison with the plasma model which disregards dissipation. This gave rise to a rather long discussion (see papers cited above and references therein). As a result it became clear that it is needed a thorough examination of the physical and theoretical grounds used in derivation of the Lifshitz formula. It is this problem that is a subject of the present article.
One can distinguish two basic approaches to obtaining the Lifshitz formula. First of all the original approach used by Lifshitz himself should be noted [12,13]. Below we argue that it is based substantially on the Callen-Welton fluctuation-dissipation theorem or simply FDT [14][15][16][17]. Many authors have derived the Lifshitz formula by summing up the natural mode contributions to vacuum energy of the macroscopic electromagnetic filed in material media (the mode-by-mode summation or spectral summation, see, for example, [18] and references therein). This approach originates in the Casimir pioneering paper [19]. Such a summation wittingly assumes the real-valued frequencies. Therefore in this approach the media with dissipation cannot be treated because of their complex-valued permittivity, for example, the Drude model permittivity. Thus we have to concentrate ourselves on analysing the Lifshitz derivation of the formula in question.
In the Lifshitz approach there is an extremely complicated point, namely, the transition to the Matsubara imaginary frequencies ω → iζ m , ζ m = 2πm(k B T / ), m = 0, ±1, ±2, . . . and k B is the Boltzmann constant (the rotation of the integration path in the complex ω plane). Later on Lifshitz with coauthors [26,27] repeated the derivation of the formula in question by making, from the very beginning, use of the very complicated formalism of quantum field theory in the statistical physics [28]. The central objects in this approach are the Matsubara Green functions at finite temperature, D ik (ζ n ; r 1 , r 2 ), defined in terms of the retarded Green functions D R ik (ω, r 1 , r 2 ): D ik (ζ n ; r 1 , r 2 ) = D R ik (iζ n , r 1 , r 2 ) .
(1. 1) Obviously in this approach the rotation of the integration contour is not required, however the equality (1.1) must be proved, in other words the transition to the Matsubara frequencies in the retarded Green functions (1.1) should be substantiated. It must be noted that in deriving the Lifshitz formula in Refs. [26,27] the FDT is used also.
Thus in different derivations of the Lifshitz formula, which are of interest for the aim of the present paper, the employment of the FDT is a substantial step. In Appendix A, we accomplish transition to the imaginary frequencies when deriving the Lifshitz formula with the plasma model permittivity at finite temperature.

HAMILTONIAN SYSTEMS
The fluctuation-dissipation theorem [14,17] connects the quantities characterizing spontaneous equilibrium fluctuations (equilibrium correlators) in the system with the generalized susceptibility specifying the linear response of the system to external action. As noted in the Introduction we derive FDT in two stages. As the first stage the averaged anticommutatorcommutator relation (see below) is proved. We use the formalism of equilibrium two-times Green functions at finite temperature with real time in statistical physics [15,16,29,30].
This formalism is based on the Hamiltonian description of the systems under study.
Let us consider a system in thermodynamic equilibrium state at the temperature T which is described by a statistical operator, or density matrix, ρ: where H(q, p) is the Hamiltonian of the system and β −1 = k B T . The averaging with respect to the equilibrium state will be denoted by braces Let A j (q, p) is a classical dynamical variable and A j (t) is its Hermitian operator in the Heisenberg representation We remind the derivation of the so called averaged anticommutator-commutator relation (AC-relation from now onwards) which expresses the symmetrized correlation function of via the averaged values of the retarded commutator of these operators, i.e., in terms of the retarded Green function. For the sake of completeness we define at once the retarded and advanced Green functions In this formula θ(t) is a step function, θ(t) = 1 for t > 0 and θ(t) = 0 at t < 0; the upper sings (+) apply to the retarded Green function, G r (t), and the lower signs, (−), concern the advanced Green function, G a (t). Below it will be shown that the Green functions (2.5), as well as the pair correlators, depend on the deference t i − t j = t. From definition (2.5) we obtain in a straightforward way the following relations between the Green functions and their Fourier transforms In addition it is easy to show that the Green functions (2.5) for the Hermitian operators Obviously, C µ do not depend on t. By means of Eqs. (2.1)-(2.3) and (2.10) we get the spectral representation for the correlator Thus the correlator function A i (t i )A j (t j ) and, consequently, the Green function (2.5) depend on the deference t i − t j = t. Passing on to the Fourier transform we deduce from (2.11) and (2.12) Doing in the same way we obtain the spectral representation for the correlator with the transposed operators A j (t j )A i (t i ) : The last equality (2.14b) is derived by interchanging the summation indexes ν ↔ µ in (2.14a).
We define the Fourier transform of A j (t j )A i (t i ) just as in Eq. (2.12) In view of Eqs. (2.15) and (2.14) we deduce the spectral density J ji (ω) Taking into account the δ-function in Eq. (2.16) we can obviously do the substitution With the help of (2.17) it is easy to derive the spectral density for the symmetrized correlator (2.4) and for the commutator The spectral density J ij (ω) (2.13) for the Hermitian operators A = A † possesses the which can be proved just as the equality (2.17). We shall not represent here this trivial proof.
Now we turn to the construction of the spectral density for the Green functions (2.5).
For this it will be necessary the integral representation of the step function (2.5) (see, for example, [29][30][31]): In this formula, as well as in (2.5), the upper sign (−) corresponds to q = r and the lower sign, (+), corresponds to q = a.
Previously we have noted that the Green functions for the Hermitian operators A = A † are real functions of time (see Eq. (2.9)), therefore their Fourier transforms obey the relation For our consideration the case is of a special interest when The point is the AC-relation can be derived in this case easily. Indeed, utilizing Eq. (2.13) one can show the spectral density J AA † (ω) to be positive By means of the known symbolic equality [31] 1 The spectral density J AA † (ω) is real (see Eq. (2.26)), therefore it follows from Eq. (2.28) Both signs in the right hand side of Eq. (2.30) lead to the same result because from Eq.
Therefore we have in the general case Further we shall use the retarded Green function (q = r) to conform our consideration to the general principle of causality, namely, cause precedes action.

Now we invert definition (2.19)
and substitute the spectral density J {AA † } (ω) from (2.31) into (2.32) with t i = t j = τ . As a result we get At the last step in this equation we have taken into account that Im G r AA † (ω) is an odd function (see Eq. (2.24)). In the first line of (2.33) averaging is carried out over the system at equilibrium. Hence the result does not depend on time τ and the correlator here can be denoted simply by A 2 .

Equation (2.33) expresses the symmetrized equilibrium correlator of the operators A(τ )
and A † (τ ), Eq. (2.4), in terms of the Fourier transform of the retarded Green function, i.e., through the retarded commutator of these operators. This formula is anticommutatorcommutator relation mentioned before Eq. (2.4) as a preliminary task in deriving FDT. It is important to note that this equation is exact and it is applicable only to the systems that admit the Hamiltonian description. In the course of derivation of this relation the system has been considered which is at the thermodynamical equilibrium.
Theoretically the exact AC-relation (2.33) may be used in both directions, namely, for finding the correlator of equilibrium fluctuations through retarded Green function as well as for obtaining this function via the respective correlator. Here it is assumed certainly that the Hamiltonian of the system under study is known. However the exact calculation of the left hand side and the right hand side in AC-relation (2.33) are the tasks of the same complexity level. Therefore the employment of this relation in such calculations does not give advantage. These difficulties are overcome in the physical FDT.
In order to proceed to derivation of the physical FDT we need the connection of the retarded Green function with the generalized susceptibility. In its turn this requires calculation of the linear response of the system under consideration to external action.

III. THE LINEAR RESPONSE OF THE HAMILTONIAN SYSTEM TO EXTERNAL ACTION
In the preceding Section we have derived the AC-relation Eq. (2.33) considering equilibrium system at the temperature T which is described by the Hamiltonian H and the statistical operator ρ (see Eqs. (2.1) and (2.2)). Now we find the linear response of this system to the external action [15,16,29,30]. To this end we add to the Hamiltonian H the following term where A j are, as before, the dynamical variables (operators) pertaining to the initial system and F j (t) are non-operator real functions describing external forces or fields. It is assumed that the fields F j (t) are adiabatically turned on in the remote past Thus we are considering the system with a complete Hamiltonian In this equation the subscript t by the operator H t points to the explicit dependence of this Schrödinger operator on time due to the external fields. The corresponding statistical operator ρ (t) is determined by the equation with the initial condition where ρ is the equilibrium statistical operator (2.1). All the operators, H, A j , H t , ρ, and ρ (t), are in the Schrödinger representation. In the general case the statistical operator ρ (t) in the Schrödinger representation depends on time due to construction (see, for example, Ref. [32]). Now we fulfil the unitary transformation defined by the formulae The parameter t in transformation (3.6) is chosen to be equal to the time t in the functional dependence F j (t) in (3.1).
Strictly speaking the transformation (3.6) implements transition to the interaction representation for the system described by the Hamiltonian (3.2). However in the preceding Section we have defined the Heisenberg operators A j (t) by the same formulae (see Eq. (2.3)).
In Eq. (3.6) the use of the tilde is required, in fact, only for the statistical operator ρ because this operator depends on t in the Schrödinger representation too. For brevity and when this does not lead to misunderstanding we shall call all the operators, depending on t, the Heisenberg operators and use for them the notations A j (t) ≡ A j (t) on the same footing.
The equation for ρ (t) follows from Eq. (3.4) with allowance for (3.6) Here we have taken into account the sign minus in Eq. (3.1). The initial condition for ρ (t) is obtained from Eq. (3.5) Equation (3.7) and condition (3.8) can be united in a single integral equation For the statistical operator in the Schrödinger representation, ρ (t), we obtain from (3.9) Till now our consideration was exact. Further we suppose that the external action (3.1) is sufficiently weak, so that the method of successive approximation is applicable and we can confine ourselves to the first iteration, i.e. substitute ρ for ρ (τ − t) in the right hand side of (3.10) We define the response of the system to external action in the usual way [15,16,29,30] Here we have taken advantage of the equality A i = A i (0) that follows from (3.6). It is worthy to remind that the braces . . . mean, as before, the averaging with respect to equilibrium state of an isolated system with statistical operator ρ (see Eq. (2.2)). With the help of the θ-function we firstly extend the integration in (3.12) to +∞ and then use the definition of the retarded Green function (2.5) Thus the retarded Green function G r ij (t) completely defines the response of the system with the Hamiltonian H to the external influence (3.1) which is taken into account in the linear approximation. (3.14) In this case the generalized susceptibility α(τ ) is defined in the following way (see, for instance, Ref. [17, §123, Eq. (123,2)]: It is obvious that in our considerationx(t) = ∆A(t). Juxtaposing equalities (3.14) and (3.15) we infer that the retarded Green function taken with the sign minus, −G r AA † (τ ), is the generalized susceptibility for the systems in question, i.e. for the Hamiltonian systems.
In other words this statement holds only for systems with reversible dynamics.

IV. FLUCTUATION-DISSIPATION THEOREM AND ITS PRACTICAL USE
In previous Section it was ascertained the relationship of the retarded Green function entering the AC-relation (2.33) with the generalized susceptibility α(τ ) defined in Eq. (3.15).
Proceeding from this we can modify the AC-relation in such a way that its field of applicability extends beyond the Hamiltonian systems including also slightly irreversible systems.
For that we simply replace the retarded Green function −G r AA † (ω) in the AC-relation (2.33) by the generalized susceptibility α(ω) a substitution is permissible when the generalized susceptibility α(ω) is calculated for the system with a reversible dynamics (i.e., for the Hamiltonian system or for the system with slightly irreversible dynamics that is for slightly dissipative system. Herewith in the latter case it is possible to anticipate that the physical FDT Eq. (4.2) will hold with the accuracy proportional to the extent of irreversibility or dissipation. In the next Section V it will be demonstrated by making use of the Drude model permittivity.
Calculation of α(ω) for slightly dissipative system and its substitution into (2.33) for G r AA † (ω) results in appearance in the FDT (4.2) of the phenomenological parameters responsible for dissipation. This is to be noted in view of the following reason.
Practically at once after appearance of the Callen-Welton paper the attempts have begun with the aim to answer the question: whether the FDT takes into account the dissipation and if yes then in what way. Here are a few typical papers on this subject [33]. As far as we know a complete clearness in this problem is absent till now.
All this is in conformity with the results of a consistent application of the exact ACrelation (2.33) when the correlator or the Green function are obtained by making use of the solutions to the pertinent equations generated by the initial Hamiltonian H [29,30]. These equations are infinite sequences of coupled equations containing the correlators and the Green functions of arbitrary high order. The lack of regular procedure for decoupling theses equations or cut-off them does not, in particular, allow one to evaluate the accuracy of the obtained results [34]. In this approach nobody has succeeded in convincing demonstration of a consistent description of dissipation effects starting from the Hermitian Hamiltonians (see the Table at the end of Ref. [29]). Of course, in this approach it is impossible to calculate the parameters describing dissipation. This assertion becomes evident if we remind the fact that the retarded Green function is, up to sign, the generalized susceptibility, which describes the linear response of the system to external action (see Sec. III). A notable advantage of the approach in question is that it provides an opportunity to avoid the explicit quantization of electromagnetic field in a medium, i.e., there is no need to settle the Hamiltonian, to postulate the canonical commutation relations and so on [23,25].
The weak features of this approach are to be noted too. It is obvious that one can judge its applicability to a given problem only after comparison of the calculations with the experiment. Further in the Lifshitz approach the photon Green function in medium at finite temperature turns out to be independent of the temperature. 1 In this approach the temperature dependence emerges only due to the fluctuation-dissipation theorem. Now we ascertain whether it is possible to use the plasma model and Drude model permittivities in the exact AC-relation (2.33) and in the physical FDT (4.2). These permittivities appear in the stated equations via the classical Green functions of the macroscopic Maxwell equations. As it was shown in the Introduction the inferences obtained here will be true for application of the considered permittivities in the Lifshitz formula also.
In the Lifshitz treatment of the electromagnetic field in a medium the FDT is used in the field version (see Eq. (76.6) in Ref. [27]). For the sake of simplicity and for a direct connection with our calculations in preceding Sections we consider, in place of quantum-field retarded Green function the quantum-mechanical analog of (5.1) omitting the dependence on k: G r (ω) = g ω 2 ε(ω) The substitution of the dielectric permittivity calculated in the plasma model for ε(ω) in Eq. (5.2) gives the respective Green function In order that to get in Eq.  Having obtained Im G r pl (ω) we can get, utilizing Eq. (2.29) with q = r, the spectral density of the correlator J AA † (ω): Thus we have here an obvious example of a correct application of FDT in the case when the dissipation is absent (ε pl (ω) in (5.3) is a real function). So the name itself of equality (4.2), the fluctuation-dissipation theorem, is rather conditional. Thereby it is proved also that explicit assertion according to which this theorem is inapplicable to systems without dissipation is erroneous [35].
Exactly to elucidate this point we have divided the derivation of FDT in two steps: at first the AC-relation (2.33) was obtained and after that the FDT (4.2).
Here the following note should be taken also. In our paper [18] we have rigorously derived the Lifshitz formula at zero temperature by making use of the plasma model permittivity.
In the present article, in the Appendix A, we do the same at finite temperature by utilizing the integration contours such as in [18] without using FDT. By the way, as far as we know Lifshitz himself, when calculating the Casimir force, has employed only real dielectric permittivity ε(ω).

Now we turn to the Drude model permittivity
For the sake of formula simplification in what follows we have denoted the relaxation parameter in (5.7) by 2γ.
The substitution of (5.7) in (5.2) gives the retarded Green function with the Drude model Decomposition (5.8) is made in such a way that when γ → 0 the first term, G 1 (ω), turns into G r pl (ω) and the second term vanishes G 2 (ω) → 0. The analytical properties of the Green functions G r D (ω) and G r pl (ω) prove to be substantially different. According to (5.8)-(5.11) the Green function G r D (ω) has three poles (Fig. 1,  right panel). Two of them are located in the lower half-plane ω = ± ω pl − iγ, and when γ → 0 these poles pass into the poles of the Green function G r pl (ω) (Fig. 1, left panel). Evidently for these poles there is no need to specify the bypass rules when integrating with respect to ω because, at positive non-zero values of γ, these poles are located in the lower half-plane of the complex variable ω, certainly outside of the real axis. Hence the symbolic equality (2.27) is not applicable for taking into account the contribution of these poles to the imaginary part of the Green function G r D (ω). In addition this Green function has a 'superfluous' pole at the origin [the first term in the right hand side of (5.10)] which has no prototype in the Green function G r pl (ω). In this case the equality (2.27) works however its application gives contribution to the real part of the Green function G r D (ω) but not to its imaginary part because of the imaginary coefficient (−2iγg/ω 2 pl ) in the front of the pole factor (ω + iε) −1 in the term under consideration. It is to be noted that this irregular finite contribution emerges for arbitrary, no matter how small, values of the constant γ. Hence this contribution is non-analytic in γ when γ → 0. Such an nonanalyticity has been revealed previously in Ref. [11].
This analysis of analytical properties of the Green function G r D (ω) shows that the test of the important condition (5.6) is impracticable in this case. and ω pl . In our consideration the typical values of these parameters are provided by gold [36] 2γ = 34.5 meV, ω pl = 9.03 eV , (5.12) and the aforesaid ratio is 2γ ω pl = 3.83 × 10 −4 . Thus we infer that both permittivities under consideration can be used in the Lifshitz formula practically on the same physical grounds. Therefore it is for comparison with the experiment to decide which permittivity is more suitable to apply in this formula.
It is worthy to note that this inference does not imply that the plasma model and the Therefore it is clear from the physical point of view that the limits of applicability of these two models should be. For example, it is well known that the plasma model permittivity cannot be applicable at low frequencies because it is derived in the limit of high frequencies [37, Chap. IX, § 78]. Theoretically it is reasonable to suggest that the Drude model does not work within some frequency range where it was not independently tested yet.
Besides these two rather simple models it is also interesting to use more complicated models for dielectric permittivity. In Ref. [4] the spatially nonlocal dielectric model is considered which takes dissipation into account and simultaneously agrees with all the measurement data for the Casimir force.
In addition to the Drude model permittivity it is appropriate here to consider briefly the oscillator dielectric permittivity describing isotropic dielectric materials where g j are the oscillator strengths, ω j = 0 are the oscillator frequencies, and γ j are the dissipation parameters. In this model the dielectric atoms are considered to be weakly damped harmonic oscillators (see, for example, [39], [ In the course of study it is explained how the fluctuation-dissipation theorem works in the case of the system with reversible dynamics, that is when dissipation is absent. Thereby it is proved that explicit assertion according to which this theorem is inapplicable to systems without dissipation is erroneous.
The research is based on making use of the rigorous formalism of equilibrium two-time Green functions in statistical physics at finite temperature and real time. This formalism is grounded on the Hamiltonian approach and hence can be considered as a first principle.
In this connection it is also worth noting the recent works where the dissipation properties in graphene are described on the basis of the first principles of QED at non-zero temperature in the framework of the Dirac model [41]. The respective results were used to calculate the Casimir force in the framework of the Lifshitz theory and found to be in excellent agreement with the measurement data (please see [42,43]) and with the requirements of thermodynamics [44].
Appendix A: Transition to the imaginary frequencies in derivation of the Lifshitz formula with plasma model permittivity at finite temperature The imaginary frequencies in the Lifshitz formula at finite temperature T = 0 can be introduced in the same way as it has been done in our work [18] for T = 0. We explain here the main points of this transition for the plasma model at first.
The Casimir energy at T = 0 was defined by us as the spectral sum 3 of the zero point energy of electromagnetic field excitations ω/2: This energy is per unite area of the boundaries. In Eq. (A1) we use the notations which are generally accepted in the problem in question (see Sec. III of our paper [18]). The first term in square brackets in Eq. (A1) is the contribution of the discrete part of the spectrum and the second term represents the contribution of the continuous part of the electromagnetic excitations spectrum in the problem under consideration. In Eq. (A1) and further σ = TE,TM.
The discrete eigenvalues ω n in Eq. (A1) are the frequencies of the surface modes ω sm and the waveguide modes ω wg . They are given by the real positive roots of the frequency where r σ (ω) is the reflection amplitudes [18] Just as in Ref. [18], we assume that ε 2 = µ 1 = µ 2 = 1 and the permittivity ε 1 (ω) is defined in Eq. (5.3). The wave vectors k 1 and k 2 are specified below.
The spectral density shift in Eq. (A1), ∆ρ(ω), as the function of the real frequency ω, is given by a known expression where δ(ω) is the phase shift. Thus Eq. (A1) involves only real positive frequencies ω.
In order to pass to imaginary frequencies the outcome to complex plane ω is needed.
For this aim in Ref. [18] two cuts have been made on this plane, namely, the first cut between the points −ω − (k) and ω − (k) and the second one on the interval (ω + (k), ∞) [see Fig. 2].
The points ω − (k) and ω + (k) are the branch points of the functions and, at the same time, they determine the boundaries between the different branches of the spectrum, namely, on the interval 0 < ω sm (k) < ω − (k) the surface mode frequencies are located, on the interval ω − (k) < ω wg (k) < ω + (k) the waveguide mode frequencies lie, and the interval ω + (k) < ω(k) < ∞ belongs to the continuous spectrum. iii) Along the imaginary positive semi-axis the function D σ (ω) assumes real values.
iv) Equation (A2) has no complex roots in Ω, that is the frequencies of the quasi-normal modes, if they exist, have negative imaginary part. Otherwise the macroscopic electrodynamics would admit the natural modes increasing in time without limits. Recall that the time dependence is taken, as usual, in the form e −iωt . These properties of the quasi-normal modes are explicitly ascertained for electromagnetic oscillations of a perfectly conducting sphere [47] and a dielectric ball without dispersion [48][49][50].
v) As follows from Eq. (A2), the poles of the function D σ (ω) may be generated only by the reflection amplitude, r σ (ω), and their position on the complex ω-plane is independent of the gap width 2a. Therefore the contribution of these poles to the vacuum energy (A1) will be canceled when subtracting in this formula the respective limiting expression obtained for a → ∞. Hence these poles can be disregarded.
vi) The function D σ (ω) possesses the same properties in Ω .
Employing the argument principle theorem (see Sec. 3.4 of Ref. [51]) we represent the contribution of the discrete spectrum into (A1) in the following form where C ε is a circle of the radius ε with ε → 0, see Fig. 2. The spectral density shift ∆ρ, as the function of the complex variable ω, is defined ultimately by the function D σ (ω) (see details in Ref. [18]). In Eq. (A8) we explicitly note that the function D(ω + i0) is calculated at the upper edge of the cut (ω + , ∞) and the function D(ω + i0) is calculated at the lower edge of this cut.
Taking into account Eqs. (A7) and (A8), we can now represent the vacuum energy (A1) in the following form The properties of the functions D σ (ω) and D σ (ω) enumerated above enable us to write the following equalities where the contours C + and C − enclose, respectively, the first and the fourth quadrants of the ω-plane as depicted in Fig. 2.
Further we proceed in the following way: the integral in Eq. (A9), containing the function D σ (ω + i0) on the upper edge of the cut (ω + , ∞), we express from Eq. (A10) and the integral including the function D σ (ω + i0), evaluated on the lower edge of this cut, we express from Eq. (A11). Doing in this way we can disregard, in virtue of the property i), the contributions due to the arcs of the big circle. One can easily see from Fig. 2 that along the contours C ε the contributions of the discrete spectrum and the continuous spectrum are mutually canceled. In addition, the integrals with the function D σ (ω + i0) and D σ (ω + i0) between the origin and the point ω + are reciprocally canceled too. As the result the integration only along the imaginary axis (ω = iζ) survives:  73) and (75) in Ref. [18]).
We shall not write out these well known formulas but at once go over to derivation of the Lifshitz formula at finite temperature T . To this end the free energy F of electromagnetic field in the problem at hand should be found. We accomplish this task by summing up, with respect to the spectrum of electromagnetic excitations, the free energy of quantum oscillator The introduction of respective cuts in the ω-plane can be avoided by using formal and nevertheless rather rigorous method [52,53], namely, the function f (ω) in Eq. (A14) should be represented by the series 4 It is easy to see that after substitution of ω/2 for f (ω) in Eq. (A12) the expression will be here in place of ζ. In this notation Eq. (A12) assumes the form Taking into account the property iii), we can again join two integrals in Eq. (A18) into one The integration by parts with allowance for the property i) gives sin m ζk B T ln D σ (iζ, k) .

(A20)
By virtue of the property ii) the second sum over m in Eq. (A20) does not contribute to the integral over dζ. Now we take advantage of the Fourier series representation for the "comb" The substitution of (A21) into (A20) replaces the integration over dζ by the summation over the Matsubara frequencies (A15) The primed sign of the sum implies that the term with m = 0 should be multiplied by 1/2.