# A Generalization of Electromagnetic Fluctuation-Induced Casimir Energy

Authors: Yi Zheng
Journal: Advances in Condensed Matter Physics
 (2015)

---

## Abstract

Intermolecular forces responsible for adhesion and cohesion can be classified according to their origins; interactions between charges, ions, random dipole—random dipole (Keesom), random dipole—induced dipole (Debye) are due to electrostatic effects; covalent bonding, London dispersion forces between fluctuating dipoles, and Lewis acid-base interactions are due to quantum mechanical effects; pressure and osmotic forces are of entropic origin. Of all these interactions, the London dispersion interaction is universal and exists between all types of atoms as well as macroscopic objects. The dispersion force between macroscopic objects is called Casimir/van der Waals force. It results from alteration of the quantum and thermal fluctuations of the electrodynamic field due to the presence of interfaces and plays a significant role in the interaction between macroscopic objects at micrometer and nanometer length scales. This paper discusses how fluctuational electrodynamics can be used to determine the Casimir energy/pressure between planar multilayer objects. Though it is confirmation of the famous work of Dzyaloshinskii, Lifshitz, and Pitaevskii (DLP), we have solved the problem without having to use methods from quantum field theory that DLP resorted to. Because of this new approach, we have been able to clarify the contributions of propagating and evanescent waves to Casimir energy/pressure in dissipative media.


---

## Body

## 1. Introduction

The phenomena of adhesion and cohesion play an important role in many areas of science of technology; they are responsible for stiction in microelectromechanical (MEMS) devices, leading to their failure; microbial adhesion is responsible for the formation of biofilms, which have beneficial as well as detrimental effects; they contribute to friction and wear between surfaces; other areas include food science, pharmacology, and locomotion and prey capture by some animals. Adhesion and cohesion can be loosely defined as the molecular attraction that holds together surfaces of two different substances or two identical substances, respectively. Intermolecular forces responsible for adhesion and cohesion can be classified according to their origins; some forces are of pure electrostatic origin arising out of interactions between charges, ions, random dipole—random dipole (Keesom), random dipole—and induced dipole (Debye); quantum mechanical effects give rise to covalent bonding, London dispersion forces between fluctuating dipoles, and Lewis acid-base interactions; pressure and osmotic forces are of entropic origin [1]. Of all these interactions, the London dispersion (it is so-named because of its relation to the dispersion of light in the visible and UV portion of the spectrum) interaction is universal and exists between all types of atoms as well as macroscopic objects. The focus of this paper is on the contribution of dispersive interactions to cohesion and adhesion between macroscopic objects.Dispersion forces between macroscopic objects, resulting from alteration of the quantum and thermal fluctuations of the electrodynamic field due to the presence of interfaces, play a significant role in the interactions at micrometer and nanometer length scales. Hamaker was the first to extend the concept of London dispersion forces between two atoms to forces between macroscopic spheres separated by vacuum by pairwise summation of the interaction energy between the atoms that constitute the spheres [2]. Hamaker established that the unretarded Casimir force between two half-spaces separated by a film of thickness l is given by AvdW/6πl3, where AvdW is a constant that is now referred to as the Hamaker constant. Since the assumption of pairwise additivity is valid only in the low density limit, alternate theories are necessary for condensed media. Lifshitz, in his seminal work [3], outlined a method based on Rytov’s theory of fluctuational electrodynamics [4] for computing the Casimir forces between two semi-infinite regions separated by a vacuum gap. It required calculation of the average value of the Maxwell stress tensor in the vacuum gap. Lifshitz theory takes many body effects into account in the continuum limit and expressions for AvdW can be derived in terms of the frequency dependent relative dielectric function, ɛ(ω), of the materials involved. For magnetic materials frequency dependent relative magnetic permeability, μ(ω) of the materials also play a role. The generalization of Lifshitz’ theory to the case when the gap between the half-spaces is filled with any dissipative medium is made surprisingly difficult because of the lack of definition of the electromagnetic stress tensor for arbitrary time-varying fields in dissipative media (see p. 161 of [5] and p. 263-264 of [6] for discussions on this topic). It was eventually solved by Dzyaloshinskii et al. (DLP from now on) [5]. They used the Matsubara diagram technique which was developed for the solution of thermal equilibrium problems in quantum many-body theory [5, 7]. This assumption is most practically reflected in the usage of the so-calledMatsubara frequencies in calculating Casimir energy and pressure. The Matsubara frequencies take on the form ξn=2πnkBT/ħ (n=0,1,2,…), where kB is the Boltzmann constant, 2πħ is Planck’s constant, and all objects are at temperature T. Van Kampen et al. [8] and Parsegian and Ninham [9] circumvented the complications of the DLP method but, in doing so, they had to postulate that the free energy of an electromagnetic mode at frequency ωm is given by kBTlog⁡sinh⁡ħωm/2kBT, even though ωm is, in general, complex for dissipative media. Barash and Ginzburg [6] argued that the above-mentioned form of the free energy is the right one for electromagnetic fields in thermal equilibrium with matter [6, 10]. While many authors have attempted to generalize Lifshitz theory to determe Casimir pressure in dissipative media, they do so by assuming an expression for the electrodynamic stress tensor [11] or by defining a Lagrangian density for the electrodynamic field [12], both of which are debatable for media with dissipation [13].A quick survey of chapters 2 and 3 of [14] should convince the reader that a lot has to be learned, in comparison to what is necessary to understand Lifshitz theory of Casimir energy/pressure in vacuum, before one can truly understand the nuances of DLP’s method to calculate Casimir energy/pressure in dissipative media. Yet the expression for AvdW for Casimir/van der Waals forces between two objects separated by vacuum, obtained by Lifshitz’ method, is strikingly similar to the expression derived by DLP. Given the similarity between the expressions of Casimir forces via the two techniques, the disparity between the relative simplicity of Lifshitz’ method and the complications of DLP’s generalization to dissipative media prompted us to look at this problem more closely. The question we asked ourselves was as follows: is it not possible to obtain the Casimir force between two objects separated by a dissipative medium without having to rely on DLP’s method? We answered this question in the affirmative and obtained a more, in our opinion, transparent method for calculating Casimir energy/pressure in dissipative media [15, 16]. The hallmark of the method we developed was that we restricted all calculations of electromagnetic stress tensor to locations in vacuum, even though the eventual goal was to compute the Casimir energy/pressure in a dissipative medium.The outline of the rest of this paper is as follows: in Section2, the principle of conservation of energy is used to relate the Casimir energy and pressure of a multilayer system of thin films to the Casimir energy/pressure of smaller units that comprise the multilayer system. In Section 2.1, the stress tensor is related to the DGF through the fluctuation-dissipation theorem. A commonly used technique of replacing integration along the real frequency axis (ω-axis) by a summation along the imaginary axis (ξ-axis) is described in Section 2.2, along with a discussion of the pros and cons of both representations. In Section 3, the method developed in Section 2 is applied to the case of a thin dissipative film sandwiched between two half-spaces. The contributions of propagating waves and evanescent waves to the Casimir energy/pressure in a multilayer system of thin films are discussed in Section 4. The main points of this paper are summarized in Section 5.


