Canonical quantization of macroscopic electromagnetism

Application of the standard canonical quantization rules of quantum field theory to macroscopic electromagnetism has encountered obstacles due to material dispersion and absorption. This has led to a phenomenological approach to macroscopic quantum electrodynamics where no canonical formulation is attempted. In this paper macroscopic electromagnetism is canonically quantized. The results apply to any linear, inhomogeneous, magnetodielectric medium with dielectric functions that obey the Kramers-Kronig relations. The prescriptions of the phenomenological approach are derived from the canonical theory.


Introduction
Macroscopic electromagnetism is an indispensable part of theoretical physics, providing an accurate description of a large variety of light-matter interactions [1,2]. The ability to encompass the electric and magnetic responses of media within macroscopic dielectric functions, without reference to the microscopic material structure, is of great practical importance. A microscopically complete treatment of the interaction of light with most real materials is currently impossible, and even if the computer resources necessary for such a task were available, simple analytical calculations using the macroscopic theory would still be adequate in many circumstances of interest.
But the source of the strength and versatility of macroscopic electromagnetism, namely the removal of microscopic material degrees of freedom, has been seen as a disadvantage when it comes to quantizing the theory [3,4]. A rigorous quantization procedure would require a Lagrangian and Hamiltonian formulation of the theory, followed by the standard canonical quantization rules. Such a canonical formulation of macroscopic electromagnetism has so far been thwarted by the presence in all media of dispersion and absorption, an obstacle that cannot be avoided if one considers all frequencies of electromagnetic waves for which the macroscopic theory is valid. This difficulty has led to a phenomenological approach to macroscopic quantum electrodynamics (QED), in which no canonical formulation is attempted (some details of this scheme will be described below; for a comprehensive description see [5,6] and references therein). The lack of a canonical basis for macroscopic QED is ameliorated by the results of calculations that introduce a simple microscopic model for a dielectric [4,7,8,9]. Huttner and Barnett [4] showed that the canonical quantization of the resulting coupled light-matter system leads to results that can be expressed in terms of a homogeneous electric permittivity, derived from the matter degrees of freedom, at which point the model reproduces relations introduced in macroscopic QED on a phenomenological basis [5,6]. The simple microscopic model has been extended to the case of an inhomogeneous dielectric [7,8] and to a magnetodielectric medium [9], again producing results at the level of the dielectric functions that agree with the prescriptions of the phenomenological theory [5,6]. An essential ingredient of the microscopic models is the introduction of a reservoir in addition to the electromagnetic and matter degrees of freedom, in order to incorporate the dissipation that is necessarily present due to the Kramers-Kronig relations.
The drawback of the microscopic models as a justification of the phenomenological approach is, of course, their simplicity.
As a classical theory, macroscopic electromagnetism is experimentally proven to be valid in the case of highly complicated materials, and a rigorous quantization would ideally also proceed on a macroscopic level.
This requires a canonical formulation and quantization of macroscopic electromagnetism, a task that is performed here. The resulting theory is valid for all linear, inhomogeneous, magnetodielectric media with dielectric functions that obey the Kramers-Kronig relations. The treatment is macroscopic, in that the electromagnetic fields are shown to obey the macroscopic Maxwell equations for an arbitrary linear medium, but a reservoir is necessary to account for dissipation, as in the microscopic models. Indeed, it is clear from the outset that a reservoir must be included in addition to the electromagnetic fields, since otherwise one would obtain a Hamiltonian involving only the latter, and thereby a conserved electromagnetic energy in dissipative media. But the necessity of including a reservoir does not imply that additional dynamical degrees of freedom for the medium must also be introduced; only the electromagnetic fields and the reservoir are required. This last point has been recognized in some previous treatments of QED in media [10,11,12,13,14,15]; we here give a comprehensive derivation of general, linear macroscopic QED from electromagnetic fields coupled to a reservoir. Due to the common feature of a reservoir, the details of the canonical macroscopic theory have much in common with the microscopic models [4,7,8,9]. In fact, it is no more difficult to canonically quantize macroscopic electromagnetism than to quantize the simple microscopic models. But the canonical theory is a greater reward for a similar amount of labour, as it is applicable to any medium where the classical macroscopic theory is valid.
The postulates of the phenomenological theory [5,6] that serve as its starting point are here derived from canonical quantization. We thereby remove the need for a phenomenological approach since the canonical theory has the same range of validity (in contrast to the microscopic models). As well as providing a foundation for the wide range of important results that have been obtained from the phenomenological prescriptions [5,6], the canonical theory will allow investigation of issues that require a Hamiltonian and the related canonical machinery to be properly addressed. An important example where the benefits of a canonical theory are desirable is the Casimir-Lifshitz effect [16,17,18,19], the phenomenon of forces on bodies due to electromagnetic zero-point and thermal fields. The standard Lifshitz theory [17,18,19] of this effect, being macroscopic and lacking a canonical foundation, suffers from the drawbacks of the phenomenological approach to macroscopic QED in general. This has led to doubts about the status of Lifshitz theory as a proper quantum treatment; the lack of an intelligible Hamiltonian has been criticized ‡ and it has recently been argued that the theory is essentially classical [20]. The canonical macroscopic QED developed here can provide a rigorous foundation for Lifshitz theory.
The action principle from which all the results follow is written in Section 2 and shown to lead to the macroscopic Maxwell equations. In Section 3 the standard rules of canonical quantization are applied and the Hamiltonian is derived. The Hamiltonian is diagonalized in Section 4 and the electromagnetic field operators are written in terms of the diagonalizing operators; at this point the postulates of the phenomenological approach emerge.