## 2. General Formulation of Casimir/Van Der Waals Energy and Pressure in Multilayered Media

The concepts of work of adhesion and cohesion are usually explained with the aid of thought experiments involving cleaving of two contiguous half-spacesA and B (to be referred to as AB from now on) into two half-spaces A and B at infinite separation with vacuum between them [1, 17, 18], as shown in Figure 1(a). The work of adhesion of AB, WABvac, is the energy required to separate AB into two half-spaces A and B at infinite separation with vacuum between them. The superscript vac is used to indicate that this is the work done in separating A and B across a vacuum gap. The work of adhesion can be related to the three free energies of the three interfaces AB, AV (AV refers to t half-space of A adjacent to a half-space f vacuum), and VB as follows:(1)UAV+UVB=UAB+WABvac,where UAV, UVB, and UAB are free energies of AV, VB, and AB, respectively. The reverse of the procedure shown in Figure 1(a), that is, starting from AV and VB and bringing them together to form AB, is shown in Figure 1(b). Free energy balance for this process corresponds to reorganization of (1) to the form UAB=UAV+UVB-WABvac, so that UAB can be determined if UAV, UVB, and WABvac are known. While Casimir/van der Waals interactions contribute to the energy of half-spaces such as AV, VB, and AB, what is measurable in an experiment is the force between two objects due to Casimir/van der Waals interactions. The Casimir force between two objects exists because of the pressure due to fluctuational electromagnetic waves in the medium between the two objects. It can be shown that the Casimir pressure in a half-space of homogeneous material, such as A or B, is identically equal to zero. So the Casimir force between half-spaces A and B can only be measured when they are separated by a film of finite thickness between them. The main insight of our method is that WABvac can be calculated exactly by using Rytov’s theory of fluctuational electrodynamics and that from such calculations the Casimir pressure in a dissipative material can be determined.Figure 1
Configurations of multilayer structures. (a)AB split into two half-spaces A and B separated by vacuum. (b) Reverse of process described in (a). (c) A multilayer system with N layers between two semi-infinite regions L and R. (d) Splitting N-layer multilayer system into two subcomponents separated by vacuum. The z-axis is perpendicular to the interfaces.
(a)
(b)
(c)
(d)Let us now extend the thought experiments performed in Figures1(a) and 1(b) to a general multilayer system consisting of N layers between half-spaces L and R, as shown in Figure 1(c). The goal is to express the Casimir free energy of the system in terms of combinations of Casimir free energy of smaller units and the work of adhesion necessary to bring those two smaller units into contact. The Casimir free energy per unit area of the configuration shown in Figure 1(c) is represented by ULRz1,…,zN. Each layer is characterized by not only the thickness zk but also the frequency dependent complex permittivity, εk(ω), and permeability, μk(ω) (both relative to that of vacuum). We use the aforementioned notation for free energy for its efficiency. If one of the semi-infinite media is vacuum, the subscript V is used instead of L or R. ULRz1,…,zN can be written as a combination of three terms: (1) the free energy of the first k layers sandwiched by semi-infinite medium L to the left and vacuum to the right of the kth layer, ULVz1,…,zk, (2) the free energy of the remaining N-k layers sandwiched by semi-infinite medium R to the right and vacuum to the left of the (k+1)th layer, UVRzk+1,…,zN, and (3) the work done in bringing the two systems from infinite separation to a separation δ→0, that is, the work of adhesion. This statement can be written as(2)ULRz1,…,zN=ULVz1,…,zk+UVR(zk+1,…,zN)-limδ→0∫δ∞Tzzavgzvdzv,where Tzzavgzv≡Tzzavgz1,…,zk,zv,zk+1,…,zN is the Casimir pressure in the vacuum region in Figure 1(d) and limδ→0∫δ∞Tzzavgzvdzv=WLRvac is the work of adhesion to create the N layer system from the two subsystems. The partial derivative ∂ULRz1,…,zN/∂zr=pLR(r)z1,…,zN gives the Casimir pressure in the rth layer of the N layer system bounded by L and R. For a thin film bounded by two semi-infinite regions, we drop the superscript (r) and denote the pressure simply as pLR. By differentiating (2) with respect to zr, we obtain the following equation for pLR(r):(3)pLR(r)z1,…,zN=∂ULV∂zrz1,…,zk+∂UVR∂zr(zk+1,…,zN)-∫0∞∂Tzzavg∂zrzvdzv.One of the first two terms on the rhs of (3) is zero, depending on whether 1≤r≤k or k+1≤r≤N. Though Tzzavgzv diverges as zv-3 for zv→0, the quantity ∂Tzzavg/∂zr is finite as zv→0∀1≤r≤N, allowing us to define the partial derivative of the last term in (2) as the integral ∫0∞∂Tzzavg/∂zrdzv (see [15] for details). Tzzavgzv is obtained simply by determining the zz component of the Maxwell stress tensor in the vacuum region. Using the procedure described above, we can write the Casimir free energy of any N layer medium in terms of UVV(z1), UVVz2,…, and UVV(zN) and contributions from terms of the form ∫∞δTzzavgzvdzv, all of which involve calculation of the Maxwell stress tensor in vacuum alone. UVV(z) is nothing but the Casimir free energy to create a thin film of thickness z in free space.


### 2.1. Relation betweenTzzavg and Dyadic Green’s Functions

We rely on Rytov’s theory of fluctuational electrodynamics to determine the value ofTzzavg. Though two objects may be neutral, the charges within them undergo random vibrations due to thermal fluctuations as well as quantum (or zero-point) fluctuations. These random vibrations result in electromagnetic waves and the interactions between the two objects because of these waves are described by Rytov’s theory of fluctuational electrodynamics [4, 19]. Fluctuational electrodynamics can be thought of as a combination of quantum mechanics, statistical physics, and macroscopic electrodynamics. Readers of this journal might be more familiar with the usage of fluctuational electrodynamics to analyze near-field effects in thermal radiative energy transfer between two objects at two different temperatures or radiative energy and momentum transfer from a single object (the other object can be considered at an infinite separation) [20]. To quantify the source for electromagnetic waves in conditions of thermal nonequilibrium, the fluctuation-dissipation (FD) theorem of thesecond kind [21, 22] is used to relate the power spectral density of the fluctuation charge density to the local temperature and the imaginary parts of the relative dielectric permittivity, ɛ(ω), and relative magnetic permeability, μ(ω), of the object. The electric and magnetic fields generated because of any charge distribution can be written as integrals of G¯¯e(r,r′;ω) and G¯¯m(r,r′;ω), the electric and magnetic dyadic Green’s functions (DGFs) for the vector Helmholtz equation [23, 24]. G¯¯e(r,r′;ω) and G¯¯m(r,r′;ω) are electromagnetic duals of each other and are solutions of(4)∇×∇×G¯¯(r,r′;ω)-k2G¯¯(r,r′;ω)=I¯¯δr-r′,where I¯¯ is the identity dyad and r and r′ are the position vectors for observation and source points, respectively. G¯¯e and G¯¯m are obtained by enforcing the continuity of (1) μ(r)(n^×G¯¯e(r,r′)), (2) n^×∇×G¯¯er,r′, (3) ɛ(r)(n^×G¯¯m(r,r′)), and (4) n^×∇×G¯¯mr,r′ on either side of an interface defined by the unit normal vector n^ at point r.Thermal energy transfer between two objects occurs only when there exists a temperature difference between them; that is, thermal energy transfer is anonequilibrium phenomenon. However, a finite Casimir/van der Waals force between two objects can exist even when the two objects are at the same temperature; that is, Casimir forces have a finiteequilibrium contribution. It is this equilibrium contribution that we are concerned with in this paper. When all objects are at the same temperature, we can use FD theorem of thefirst kind to simplify the problem considerably by directly relating the cross-spectral power density of the electric and magnetic field components to the elements of the DGFs as (5a)Epr,ωEq*r,ωs=4ωμoΘω,TImGe,pqr,r;ω,(5b)Hpr,ωHq*r,ωs=4ωεoΘω,TImGm,pqr,r;ω,where p,q=x,y,z, Θ=(ħω/2)coth⁡(ħω/2kBT), and Ge,pq and Gm,pq are the pq component of G¯¯e(r,r′;ω) and G¯¯m(r,r′;ω), respectively [25, 26]. The brackets 〈〉 represent a quantum statistical mechanical averaging process over all possible realizations of the fields. The cross-spectral power density is defined such that the average equal time correlation function is given by (6a)Epr,tEqr,t=∫0∞dω2πEpr,ωEq*r,ωs,(6b)Hpr,tHqr,t=∫0∞dω2πHpr,ωHq*r,ωs.The pq component of the average Maxwell stress tensor in vacuum is given by(7)Tpqavg=ε0Epr,tEqr,t-12Epr,tEpr,t+μ0Hpr,tHqr,t-12Hpr,tHpr,t,where repeated indices imply summation over that index. For the multilayered media shown in Figure 1, only the zz component (z-axis is perpendicular to the plane of the thin films) of the stress tensor is necessary. The zz component of the Maxwell stress tensor in vacuum can be expressed in terms of G¯¯e and G¯¯m as [25, 26](8)Tzzavg=∫0∞ħω2πc2coth⁡ħω2kBTImGωdω,where G(ω)=(1/2)(Ge,zz(r,r;ω)-Ge,xx(r,r;ω)-Ge,yy(r,r;ω))+(1/2)(Gm,zz(r,r;ω)-Gm,xx(r,r;ω) − Gm,yy(r,r;ω)). Tzzavg and Gω are independent of the location, r, at which the stress tensor is calculated, in accordance with the law of conservation of linear momentum. For the multilayer structure in Figure 1(d), G(ω) at any location within the vacuum cavity can be written as G(ω)=G0(ω)+Gsc(ω). G0(ω) is the contribution due to background or source radiation that would have been present even in the absence of any boundaries. Gsc(ω) is the contribution from waves that are scattered by the various interfaces. G0(ω) and Gsc(ω) are given by (9a)G0(ω)=ic22πω2∫0∞dkρkzvkρ∑p=e,h‍1,(9b)Gscω=-ic22πω2∫0∞dkρkzvkρ∑p=e,h‍R~vLpR~vRpei2kzvzv1-R~vLpR~vRpei2kzvzv,where p=e,h refer to the transverse electric and transverse magnetic polarizations, respectively, and kzv=(ω2/c2)εvμv-kρ2. R~vL(p) (R~vR(p)) is the generalized reflection coefficient for unit amplitude p-polarized (p=e,h) waves from the vacuum cavity incident at the interface with any structure to the left (right) of the cavity. See Figure 1(d) for a pictorial description of R~vL(p) and R~vR(p). R~vL(p) is a function of the magnitude of the in-plane wavevector, kρ, the thicknesses z1,z2,…,zk, and properties, ɛ(ω) and μ(ω), of the thin films as well as L. Similarly, R~vR(p) depends on the properties of the films and half-space R to the right of the vacuum cavity. R~vL(p) and R~vR(p) can be determined using the transfer matrix method [27–29] or a recursive algorithm [30], both of which are discussed extensively in literature. Clearly, only Gsc(ω) contributes to the Casimir energy/pressure since G0(ω) is independent of dimensions or properties of the constituent films in the multilayer system of interest.