The action and the field equations
In the canonical formulation of electromagnetism [21] the dynamical fields must be taken to be the scalar potential φ and vector potential A, related to the electric and magnetic fields by The relationship between fields in the time and frequency domains is, for the example of the electric field, In limited frequency ranges where losses are negligible the action for dispersive macroscopic electromagnetism can be written as a functional of φ and A only [22]. Our concern here, however, is with the theory in its full generality, where all frequencies are included and the dielectric functions obey the Kramers-Kronig relations [2,1]. This means that in addition to φ and A the Lagrangian must include dynamical variables that account for absorption and dissipation of the electromagnetic energy. We incorporate dissipation in exactly the manner used in the microscopic models [4,7,8], by introducing a continuum of reservoir oscillator fields. As we allow for a magnetic as well as an electric response in the macroscopic medium, the reservoir will be a set of two such continua of oscillator fields, which we denote by X ω (r, t) and Y ω (r, t), where ω > 0 is a continuous label giving the frequency of oscillation. The dynamical variables are ‡ I am indebted to Gabriel Barton for properly stressing this deficiency of Lifshitz theory.
thus φ, A, X ω and Y ω , and we must write an action that is a functional only of these fields and whose dynamical equations are the macroscopic Maxwell equations in a general inhomogeneous magnetodielectric medium. As there are no dynamical variables representing the medium, it must be included only through the coupling of the electromagnetic fields with the reservoir. The required action is where S em is the free electromagnetic action: S X and S Y are the actions for the free reservoir oscillators: and S int is the interaction part of the action, coupling the electromagnetic fields to the reservoir: α(r, ω) = 2ε 0 π ωε I (r, ω) In the coupling functions (8), ε I (r, ω) is the imaginary part of ε(r, ω), the relative permittivity of the desired macroscopic medium, and κ I (r, ω) is the imaginary part of κ(r, ω) = µ(r, ω) −1 , where µ(r, ω) is the relative permeability. For dissipative media, ε I (r, ω) > 0, whereas κ I (r, ω) < 0 because µ I (r, ω) = −κ I (r, ω)/|κ(r, ω)| 2 > 0. The coupling functions (8) are thus real and so is the action (3); this will ensure that the Hamiltonian of the quantized theory is hermitian. The entire dielectric functions ε(r, ω) and κ(r, ω) of the medium, not just the imaginary parts, are specified by (8) because the real parts are given by the following Kramers-Kronig relation [1,2]: and similarly for κ(r, ω).
It should also be remembered in what follows that ε I (r, ω) and κ I (r, ω) are odd functions of ω, whereas ε R (r, ω) and κ R (r, ω) are even, so that ε * (r, ω) = ε(r, −ω) and κ * (r, ω) = κ(r, −ω), which is a consequence of the electric and magnetic susceptibilities being real [1]. We have assumed the medium is isotropic, with scalar dielectric functions ε(r, ω) and κ(r, ω); anisotropy can be included by obvious modifications, but this is not done here as it would further tax an already burdened notation in what follows. The coupling of electromagnetic fields to a reservoir has been considered previously in [10,11,12,13,14,15]. Here we prove that the specific action (3)-(8) gives the macroscopic Maxwell equations for the linear, isotropic, inhomogeneous medium with dielectric functions ε(r, ω) and κ(r, ω). Variation of φ, A, X ω and Y ω in the action (3) gives, respectively, In order to find the independent equations satisfied by the electromagnetic fields we must solve (12) and (13) for the reservoir fields X ω and Y ω in terms of E and B, and substitute the results into (10) and (11). To solve (12) we write it in the frequency domain using (2) and obtain The general solution of (14) is where Z ω (r) is an arbitrary function independent of ω ′ that gives the solution of the homogeneous equation ( (12) with E = 0), and 0 + is an infinitesimal positive number introduced to deal with the pole at ω ′ = ω that arises from dividing (14) by ω ′2 −ω 2 . The prescription chosen in (15) for dealing with the pole amounts to a choice of boundary condition for the reservoir field X ω ; in the microscopic models [4,7,8], which also include a reservoir, similar poles are encountered in solving the coupled system of field equations and a choice is required that leads to the desired class of solutions. Note that the prescription for the pole, and the form of the solution of the homogeneous equation, must respect the frequency-domain property X * ω (r, ω ′ ) = X ω (r, −ω ′ ) that follows from the reality of the field X ω (r, t) in the time domain; this property holds in (15). In transforming (15) to the time domain the integration over frequency is determined in terms of principle values [23] by The restriction to ω > 0 in the last line is appropriate because ω is positive in (15) (see the action (3)- (7)). Note again that after insertion of (16) in (15) the property X * ω (r, ω ′ ) = X ω (r, −ω ′ ) is seen to hold. Using (2) in the case of X ω , we obtain from (15) and (16) To evaluate the X ω -dependent terms in the field equations (10) and (11) we require the following integral, which we evaluate using (8) and the Kramers-Kronig relation (9): The second term on the right-hand side of the final equality in (18) is the electric displacement D(r, t): Equation (13) for Y ω can be solved in an identical manner to that used to obtain (17); denoting the arbitrary function, corresponding to Z ω (r) in (17), by W ω (r), the solution To evaluate the Y ω -dependent terms in the field equations (10) and (11) we require the following integral, which is analogous to (18) and is evaluated in an identical manner: The second term on the right-hand side of (21) is the negative of the magnetic H(r, t) field: Substitution of (18), (19), (21) and (22) into the electromagnetic field equations (10) and (11) yields the macroscopic Maxwell equations where the free charge density σ and free current density j are given by The charge and current densities automatically obey the conservation law and the remaining two Maxwell equations are also identities because of (1). We have thus proven that that (3)-(7) is the action for macroscopic electromagnetism. The electric and magnetic fields are found in the usual manner. Because of (1), a time derivative of (24) yields the following frequency-domain equation for E(r, ω): The electric field can then be written as where the Green bi-tensor G(r, r ′ , ω) is the solution of with retarded boundary conditions. Due to (1) the magnetic field is given in terms of E(r, ω) by B(r, ω) = −i∇ × E(r, ω)/ω.