### 2.2. Real and Imaginary Frequencies

Though (9b) and (8) can be used to calculate Tzzavg, they are not usually because of the rapid oscillations of Gsc(ω) due to propagating waves (0≤kρ≤ω/c). Instead, we use contour integration in the complex frequency plane (ω-ξ plane) to convert the integral in (8) into a more convenient sum over poles of the integrand. Any complex frequency can be identified as ω+iξ. G¯¯e and G¯¯m are analytic in the region ξ≥0 by virtue of being response functions. Since G(ω) is a linear combination of different components of G¯¯e and G¯¯m, it is also analytic in the region ξ≥0. In addition, G(-ω)=G*(ω), where * means complex conjugate. The only poles in the complex frequency plane of the integrand in (8) correspond to the poles of coth⁡(ħω/2kBT). The poles lie on the imaginary (ξ) axis and are given by ξn=2nπkBT/ħ, n=0,1,2,…. We can therefore use contour integration, as was done by Lifshitz [3], to replace the integral over ω along the real positive frequency axis by a summation over Matsubara frequencies on the imaginary frequency axis as(10)Tzzavg=-2kBTc2∑n=0∞‍′ξn2Giξn,where ξn=2πnkBT/ħ, Kn=-2ξn2G(iξn)/c2, K0=-limξ→02ξ2G(r,iξ)/c2, and n=0,1,2,…. The prime (′) next to ∑‍ indicates that the n=0 term is given weight 0.5.In (8), ω appears explicitly in the 1/ω2 term as well as implicitly through the ω-dependence of ɛ, μ, and kz. The transformation of relevant functions as one proceeds from the ω-axis to the ξ-axis is shown in Table 1. ɛ(ω)=ɛ′(ω)+iɛ′′(ω) and μ(ω)=μ′(ω)+iμ′′(ω) are complex numbers, the real and imaginary parts of which are related to each other through the Kramers-Kronig relations [30, 31]. However, on the imaginary axis, ɛ(iξ) and μ(iξ) are always positive numbers greater than unity. Since ɛ and μ satisfy identical properties, the differences between the functions along the ω-axis and the ξ-axis are illustrated by plotting two hypothetical dielectric functions, ε1(ω)=1-(1/ω2-1+i0.2ω) and ε2(ω)=1-(1.5/ω2-1.5+i0.2ω), that obey Kramer-Kronig relations along the ω-axis and along the ξ-axis. Along the ξ-axis, the dielectric functions take the form ε1(iξ)=1+(1/ξ2+1+0.2ξ) and ε2(iξ)=1+(1.5/ξ2+1.5+0.2ξ). The qualitative differences between the functions ε1(ω) and ε2(ω) along the ω-axis and ε1(iξ) and ε2(iξ) along the ξ-axis are apparent. The differences between ε1(ω) and ε2(ω) along the ω-axis are compressed to a small portion of the ξ-axis near ξ→0.Table 1
Common function in real and imaginary space.
Real axis
Imaginary axis
Comments
Frequency
ω
i
ξ
n
ω
∈
(
0
,
∞
), ξn=2πnkBT/ℏ,n=0,1,2,…
Electrical permittivity
ε
(
ω
)
ε
(
i
ξ
n
)
ε
′
′
(
ω
)
≥
0, -∞<ε′(ω)<∞, ε(iξn)≥1∀ξn, εv(ω)=εv(iξn)=1∀ω,ξn
Magnetic permeability
μ
(
ω
)
μ
(
i
ξ
n
)
μ
′
′
(
ω
)
≥
0, -∞<μ′(ω)<∞, μ(iξn)≥1∀ξn, μv(ω)=μv(iξn)=1∀ω,ξn
In-plane wavevector
k
ρ
k
ρ
Phase variation inxy plane is eikρ·rin
Out-of-plane wavevector
k
z
(
ω
)
β
z
(
i
ξ
)
k
z
(
ω
)
=
ω
2
c
2
ε
(
ω
)
μ
(
ω
)
-
k
ρ
2   βz(iξ)=ξ2c2ε(iξ)μ(iξ)+kρ2
Phase variation inz-direction
e
i
k
z
z
e
-
β
z
z
In vacuum,kz(ω)is a real number for PWa and imaginary for EWb. βz(iξ) is always positive ∀ξ
aPropagating wave.
bEvanescent wave.Thez-component of wavevector, kz(ω), is given by kz(ω)=(ω2/c2)ɛ(ω)μ(ω)-kρ2. Along the ξ-axis, kz(iξ) is given by kz(iξ)=-(ξ2/c2)ɛ(iξ)μ(iξ)-kρ2=iβz(iξ) and βz(iξ)=ξ2/c2ɛiξμiξ+kρ2≥(ξ/c)ɛ(iξ)μ(iξ). The phase variation along the z-direction along the ω-axis, eikz(ω)z, is transformed to e-βz(iξ)z along the ξ-axis.Taking the transformation of relevant functions from theω-axis to the ξ-axis, (9b) is transformed to(11)Gsciξn=-c22πξn2∫0∞dkρkρβzv∑p=e,h‍R~vLpR~vRpe-2βzvzv1-R~vLpR~vRpe-2βzvzv.Because of the e-βzv(iξ)zv factor, all waves, irrespective of being propagating or evanescent waves, decay with zv. The oscillations that pose numerical difficulties in the integration of (8) are eliminated, making computations along the ξ-axis preferable to the usual computations along the ω-axis. Despite the computational simplicity, a disadvantage is that while the power spectral density of Tzzavg can be obtained directly from the integrand of (8), to do so from (10) is not straightforward. Possibly more important is that the distinction between propagating waves and evanescent waves, which is clear because kzv is real for propagating waves and imaginary for evanescent waves, is eliminated because all waves decay with z along the ξ-axis. How we can allocate Casimir pressure/energy between propagating and evanescent waves is unclear when we rely on computations along the ξ-axis. In Section 3, we will apply this method to calculating Casimir pressure in a thin film (indicated by m) bounded by two semi-infinite objects, L and R, as shown in Figure 2. In Section 4, we will use (8) and (9b) to investigate the contributions of propagating waves and evanescent waves to the Casimir energy/pressure in planar dissipative media.Figure 2
Method proposed by Zheng and Narayanaswamy [15]. To solve the equilibrium Casimir energy/pressure in the thin film surrounded by two half planes, we proceed from problem P1 to problem P2, at the end of which we arrive at the thin film surrounded by two half planes. Problem P3 needs to be solved to determine UVV(zm) and pVV(zm). The “initial configuration” for all three problems corresponds to zv→∞.


## 2.1. Relation betweenTzzavg and Dyadic Green’s Functions

We rely on Rytov’s theory of fluctuational electrodynamics to determine the value ofTzzavg. Though two objects may be neutral, the charges within them undergo random vibrations due to thermal fluctuations as well as quantum (or zero-point) fluctuations. These random vibrations result in electromagnetic waves and the interactions between the two objects because of these waves are described by Rytov’s theory of fluctuational electrodynamics [4, 19]. Fluctuational electrodynamics can be thought of as a combination of quantum mechanics, statistical physics, and macroscopic electrodynamics. Readers of this journal might be more familiar with the usage of fluctuational electrodynamics to analyze near-field effects in thermal radiative energy transfer between two objects at two different temperatures or radiative energy and momentum transfer from a single object (the other object can be considered at an infinite separation) [20]. To quantify the source for electromagnetic waves in conditions of thermal nonequilibrium, the fluctuation-dissipation (FD) theorem of thesecond kind [21, 22] is used to relate the power spectral density of the fluctuation charge density to the local temperature and the imaginary parts of the relative dielectric permittivity, ɛ(ω), and relative magnetic permeability, μ(ω), of the object. The electric and magnetic fields generated because of any charge distribution can be written as integrals of G¯¯e(r,r′;ω) and G¯¯m(r,r′;ω), the electric and magnetic dyadic Green’s functions (DGFs) for the vector Helmholtz equation [23, 24]. G¯¯e(r,r′;ω) and G¯¯m(r,r′;ω) are electromagnetic duals of each other and are solutions of(4)∇×∇×G¯¯(r,r′;ω)-k2G¯¯(r,r′;ω)=I¯¯δr-r′,where I¯¯ is the identity dyad and r and r′ are the position vectors for observation and source points, respectively. G¯¯e and G¯¯m are obtained by enforcing the continuity of (1) μ(r)(n^×G¯¯e(r,r′)), (2) n^×∇×G¯¯er,r′, (3) ɛ(r)(n^×G¯¯m(r,r′)), and (4) n^×∇×G¯¯mr,r′ on either side of an interface defined by the unit normal vector n^ at point r.Thermal energy transfer between two objects occurs only when there exists a temperature difference between them; that is, thermal energy transfer is anonequilibrium phenomenon. However, a finite Casimir/van der Waals force between two objects can exist even when the two objects are at the same temperature; that is, Casimir forces have a finiteequilibrium contribution. It is this equilibrium contribution that we are concerned with in this paper. When all objects are at the same temperature, we can use FD theorem of thefirst kind to simplify the problem considerably by directly relating the cross-spectral power density of the electric and magnetic field components to the elements of the DGFs as (5a)Epr,ωEq*r,ωs=4ωμoΘω,TImGe,pqr,r;ω,(5b)Hpr,ωHq*r,ωs=4ωεoΘω,TImGm,pqr,r;ω,where p,q=x,y,z, Θ=(ħω/2)coth⁡(ħω/2kBT), and Ge,pq and Gm,pq are the pq component of G¯¯e(r,r′;ω) and G¯¯m(r,r′;ω), respectively [25, 26]. The brackets 〈〉 represent a quantum statistical mechanical averaging process over all possible realizations of the fields. The cross-spectral power density is defined such that the average equal time correlation function is given by (6a)Epr,tEqr,t=∫0∞dω2πEpr,ωEq*r,ωs,(6b)Hpr,tHqr,t=∫0∞dω2πHpr,ωHq*r,ωs.The pq component of the average Maxwell stress tensor in vacuum is given by(7)Tpqavg=ε0Epr,tEqr,t-12Epr,tEpr,t+μ0Hpr,tHqr,t-12Hpr,tHpr,t,where repeated indices imply summation over that index. For the multilayered media shown in Figure 1, only the zz component (z-axis is perpendicular to the plane of the thin films) of the stress tensor is necessary. The zz component of the Maxwell stress tensor in vacuum can be expressed in terms of G¯¯e and G¯¯m as [25, 26](8)Tzzavg=∫0∞ħω2πc2coth⁡ħω2kBTImGωdω,where G(ω)=(1/2)(Ge,zz(r,r;ω)-Ge,xx(r,r;ω)-Ge,yy(r,r;ω))+(1/2)(Gm,zz(r,r;ω)-Gm,xx(r,r;ω) − Gm,yy(r,r;ω)). Tzzavg and Gω are independent of the location, r, at which the stress tensor is calculated, in accordance with the law of conservation of linear momentum. For the multilayer structure in Figure 1(d), G(ω) at any location within the vacuum cavity can be written as G(ω)=G0(ω)+Gsc(ω). G0(ω) is the contribution due to background or source radiation that would have been present even in the absence of any boundaries. Gsc(ω) is the contribution from waves that are scattered by the various interfaces. G0(ω) and Gsc(ω) are given by (9a)G0(ω)=ic22πω2∫0∞dkρkzvkρ∑p=e,h‍1,(9b)Gscω=-ic22πω2∫0∞dkρkzvkρ∑p=e,h‍R~vLpR~vRpei2kzvzv1-R~vLpR~vRpei2kzvzv,where p=e,h refer to the transverse electric and transverse magnetic polarizations, respectively, and kzv=(ω2/c2)εvμv-kρ2. R~vL(p) (R~vR(p)) is the generalized reflection coefficient for unit amplitude p-polarized (p=e,h) waves from the vacuum cavity incident at the interface with any structure to the left (right) of the cavity. See Figure 1(d) for a pictorial description of R~vL(p) and R~vR(p). R~vL(p) is a function of the magnitude of the in-plane wavevector, kρ, the thicknesses z1,z2,…,zk, and properties, ɛ(ω) and μ(ω), of the thin films as well as L. Similarly, R~vR(p) depends on the properties of the films and half-space R to the right of the vacuum cavity. R~vL(p) and R~vR(p) can be determined using the transfer matrix method [27–29] or a recursive algorithm [30], both of which are discussed extensively in literature. Clearly, only Gsc(ω) contributes to the Casimir energy/pressure since G0(ω) is independent of dimensions or properties of the constituent films in the multilayer system of interest.