Quantization
The canonical formalism [21] is readily applied to the action (3)- (7). A degree of choice is present in writing the Lagrangian density L, associated with the action (3) by since integrations by parts in (32) change the form of L without affecting the dynamics. We take L to be given by the integrands (3)-(7), without any integrations by parts. The canonical momenta are then The familiar complications of canonical QED [21] are present here: the canonical momentum associated with the scalar potential φ vanishes and the gauge freedom also reduces the number of independent electromagnetic degrees of freedom. These difficulties are dealt with [21] by eliminating φ as a dynamical variable, replacing it by its solution obtained from the field equations, and by making a choice of gauge. Even after φ is removed from the list of canonical variables, the vanishing of Π φ is responsible [21] for a constraint on Π A : from the field equation (10) and the expression for Π A in (33) we obtain The general method for quantizing field theories with constraints, such as (35) and a choice of gauge, is that of Dirac brackets [21]. For unconstrained theories the classical Poisson brackets of the canonical variables are to be replaced by −i/ times quantum commutators, but when there are constraints on the canonical variables it is a generalization of the canonical Poisson brackets, known as Dirac brackets, that become commutators in the quantized theory. The Dirac brackets are computed from the Poisson brackets of the quantities that vanish due to the constraints [21]; this requires that the constraints be written in terms of the canonical variables. Here (35) is one constraint, and a choice of gauge will give another. For free-space QED the Coulomb gauge is a convenient choice, but the macroscopic Maxwell equation (23) suggests a generalization that has the frequency-domain form In homogeneous media (37) reduces to the Coulomb gauge (36), but the gauge choice (37) is not appropriate for inhomogeneous media because in the time domain it cannot be written in terms of the dynamical fields and their canonical momenta. As a result, Dirac brackets involving the constraint (37) cannot be calculated and there is no obvious way to proceed with the quantization. We therefore choose the Coulomb gauge (36), even in the general case of inhomogeneous media.
After the elimination of φ as a dynamical field and imposition of the Coulomb gauge, the electromagnetic canonical variables are A and Π A , subject to the constraints (35) and (36). These are the same constraints as in free-space QED in the Coulomb gauge [21]; an evaluation of the Dirac brackets arising from these constraints shows [21] that the quantum canonical operatorsÂ andΠ A should be defined by the equal-time commutation relations where δ Tij (r − r ′ ) is the transverse delta function [21,24] δ Tij (r − r ′ ) = δ ij δ(r − r ′ ) The canonical operators of the reservoir are subject to no constraints and therefore satisfy the commutation relations BothÂ andΠ A are transverse due to the constraints (36) and (35) and this is consistent with the commutation relations (38). From (1), (33) and (35) we see that the transverse part of the electric fieldÊ is equal to −∂ tÂ , while the longitudinal part ofÊ determines the scalar potentialφ and is given by the longitudinal part of ∞ 0 dω α(r, ω)X ω : To write the Hamiltonian we must express the time derivatives of the dynamical fieldsÂ,X ω andŶ ω in terms of their canonical momenta; this is done for the reservoir fields by (34), while for the vector potential we obtain from (33) and (1) The Hamiltonian densityĤ is given bŷ To express the HamiltonianĤ = d 3 rĤ in terms of the canonical variables, all time derivatives in (47) must be eliminated using (34) and (46). After substitution of (46) into the first term of the Hamiltonian density (47) we encounter a term −Π A · ∇φ; upon taking the volume integral d 3 rĤ this term becomes, after an integration by parts, d 3 rφ∇ ·Π A , which vanishes by (35). Alternatively, the term − d 3Π A · ∇φ can be seen to vanish by noting that ∇φ is a longitudinal vector, whereasΠ A is transverse by (35); the volume integration of the scalar product of a transverse with a longitudinal field necessarily vanishes [24]. In the Lagrangian-density term in (47) the electric field can be expressed in term of the canonical variables by (33) and the magnetic field by (1). The resulting Hamiltonian iŝ The Hamiltonian (48) gives dynamical equations for the field operators that are the quantum versions of the classical equations (10)- (13). First note that (10) is now a statement of the constraint (35) (see (33)), rather than a dynamical equation. The Hamiltonian (48) and the commutation relations (38) determine the time derivative of A as In this equation the transverse part of ∞ 0 dω α(r, ω)X ω appears because its longitudinal part does not contribute the volume integral ofΠ A · ∞ 0 dω α(r, ω)X ω in (48) (recall that Π A is transverse). The time derivative ofΠ A is (49) and (50) gives the dynamical equation forÂ, which can be written in terms of the transverse part ofÊ by (45): This is simply the quantum version of the transverse part of (11); the longitudinal part of (11) appears here as the time derivative of the second equation in (45), which arises from the constraints. The full equation (11) thus holds also in the quantum theory. It is straightforward to verify, in a similar manner, that the Hamiltonian (48) gives dynamical equations forX ω andŶ ω that are simply the quantum versions of (12) and (13), respectively. It therefore follows from the results of the last Section that the quantum theory specified by the commutation relations (38)-(43) and the Hamiltonian (48) gives the quantum macroscopic Maxwell equations (23)-(24) forD andĤ. These quantum Maxwell equations will also be explicitly derived in the following Section.