## 2.2. Real and Imaginary Frequencies

Though (9b) and (8) can be used to calculate Tzzavg, they are not usually because of the rapid oscillations of Gsc(ω) due to propagating waves (0≤kρ≤ω/c). Instead, we use contour integration in the complex frequency plane (ω-ξ plane) to convert the integral in (8) into a more convenient sum over poles of the integrand. Any complex frequency can be identified as ω+iξ. G¯¯e and G¯¯m are analytic in the region ξ≥0 by virtue of being response functions. Since G(ω) is a linear combination of different components of G¯¯e and G¯¯m, it is also analytic in the region ξ≥0. In addition, G(-ω)=G*(ω), where * means complex conjugate. The only poles in the complex frequency plane of the integrand in (8) correspond to the poles of coth⁡(ħω/2kBT). The poles lie on the imaginary (ξ) axis and are given by ξn=2nπkBT/ħ, n=0,1,2,…. We can therefore use contour integration, as was done by Lifshitz [3], to replace the integral over ω along the real positive frequency axis by a summation over Matsubara frequencies on the imaginary frequency axis as(10)Tzzavg=-2kBTc2∑n=0∞‍′ξn2Giξn,where ξn=2πnkBT/ħ, Kn=-2ξn2G(iξn)/c2, K0=-limξ→02ξ2G(r,iξ)/c2, and n=0,1,2,…. The prime (′) next to ∑‍ indicates that the n=0 term is given weight 0.5.In (8), ω appears explicitly in the 1/ω2 term as well as implicitly through the ω-dependence of ɛ, μ, and kz. The transformation of relevant functions as one proceeds from the ω-axis to the ξ-axis is shown in Table 1. ɛ(ω)=ɛ′(ω)+iɛ′′(ω) and μ(ω)=μ′(ω)+iμ′′(ω) are complex numbers, the real and imaginary parts of which are related to each other through the Kramers-Kronig relations [30, 31]. However, on the imaginary axis, ɛ(iξ) and μ(iξ) are always positive numbers greater than unity. Since ɛ and μ satisfy identical properties, the differences between the functions along the ω-axis and the ξ-axis are illustrated by plotting two hypothetical dielectric functions, ε1(ω)=1-(1/ω2-1+i0.2ω) and ε2(ω)=1-(1.5/ω2-1.5+i0.2ω), that obey Kramer-Kronig relations along the ω-axis and along the ξ-axis. Along the ξ-axis, the dielectric functions take the form ε1(iξ)=1+(1/ξ2+1+0.2ξ) and ε2(iξ)=1+(1.5/ξ2+1.5+0.2ξ). The qualitative differences between the functions ε1(ω) and ε2(ω) along the ω-axis and ε1(iξ) and ε2(iξ) along the ξ-axis are apparent. The differences between ε1(ω) and ε2(ω) along the ω-axis are compressed to a small portion of the ξ-axis near ξ→0.Table 1
Common function in real and imaginary space.
Real axis
Imaginary axis
Comments
Frequency
ω
i
ξ
n
ω
∈
(
0
,
∞
), ξn=2πnkBT/ℏ,n=0,1,2,…
Electrical permittivity
ε
(
ω
)
ε
(
i
ξ
n
)
ε
′
′
(
ω
)
≥
0, -∞<ε′(ω)<∞, ε(iξn)≥1∀ξn, εv(ω)=εv(iξn)=1∀ω,ξn
Magnetic permeability
μ
(
ω
)
μ
(
i
ξ
n
)
μ
′
′
(
ω
)
≥
0, -∞<μ′(ω)<∞, μ(iξn)≥1∀ξn, μv(ω)=μv(iξn)=1∀ω,ξn
In-plane wavevector
k
ρ
k
ρ
Phase variation inxy plane is eikρ·rin
Out-of-plane wavevector
k
z
(
ω
)
β
z
(
i
ξ
)
k
z
(
ω
)
=
ω
2
c
2
ε
(
ω
)
μ
(
ω
)
-
k
ρ
2   βz(iξ)=ξ2c2ε(iξ)μ(iξ)+kρ2
Phase variation inz-direction
e
i
k
z
z
e
-
β
z
z
In vacuum,kz(ω)is a real number for PWa and imaginary for EWb. βz(iξ) is always positive ∀ξ
aPropagating wave.
bEvanescent wave.Thez-component of wavevector, kz(ω), is given by kz(ω)=(ω2/c2)ɛ(ω)μ(ω)-kρ2. Along the ξ-axis, kz(iξ) is given by kz(iξ)=-(ξ2/c2)ɛ(iξ)μ(iξ)-kρ2=iβz(iξ) and βz(iξ)=ξ2/c2ɛiξμiξ+kρ2≥(ξ/c)ɛ(iξ)μ(iξ). The phase variation along the z-direction along the ω-axis, eikz(ω)z, is transformed to e-βz(iξ)z along the ξ-axis.Taking the transformation of relevant functions from theω-axis to the ξ-axis, (9b) is transformed to(11)Gsciξn=-c22πξn2∫0∞dkρkρβzv∑p=e,h‍R~vLpR~vRpe-2βzvzv1-R~vLpR~vRpe-2βzvzv.Because of the e-βzv(iξ)zv factor, all waves, irrespective of being propagating or evanescent waves, decay with zv. The oscillations that pose numerical difficulties in the integration of (8) are eliminated, making computations along the ξ-axis preferable to the usual computations along the ω-axis. Despite the computational simplicity, a disadvantage is that while the power spectral density of Tzzavg can be obtained directly from the integrand of (8), to do so from (10) is not straightforward. Possibly more important is that the distinction between propagating waves and evanescent waves, which is clear because kzv is real for propagating waves and imaginary for evanescent waves, is eliminated because all waves decay with z along the ξ-axis. How we can allocate Casimir pressure/energy between propagating and evanescent waves is unclear when we rely on computations along the ξ-axis. In Section 3, we will apply this method to calculating Casimir pressure in a thin film (indicated by m) bounded by two semi-infinite objects, L and R, as shown in Figure 2. In Section 4, we will use (8) and (9b) to investigate the contributions of propagating waves and evanescent waves to the Casimir energy/pressure in planar dissipative media.Figure 2
Method proposed by Zheng and Narayanaswamy [15]. To solve the equilibrium Casimir energy/pressure in the thin film surrounded by two half planes, we proceed from problem P1 to problem P2, at the end of which we arrive at the thin film surrounded by two half planes. Problem P3 needs to be solved to determine UVV(zm) and pVV(zm). The “initial configuration” for all three problems corresponds to zv→∞.


## 3. Casimir/Van Der Waals Pressure in Thin Films

We start with the assertion that the Casimir pressure in any infinite or semi-infinite planar medium is zero. To solve for the Casimir pressure in the dissipative thin filmm, we will perform a free energy balance for the three problems marked P1, P2, and P3 in Figure 2. Each problem has an initial and a final configuration. Equation (8) or (10) can be used to determine the pressure against which work needs to be done and the work of adhesion (or cohesion) required to affect the change from the initial to final configuration in each of the problems. The conservation of energy for each of the problems is given below: (12a)P1:ULVzm=UVVzm-limδ→0∫δ∞Tzzavg,1dzv,(12b)P2:ULRzm=ULVzm-limδ→0∫δ∞Tzzavg,2dzv,(12c)P3:0=UVV(zm)-limδ→0∫δ∞Tzzavg,3dzv.The lhs and rhs of (12a) and (12b) correspond to the free energy of the final and initial configurations of P1 and P2, respectively. In P3, the initial configuration corresponds to a half-space of m and a film of m at infinite separation from each other. The final configuration is a half-space of m. Hence the free energies of the two half-spaces cancel each other, leaving only the free energy of the film, UVV(zm), and the work done in brings the two objects together as the two terms in (12c). In (12a) UVV(zm) as well as ULV(zm) are unknown; in (12b) ULV(zm) and ULR(zm) are unknown; in (12c) UVV(zm) is the only unknown. Tzzavg,1, Tzzavg,2, and Tzzavg,3 can all be calculated using (8). Differentiating (12a), (12b), and (12c) with respect to zm, we obtain the following equations for the Casimir pressure (rearranged in the order in which the three problems should be solved): (13a)P3:pVV(zm)=∫0∞∂Tzzavg,3∂zmdzv,(13b)P1:pLV(zm)=pVV(zm)-∫0∞∂Tzzavg,1∂zmdzv,(13c)P2:pLR(zm)=pLV(zm)-∫0∞∂Tzzavg,2∂zmdzv.Tzzavg,1, Tzzavg,2, Tzzavg,3 can be evaluated using (10) and the details are given in [15]. We give only the final result here for the Casimir pressure in the thin film surrounded by L and R:(14)pLRzm=kBTπc3∑n=0∞‍′εmμm3/2ξn3∫1∞dqq2×∑p=e,h‍RmL(p)RmR(p)e-2qξnεmμmzm/c1-RmL(p)RmR(p)e-2qξnεmμmzm/c,where q=βzm/(ξnεm(iξn)μm(iξn)/c). Equation (14) agrees with the expression for Casimir pressure in a thin film according to DLP [5, 32]. The only complication in extending this method to multilayered media is to determine the appropriate reflection and transmission coefficients [30, 33] and solve more problems instead of just three for the configuration shown in Figure 2. In general, for a multilayer structure with N films between L and R (see Figure 1(c)), we will need to solve 2N+1 problems (N problems of finding work of cohesion similar to P3 and N+1 problems of finding work of adhesion similar to P1 and P2).


## 4. Role of Propagating and Evanescent Waves in Casimir Energy/Pressure in Dissipative Media