Diagonalization of the Hamiltonian
We now proceed to diagonalize the Hamiltonian (48), which essentially amounts to solving the Hamilton equations for the quantum canonical operators. This diagonalization is a familiar procedure from quantization of the microscopic models [4,7,8,9], and we follow the method used in [8], where the model was of an inhomogeneous dielectric. The field operators in the diagonalized Hamiltonian will be bosonic creation and annihilation operators for the eigenmodes of the system. We will show that the Hamiltonian can be written whereĈ e (r, ω) andĈ m (r, ω) are annihilation operators associated with the electric and magnetic responses of the medium, respectively. The diagonalizing operators obey the usual commutation relations of bosonic annihilation and creation operators: From (52) and (53) we obtain the commutators Ĉ λi (r, ω),Ĥ = ωĈ λ (r, ω), which imply that the time-dependent operators are given bŷ The canonical operators in the Hamiltonian (48) must be expressible as a linear combination of the diagonalizing operators; for example, it must be possible to write the vector potentialÂ and it canonical momentumΠ A aŝ where f λ A (r, r ′ , ω) and f λ Π A (r, r ′ , ω) are c-number bi-tensors. Relations similar to (56) and (57) must also hold for the remaining canonical operators in (48); thus forX ω andΠ Xω we requireX with analogous formulae forŶ ω andΠ Yω in terms of bi-tensorial coefficients f λ Y and f λ Π Y . The bi-tensors f λ A and f λ Π A are transverse because of (36) and (35). From (33) we see that (56) and (58) imply an expansion of the electric fieldÊ similar to (56), in terms of a bi-tensorial coefficient f λ E related to f λ A and f λ X by To establish the diagonalization we must show that there exist bi-tensorial coefficients in (56)-(59) (and in the corresponding equations forŶ ω andΠ Yω ) such that the Hamiltonian (48) takes the form (52). From (53) and (56)-(59) it follows that the desired bi-tensorial coefficients must be given by commutators of the canonical operators with the diagonalizing operators; we give two examples: The diagonalizing operators must in turn be expressible as linear combinations of the canonical operators (the inverse of the relations (56)-(59)) and the collection of commutators of type (61), together with (38)-(43), imply thatĈ λ (r, t, ω) has the expansion where the asterisks on the bi-tensorial coefficients denote complex conjugation. We can now substitute (62) and (48) into (54) and use (38)-(43) to obtain an expansion of ωĈ λ (r, t, ω) in terms of the canonical operators; when the coefficients in this expansion are compared with those in the expansion (62) we find the following conditions on the bi-tensorial coefficients: The last term in (64) is the transverse part in respect of the first index on this bitensor; the transverse apart appears because this quantity was the coefficient of a transverse operator in a volume integral. Note that (63) and (64) are consistent with the transverseness of f λ A and f λ Π A , mentioned after (59). We must now solve (63)-(68) to find the required bi-tensorial coefficients.
Equations (65) and (66) imply where (60) has been used in the second equality. Equation (69) is essentially the same as (14), and we write its general solution with the same pole prescription as in (15): The last term in (70) is the solution of the homogeneous equation ((69) with f λ E = 0) featuring an arbitrary bi-tensor h λ X . We do not need to include a term proportional to δ(ω + ω ′ ) in (70), as is done in (15), because we need only consider positive frequency arguments in the bi-tensorial coefficients. From (67) and (68) we find a solution for f λ Y similar to (70): Elimination of f λ Π A from (63) and (64) gives The integrations in (72) can be performed using (70) and (71), in the same way as the integrals (18) and (21) were evaluated in Section 2. We use the principle-value relation (16) in (70) and (71), where here the term proportional to δ(ω + ω ′ ) does not contribute as we take only positive frequency arguments in the bi-tensorial coefficients. The definitions (8) and the Kramers-Kronig relation (9) are also required, and we obtain Equations (60) and (64) show that which implies Insertion of (73), (74) and (75) in (72) and application of (76) produces This equation is essentially the same as the wave equation (28)-(29) for the electric field in the frequency domain, and its solution can be written as in (30) using the Green bi-tensor (31): All of the bi-tensorial coefficients can now be written in terms of h λ X and h λ Y ; to avoid overly lengthy expressions we will write all the coefficients in terms of h λ X , h λ Y and f λ E , bearing in mind that f λ E is itself determined by h λ X and h λ Y through (78)-(79). The right-hand sides of (63)-(68) are written in terms of h λ X , h λ Y and f λ E using (70), (71), Equation (77) shows that the right-hand side of (80) is transverse, a fact that is used in deriving (80) and (81). The pole prescription in (70) and (71) has been rewritten in (82) and (84) as shown in the first equality in (16). All that remains is to deduce the bi-tensors h λ X and h λ Y ; once these have been specified the bi-tensorial coefficients are found from (78)-(85). The conditions that determine the form of h λ X and h λ Y follow from substituting the expansion (62) and its hermitian conjugate into the commutation relations (53). When this substitution is made in the first of (53), and the relations (80), (81), (83) and (85) are used, we obtain the requirement The quantities in square brackets in the first two lines of (86) are transverse because of (77); they therefore project out the transverse parts of the bi-tensors f λ ′ E and f λ * E with which they are combined in a scalar product. For this reason the transverse restriction on f λ ′ E and f λ * E , which was present because of (81), has been removed in (86). To find the condition on h λ X and h λ Y that results from (86), we must insert (82) and (84) and perform a simplification that makes use of (77); the details of this tedious calculation are described in Appendix A and the result is We must also impose the condition arising from the substitution of (62) into the second commutation relation in (53). But this condition is found to be automatically satisfied when (80)-(85) and (77) are utilized, and we relegate the details of this calculation to Appendix B. The diagonalization (52) is achieved once we find bi-tensors h λ X and h λ Y that satisfy (87). The diagonalization operators are of course not unique since, for example, a unitary transformationĈ λj (r, ω) → U ijĈλj (r, ω), U † ik U kj = δ ij does not change the form (52), and similarly a unitary transformation can be applied to h λ X and h λ Y without affecting the condition (87). A suitable solution of (87) is which completes the diagonalization procedure.
We can now write the electromagnetic field operators in terms of the diagonalizing operatorsĈ λ (r, t, ω). The bi-tensor s λ in (78)-(79) is seen from (88) to be s e (r ′ , r, ω) = ω 2 ε 0 π ε I (r ′ , ω) which determines through (78) the bi-tensor f λ E . Having obtained f λ E , the electric-field version of (56) together with (55) show that The current operator defined in (92) is easily seen from (53) to obey the commutation relations ĵ (r, ω),ĵ(r ′ , ω ′ ) = 0, where the notation × ← ∇ ′ denotes a curl with respect to the right-hand index, so that V(r)× ← ∇ = ∇ × V(r) for a vector V(r). The transverse and longitudinal parts ofÊ(r, t) are related to the electromagnetic potentials by (45) and the magnetic field operator can be found from which is an identity in terms of the electromagnetic potentials.
Equations (91) and (92) are the quantum versions of the classical results (30) and (29). The diagonalizing operatorsĈ e (r, ω) andĈ m (r, ω) are thus seen to be closely related to the classical fields Z ω (r) and W ω (r) in Section 2. It is easy to show that (91) is the solution of the quantum macroscopic Maxwell equations. From (91) and (95), and the defining equation (31) of the Green bi-tensor, it follows that with use of the definitions (19) and (22). Similarly, (91), (19) and (31) imply where in the frequency domain the charge density iŝ so that the quantum version of the conservation law (27) holds.
An interesting issue raised by Huttner and Barnett [4] is the possibility of negativeenergy eigenstates in their microscopic model if the coupling constants are too large, since this can happen in the coupling of harmonic oscillators. Huttner and Barnett [4] found only positive frequencies in the diagonalized Hamiltonian, without restrictions on the size of the coupling constants in their theory, in which only the transverse part of the vector potential is quantized. This conclusion was supported by verifying that operator expansions analogous to (56)-(59) can be inverted to give expansions analogous to (62), with only positive frequencies present in the expansions. Here we have similarly found that only positive frequencies are required in the diagonalized Hamiltonian (52). The relation of these results to the case of general coupled harmonic oscillators and the occurrence of negative energies may be interesting to investigate.

Conclusions
The absorption of electromagnetic energy in media is a necessary consequence of the restriction to retarded solutions of Maxwell's equations, which leads to the Kramers-Kronig relations. This means that the canonical quantization of electromagnetism in media cannot proceed by treating only the electromagnetic fields as dynamical variables, unless restrictions are made in the frequency range of applicability. By including a reservoir in addition to the electromagnetic fields, we developed macroscopic QED using the standard canonical formalism of quantum field theory. The theory applies to any linear, inhomogeneous, magnetodielectric medium with dielectric functions satisfying the Kramers-Kronig relations. The starting point of the phenomenological approach to macroscopic QED was derived from the underlying canonical theory. Important results such as the Lifshitz theory [17,18,19] of vacuum and thermal electromagnetic forces can therefore be provided with a canonical foundation.