The dependence of the work of adhesion,WLRvac=limδ→0∫δ∞Tzzavgdzv, between any two subsystems of a planar multilayer object on the properties and thicknesses of the constituent films arises entirely through Δsc(ω)=limδ→0∫δ∞Gsc(ω)dzv. The implication of splitting G(ω) into G0(ω) and Gsc(ω) is illustrated in Figure 3 by referring to the cavity in problem P1 of Figure 2. The work of adhesion in problem P1 is the work done in translating the film in Figure 3 from zv→0 to zv→∞. In this case, the net force per unit area on the film is obtained by integrating the zz component of the stress tensor on either side of the thin film. Within the cavity, both G0(ω) and Gsc(ω) contribute to the stress tensor. To the right of the thin film, scattered waves, though present, do not contribute to the stress tensor and G(ω)=G0(ω). The contributions to the stress tensor from G0(ω) on either side of the thin film cancel each other, and the net force on the thin film is dependent only on Gsc(ω). Should the thin film be replaced by a half-space, the force per unit are can be determined in two ways: (1) assume that the half-space is a film of finite thickness and take the limit as that thickness approaches infinity or (2) consider only the stress tensor within the cavity, keeping in mind that G0(ω) does not contribute to Casimir energy/pressure because it is not a function of thickness or optical properties of any constituent materials. In either case, the force that is responsible for Casimir energy/pressure arises only from Gsc(ω).Figure 3
Contributions to electromagnetic stress tensor in the vacuum cavity and the vacuum half-space to the left and the right, respectively, of the thin film. Within the cavity, source radiation as well as scattered radiation contribute to stress tensor. To the right of the thin film, only source radiation contributes to the stress tensor.The integral expression forGsc(ω) can be split into GscPW(ω), due to propagating waves in vacuum (0≤kρ<ω/c) and GscEW(ω), due to evanescent waves in vacuum (ω/c≤kρ<∞). The contributions of propagating and evanescent waves are proportional to ΔscPW(ω)=limδ→0∫δ∞ImGscPW(ω)dzv and ΔscEW(ω)=limδ→0∫δ∞ImGscEW(ω)dzv, respectively. For propagating waves kzv is a real number whereas it is an imaginary number given by kzv=ikzv=iKzv for evanescent waves. The expressions for GscPW(ω) and GscEW(ω) are given by(15a)GscPWω=-ic22πω2∫0ω/cdkρkzvkρ∑p=e,h‍R~vLpR~vRpei2kzvzv1-R~vLpR~vRpei2kzvzv,(15b)GscEW(ω)=c22πω2∫ω/c∞dkρKzvkρ∑p=e,h‍R~vLpR~vRpe-2Kzvzv1-R~vLpR~vRpe-2Kzvzv.Since GscEW(ω) is a decaying function of zv, (15b) can be readily integrated with respect to zv to yield the following equation:(16)ΔscEWω=c24πω2limδ→0∫ω/c∞dkρkρ×Im∑p=e,h‍ln⁡1-R~vLpR~vRpe-2Kzvδ.To compute ΔscPW(ω), we use the multiple-reflection expansion of the factor (1-R~vL(p)R~vR(p)ei2kzvzv)-1 to rewrite (15a) as(17)GscPW(ω)=-ic22πω2∫0ω/cdkρkzvkρ∑p=e,h‍∑n=1∞‍R~vLpR~vRpei2kzvzvn.For propagating waves, limδ→0 can be replaced by δ=0 without incurring any error. Since ∫0∞eixdx=0, it can be seen from (17) that ΔscPW(ω)≡0. The implication of the result ΔscPW(ω)≡0 is thatpropagating waves do not contribute to the Casimir energy or pressure in a multilayer system with only dissipative thin films. However, this claim is invalid if any of the thin films in the multilayer system, say, the film k+1 in Figure 1(c), is a vacuum layer. Though the method described until now is applicable to this case too, (2) can be simplified to take into account the presence of the vacuum layer to yield the following:(18)ULRz1,…,zN=ULVz1,…,zk+UVR(zk+2,…,zN)-∫zk+1∞Tzzavgzvdzv.The integration of Tzzavg now proceeds only from zk+1 to ∞, instead of from zv→0 to ∞, because of which there is a finite contribution from propagating waves too.Propagating waves contribute to Casimir energy/pressure only if at least one of the thin films is a vacuum layer. It should be kept in mind that classifying a plane wave as propagating or evanescent is based on the nature of the plane wave in vacuum.


## 5. Summary

Casimir/van der Waals energy and pressure in macroscopic objects arise due to the modifications of zero-point and thermal fluctuations of the electromagnetic field due to the presence of interfaces. Though Lifshitz’ original work on Casimir/van der Waals forces between two half-spaces separated by a vacuum gap relies on Rytov’s theory of fluctuational electrodynamics, the generalization to Casimir energy/pressure in dissipative media by Dzyaloshinskii, Lifshitz, and Pitaevskii was possible only by using techniques from quantum field theory. Possibly because of this high bar imposed by DLP’s method, not many have explored the limits or alternatives to DLP’s method. The condition of thermal equilibrium has to be satisfied for the applicability of DLP’s method. Thermal nonequilibrium contributions to Casimir pressure, which could be important in applications like heat assisted magnetic recording because of the temperature differences between the region heated by the focused laser spot and the rest of the body, cannot be determined using DLP’s method.In a multilayer system composed of planar films Casimir energy/pressure in a dissipative thin film can also be determined by computing the works of cohesion or adhesion in creating the multilayer system from component thin films. Unlike DLP’s method, in which the total Casimir energy/pressure is calculated by summation over discrete Matsubara frequencies along the imaginary frequency axis (ξ-axis), the method described in this paper can proceed as integration over the real frequency axis (ω-axis) or summation over the imaginary frequency axis. The advantage of working with real frequencies is that the division into propagating and evanescent waves based on the magnitude of the in-plane wavevector is unambiguous. By dividing the frequency spectrum of the average stress tensor in the vacuum cavity into contributions from propagating and evanescent waves, it can be shown that the Casimir energy/pressure in a multilayer system of planar films with only dissipative materials is only due to evanescent waves. Propagating waves contribute to the Casimir energy/pressure only if at least one of the constituent films is vacuum. Extension of the technique described in this paper to objects of general shapes and to thermal nonequilibrium conditions will help us better understand Casimir energy/pressure in dissipative materials [34, 35] and the Casimir-Lifshitz forces between coated nanosurfaces [36].



---

*Source: 198657-2015-03-25.xml*