Quantum Electromagnetic Stress Tensor in an Inhomogeneous Medium

Continuing a program of examining the behavior of the vacuum expectation value of the stress tensor in a background which varies only in a single direction, we here study the electromagnetic stress tensor in a medium with permittivity depending on a single spatial coordinate, specifically, a planar dielectric half-space facing a vacuum region. There are divergences occurring that are regulated by temporal and spatial point-splitting, which have a universal character for both transverse electric and transverse magnetic modes. The nature of the divergences depends on the model of dispersion adopted. And there are singularities occurring at the edge between the dielectric and vacuum regions, which also have a universal character, depending on the structure of the discontinuities in the material properties there. Remarks are offered concerning renormalization of such models, and the significance of the stress tensor. The ambiguity in separating"bulk"and"scattering"parts of the stress tensor is discussed.


I. INTRODUCTION
Most studies of the Casimir effect deal with quantum fluctuation forces between rigid bodies separated by vacuum. Such forces are finite and can be calculated exactly, in principle. (For reviews, see, for example, [1][2][3].) Casimir's original configuration was that of perfectly conducting plates in otherwise empty space [4]. This was generalized by Lifshitz to dielectric slabs, but again they were separated by vacuum [5]. The addition of Dzyaloshinskii and Pitaevskii was essential to the replacement of the intervening vacuum by a homogeneous medium [6]. The resulting theory has been remarkably successful, and was confirmed by the verification of the attractive force of a helium film by a substrate [7,8], well before the modern demonstration of the vacuum Casimir force [9]. The theory has been applied to a wide variety of fields [10][11][12][13][14][15].
The local Casimir energy density and other components of the stress tensor have also been intensively investigated. These exhibit well-known behaviors near the surfaces of the bodies. (For a review of some of the literature on this, see Ref. [16].) This is relevant, not only for a deeper understanding of the Casimir force, but fundamentally for the coupling to gravity; in simple contexts, the local Casimir stress tensor has been shown to be consistent with the equivalence principle, including the divergent contributions [17]. Consistent results for finite Casimir stress tensor components were earlier obtained in Refs. [18,19].
At least formally, separating rigid bodies by a uniform dielectric leads to no difficulties in computing vacuum forces, and even dispersion can be incorporated, although including dissipation may present challenges. However, the situation is much less clear when the bodies are immersed in an inhomogeneous medium. There have been various attempts to describe Casimir forces with nonuniform dielectrics [20][21][22]. The most ambitious treatment of the inhomogeneous electromagnetic Casimir problem seems to be that of Griniasty and Leonhardt [23,24], who examine the local stress tensor and propose a specific renormalization scheme to remove the divergences that occur in such circumstances. For the case of a one-dimensional slab with a dielectric response that varies smoothly except for a discontinuity in the slope as one enters the material, they find a universal singularity behavior in the normal-normal component of the vacuum expectation value of the stress tensor at the boundary between vacuum and the dielectric.
For some years we have been investigating similar issues, but in the scalar field context [25][26][27][28]. In particular, using a WKB analysis, we identified the universal Weyl divergences in the stress tensor components for an arbitrary semi-infinite slab described by a potential v(z), where z is the distance into the slab. For particular cases (a linear or a quadratic wall) we also examined how the remainder of the stress tensor, after the divergent and growing terms are removed, behaves near the edge. In this connection the work of Mazzitelli et al. should be mentioned [29,30]. (For more references, see the appendix of Ref. [27], and also Ref. [31], which should have been included there.) Very recently, we have made further progress in understanding how the divergences are to be renormalized [32].
In the present paper, inspired by the remarkable results of Ref. [24], we generalize our considerations [25][26][27][28] of the local stress tensor in one-dimensional geometries to the electromagnetic case, in which the role of the potential is played by the permittivity. More precisely, the deviation of the permittivity from its vacuum value will be referred to as the potential in this paper. In the next section, we review the difficulty of formulating the stress tensor in inhomogeneous media, and derive the non-conservation law satisfied classically by the spatial stress tensor. In Sec. III we show how the Green's dyadic for this problem breaks up into transverse electric (TE) and transverse magnetic (TM) parts. We also write down the construction of the various components of the stress tensor in terms of the TE and TM Green's functions. This also includes the correct dispersive factor for the energy density [33].
The generic set-up of the problem is given in Sec. IV, including the break-up of the Green's functions into "scattering" and "bulk" parts, referring to the contributions from the outgoing wave and incoming wave contributions. This break-up, of course, is not unique. An example, the reflectionless potential considered in Ref. [24], is treated somewhat more generally in Sec. V A. There we show, using the uniform (Debye) asymptotic expansions for the modified Bessel functions, that there are two types of singularities in the normal-normal component of the stress tensor occurring at the edge between the vacuum and dielectric region: a cubic singularity if there is a discontinuity in the permittivity, and a quadratic one (coinciding with that found in Ref. [24]) if only the derivative of the permittivity is discontinuous. We also show that the bulk term (the term independent of the reflection coefficient) contains the expected leading Weyl divergence, as well as further divergences involving the potential, which are regulated by point-splitting.
A second example for which the TE and TM Green's functions may be exactly found is given in Sec. V B. The same edge behavior is found as in Sec. V A for the continuous case. This behavior is evidently universal, as claimed by Ref. [24], and we demonstrate that explicitly in Sec. VI A, using a general perturbative expansion of the Green's functions. All of the above neglects dispersion. In Sec. VI B we discuss the more realistic plasma model, which results in the elimination of the edge singularity in the normal-normal stress, but yields the divergence structure for the bulk contribution coinciding with that for the scalar case considered in Ref. [28]. For the plasma model of dispersion, the TE Green's function is identical with the scalar one.
Other components of the stress tensor are considered in Sec. VII. Again, for the plasma model, the divergences arising from the bulk term in the Green's function coincide with those found for the scalar situation for both TE and TM modes, and the edge singularity for the TE mode for the energy density coincides with that found for the canonical scalar energy density in Ref. [28], while the TM mode has a different numerical coefficient.
The break-up into bulk and scattering parts is not unique, because we can always add an arbitrary admixture of the exponentially suppressed fundamental solution to the exponentially growing one. We attempt to explore this further in Sec. VIII, for the TE mode, which can be exactly solved for a potential that depends on the z coordinate linearly. Numerically, we show that the scattering part of the energy density and the normal-normal component of the stress tensor rapidly go to zero as the dielectric is penetrated, the former exhibiting the expected edge singularity. If an admixture of the first solution is added to the second, the edge singularities do not change, but the behavior inside the dielectric is altered, but still tending to zero as one goes deeply within the material. Only if the scattering part of the Green's function is completely suppressed (a set of measure zero in parameter space) does the qualitative (and quantitative, for the divergences and edge singularities) behavior change.
We finally consider a situation with mirror symmetry in Sec. IX. Here we consider two reflected potentials meeting at z = 0 so there is no vacuum region. In this case, not surprisingly, the edge singularity is doubled. Concluding remarks are offered in Sec. X. In Appendix A we explain the point-split regulation we use in this paper, while in Appendix B we develop the perturbation theory for a potential which is both continuous and has a continuous first derivative, but where the second derivative is discontinuous.
In this paper we use Heaviside-Lorentz electromagnetic units, and = c = 1.

II. FORCE ON DIELECTRIC
From the Maxwell-Heaviside equations we can derive the statement of electromagnetic momentum conservation. We follow Sec. 7.1 of Ref. [34]. Equation (7.10) there says that is the force density on the charged particles, and the field momentum is Here, a summation convention is used for repeated indices, and ρ and j are the free charge and current densities.
To what extent is the right side of Eq. (2.1) the negative of a total divergence, −∇ · T, which would imply a local conservation law of momentum? As usual it is convenient to do a Fourier (frequency) transform of the fields (we will here suppress the spatial coordinates), assuming a linear medium. For the electric fields where we have introduced a frequency-dependent permittivity tensor, ε(ω), which we allow to be spatially varying. Similarly for the magnetic fields, We now take the average over a time T large compared to atomic time scales but short compared to macroscopic times, so the dyadic product can be written, for example, as Then, in the absence of dissipation, we use the Hermiticity property arising from the reality of the constitutive relations in spacetime, ε ij (ω) = ε ji (−ω) = ε ji (ω) * . 1 If the permittivity and permeability were independent of position, there would be an averaged macroscopic stress tensor, However, if the electrical properties depend on position, this is not the case, but, rather, the right side of Eq. (2.1) would be For a recent review concerning electromagnetic stress tensors see Ref. [35]. For example, consider a dielectric body (µ = 1) immersed in a static classically imposed electric field. Because there is no time dependence and no free charge, we have where the trace is over the tensor indices, and the notation (∇) is a reminder that the free vector index is on the gradient operator. Suppose the body, which need not be homogeneous, is immersed in a homogeneous medium of permittivity ǫ. The force on the body is the momentum flux into the body, since the local momentum conservation law holds there, where S is a surface that entirely surrounds the body. By the divergence theorem where the spatial integral is over the interior of the body (because the permittivity is constant outside the body). This is a generalization of the familiar formula for the force on a dielectric, Eq. (11.44) of Ref. [34], to which it reduces for the isotropic case. We can immediately generalize this to the Casimir force by replacing in Eq. (2.8) 12) in terms of the Green's dyadic Γ, so that the dispersion force on the dielectric body is tr Γ(r, r; ω)(∇)ε(r, ω). (2.13) Here we have identified 2πδ(0) with the averaging time T . In particular, if the body has a homogeneous dielectric constant ε = ǫ, then ∇ε = −ŝ(ε − ǫ)δ(s − s 0 (r ⊥ )), (2.14) where the surface of the body is given by s = s 0 (r ⊥ ), in terms of a coordinate s (outwardly) normal to the surface. The other coordinates are denoted by r ⊥ . (For the case of a planar body in the x-y plane, s = z.) Thus the Casimir force on the body is given by an integral over the surface of the body, Again, this is an obvious generalization of known formulas. 2 The general form for the nonconservation of the vacuum expectation value of the electromagnetic stress tensor in a medium is This is, of course, quite analogous to the nonconservation equation satisfied by the stress tensor for a scalar field in a background potential [28].

III. GREEN'S FUNCTIONS
In this paper we will consider planar situations in which the permittivity ε(z) and the permeability µ(z) depend only on a single coordinate z. We will also allow ε and µ to depend on frequency. For simplicity, we will henceforth assume that ε and µ are isotropic. It is also convenient to make a Euclidean transformation ω → iζ. The general Green's dyadic obeys an equation which follows from the Maxwell-Heaviside equations, which breaks into two modes, TE and TM modes, denoted by two scalar Green's functions labelled by E and H, respectively. These satisfy the differential equations The spatial Fourier components of Γ, defined by are given in terms of these two scalar Green's functions, in the coordinate system where k ⊥ has only a component in the x direction (we drop the z, z ′ dependence of g E and g H ): Here ε = ε(z), ε ′ = ε(z ′ ). These are just as given in Refs. [1,36]. The Fourier-transformed electromagnetic stress tensor may also be given in simple form in terms of these two scalar Green's functions. For example, the zz component of the reduced stress tensor is simply where after differentiation, the limit z → z ′ is understood. Let us also record the other diagonal components of the reduced stress tensor. First, the energy density, which must include the dispersive factors: To preserve the symmetry between the transverse components of the reduced stress tensor, we rotate γ to a general coordinate system. Doing so does not affect t 00 and t zz , but yields after using the equations of motion (3.2) There are also off-diagonal terms, linear in k x or k y , which would vanish upon regulated integration, if that regulation respects the two-dimensional rotational symmetry of the problem. Such a regulator reduces t xx and t yy to The four-dimensional trace is zero if there is no dispersion.

IV. GENERIC PLANAR PROBLEM
To save typographical space, we use comma-separated notation, (µ, ε) and (E, H) to write the TE and TM mode expressions in the following. We can construct the Green's functions from the solutions of the homogeneous equations Here we take F to denote a solution that does not diverge for z → ∞ (typically goes to zero), while G is an arbitrary independent solution. The Wronskian of these two solutions is We want to solve the Green's function equations (3.2) in terms of these solutions, for the situation of a "soft wall", where The solutions are (κ = k 2 + ζ 2 ) Here, the constant α is related to the Wronskian by The reflection coefficients are determined by requiring that g E.H be continuous at z = 0, and that 1 µ,ε ∂ z g E,H also be continuous there. This corresponds to the continuity ofẑ × E andẑ · B, and ofẑ × H andẑ · D. (Imposing these matching conditions requires the form of the Green's function for z > > 0 > z < , not displayed here.) The consequence is Here µ =μ(0), ǫ =ε(0). In the above construction, G is completely arbitrary, save that it be a solution, independent of F , to the differential equation (4.1). Therefore, the reflection coefficientR is not unique, and indeed can be made equal to zero by the replacement G → G −RF . To have a unique reflection coefficient, we need a condition to determine the form of G. Such is supplied by imposing a boundary condition at z → −∞, even though this is outside the region z, z ′ > 0 where the construction (4.4) holds. That is, assuming the continuous functionsε(z),μ(z) hold in all space, so there is no discontinuity, we will henceforth choose G subject to the boundary condition Then the reflection coefficient is uniquely defined. (These boundary conditions as stated here are somewhat schematic; the specific conditions at ±∞ depend on the structure of ε(z).) The stress in the vacuum region, to the left of the wall (z < 0), is immediately calculated from Eq. (3.5): which is independent of z, the term involving the reflection coefficient having cancelled out. This is universally recognized as an irrelevant bulk term, since it has no contribution from the wall, and would be present if vacuum filled all space, so is to be omitted. It is the assertion of Ref. [24] that the same omission is to be done for the contribution to the stress tensor coming from the part of the Green's function in the z > 0 region that is not proportional to the reflection coefficient: In particular, they advocate omitting the stress tensor contribution arising from the term in the Green's function (4.4) 1 α F (z > )G(z < ), even though it is spatially varying, because this term would be there in the absence of the edge at z = 0. This hypothesis may be suspect, but we will follow it for the moment.

V. EXACTLY SOLVABLE EXAMPLES
Now we examine two cases where both the TE and TM modes may be explicitly given. In the first example, the permittivity has a singularity at a finite value of z, which is the natural boundary of the problem, and in the second the permittivity has an exponential behavior.

A. A first example
Let us consider a planar medium described by which has a singularity at z = a. (This is a slightly generalized version of the medium considered in Ref. [24], where the potential was continuous, so ǫ ≡ε(0) = λ/a 2 = 1.) Because of that singularity, the right side of the wall has a finite depth, 0 < z < a; the region z > a is completely disconnected from the region containing the wall. This potential has the virtue of allowing explicit solutions: Here F is chosen to be finite as z → a. Indeed, I ν (0) = 0, K ν (+∞) = 0, consistent with the criteria stated in the previous section. Then the reflection coefficients in the medium arẽ The scattering part of the zz component of the reduced stress tensor (the part proportional to the reflection coefficients) is )) . (5.5) As we wish to examine the stress just inside the wall, we can use the uniform asymptotic expansion (UAE) for the Bessel functions, because it captures the short-distance behavior [37]. That expansion is, as ν → ∞: where u k and v k are polynomials in t = (1 + Z 2 ) −1/2 . The first of these are All we need to know about the functions in the exponents is the derivative: If we retain only the leading factor in the UAE the reflection coefficients are approximatelỹ In the remaining factor of Eq. (5.5) we must keep the O(1/ν) corrections because they are of the same order in κa as the leading term in the stress tensor construction, leaving for the rest of the zz component of the reduced stress tensor . The UAE presumes that the significant values ofκ are large. If ǫ = ε(0) = 1, in the first approximation we may neglect terms of order 1/(κa) and smaller, so the reflection coefficients reduce tõ which has the form familiar from a step discontinuity in the dielectric constant. Further, near the boundary, the exponents combine: which makes use of Eq. (5.8). Finally, to carry out the integrals over frequency and transverse wavevectors we adopt polar coordinates, so that with √ ǫζ =κ cos θ, k =κ sin θ. The angle θ occurs inside the two reflection coefficients, as well as inside the formula for t s zz , Eq. (5.10), since near the wall [u 1 (ka/ν) − v 1 (ka/ν)]/ν = (1 − cos 2 θ)/(2κa), and the integrals of these dependencies for the TE and TM modes give The functions E(ǫ) and H(ǫ) are elementary, given in terms of logarithms, but are not very illuminating to display. Instead we show the plot of them in Fig. 1, and give the limits for small and large values of ǫ − 1: The remaining integral onκ is simple, so after integrating over k and ζ, we are left with the "scattering part" of the zz component of the stress tensor near the wall (z → 0+): The ǫ-dependent factors in the zz components of the stress tensor in Eq. (5.14). The small and large ǫ − 1 limits go out to third order and −7/2 order, respectively. Clearly, the TE contribution is almost insignificant, and the two asymptotic limits accurately cover the full range of ǫ.
And the total zz component of the stress is the sum of these two components, which for the case of a small discontinuity reduces to This cubic singularity disappears if there is no discontinuity, that is, ǫ = 1, whereκ = κ. Then we need to keep the order 1/ν correction in the reflection coefficients as well, so Eq. (5.9) gets modified tõ which is exactly the result found in Ref. [24]. The similarity of the coefficients in Eqs. (5.17) and (5.20), (5.21) is striking. We close this subsection by examining the omitted contribution from the "bulk" term in the interior, It is quite obvious that this does not give singular behavior near the discontinuity in ε(z) at z = 0, but it does yield divergent contributions. The corresponding reduced stress tensor has a form similar to that given in Eq. (5.5): Now, when the UAE is inserted, the cancellation observed in the reflection-dependent part does not occur, so the leading term is For the moment we examine only the leading term in the limit of this expression as z → 0, which is the obvious generalization of Eq. (4.9). When this is integrated over all frequencies and wavenumbers, and regulated by point-splitting as in Ref. [28], we obtain (see Appendix A) exactly the leading bulk divergence seen for each scalar mode in Ref. [28], apart from the expected index of refraction factor. Later we shall encounter the subleading divergences dependent on the potential; beyond them, in the exact t b,E,H there are finite terms that presumably have physical significance.

B. Exponential permittivity
Let us give another exactly solvable model. Consider the permittivity function ε(z) = 1, z < 0, e αz , z > 0. (5.27) For the two modes, the two fundamental solutions to Eq. (4.1) are for z > 0 [38] Again, the second solution is unique, according to the criteria enunciated in Sec. IV, because I ν (0) = 0. In each case, the effective Wronskian (4.5) is the same, Using the UAE, the leading bulk stress tensor component is The scattering part of the reduced stress tensor, near the wall, has the form seen before in Eq. (5.19) if we replace ζ 2 by κ 2 cos 2 θ, and a by 2/α, where α is the slope of the potential at the edge: From this follows the same result for the stress tensor as in Eq. (5.20).
We will see in the following section that this behavior is universal, as long as the potential is continuous and has a linear slope at the edge.

VI. UNIVERSAL EDGE BEHAVIOR
A. First-order perturbation theory Griniasty and Leonhardt [24] asserted that the behavior of the zz component of the subtracted stress tensor seen in Eq. (5.20) is universal. That is, it holds whenever the potential is continuous, but has a discontinuous slope at the origin, the slope being in that case α = 2/a. We will prove that assertion here, which follows from perturbation theory near the edge. We can generalize this slightly, by allowing for a discontinuity ǫ − 1 in the permittivity near the boundary. Sufficiently close to the edge, ε(z) = ǫ(1 + αz), and we will calculate the stress tensor in the approximation that α is very small compared to κ.
We start with the TE mode. The functions F and G satisfy This is easily solved perturbatively for solutions that decay exponentially fast, or that grow exponentially fast, at infinity: keeping terms out through O(α). Since the differential equation contains no first derivatives, the Wronskian remains constant, Using the "bulk" part of the Green's function in the medium, the first term in the second line of Eq. (4.4), we find for the corresponding reduced stress tensor which agrees with Eqs. (5.24) or (5.31a) when they are expanded for small α (fixed z). Integrated over frequency and wavenumbers, we obtain the full bulk stress tensor, when time-splitting, or transverse space-splitting, regulation as in Eq. (5.26) is inserted (see Appendix A), The relative factor of −3 between the linear dependencies of these two forms is the result of the identity given in Ref. [17], reproduced here in Eq. (A6). The reflection coefficient computed from Eq. (4.7) to first order in α is The first term in the parentheses refers to the scattering due to the discontinuity in ε(z) at the edge, while the second term refers to the contribution arising from the slope of the potential. If the latter effect is negligible, this agrees with the form in Eq. (5.11). A bit of algebra shows the "scattering" part of the reduced stress tensor is When Eq. (6.7) is integrated over frequency and transverse wavenumbers according to Eq. (5.13), the result for the stress coincides with that given in Eq. (5.16): which has a first-order perturbative solution Now the Wronskian of the two solutions is not constant, which is exactly what is needed to make α H a constant: The bulk term in the zz-component of the stress tensor turns out to be the same as its TE counterpart (6.5), for example, for time splitting: Now it is straightforward to calculate the scattering part of the stress tensor to O(α 2 ), in terms of the reflection coefficient:R The zz component of the scattering part of the reduced stress tensor is then Again, if for ǫ = 1 we drop the second term in the square brackets, we see the appearance of the TM reflection coefficient for a discontinuity in the permittivity, which leads to the stress tensor as z → 0+:

B. Dispersion
The above assumes that the permittivity does not depend on frequency. This is quite unrealistic. Instead, let's examine what happens if we use a plasma model, where α = α 0 /ζ 2 . This then makes the TE mode coincide with the linear scalar problem considered in Ref. [28]. There the divergent terms were isolated using a WKB approximation.
We can easily reproduce those leading divergences. To compare with the results there, we set the discontinuity ǫ − 1 equal to zero.
With the plasma dispersion relation, the bulk term (6.4) reads, before integration, 19) and then carrying out the frequency and wavenumber integrations using the formulas in Appendix A, we find which are the two leading divergent terms found in Ref. [28] for a linear potential. Perhaps surprisingly, the same holds for T b,H zz . To get the logarithmically divergent term in T b,E zz one might think we would have to work out perturbation theory to second order, which we will do in the next Section. However, the zz component of the reduced bulk stress tensor to second order can be calculated by knowing only the O(α) solutions because we easily see from the definition of the Wronskian that From this follows where we have introduced the abbreviation γ = ζ 2 ǫ/(4κ 2 ). The small δ expansion of which corresponds to the logarithmically divergent term found in Ref. [28]. Again only the first order result is necessary to give the order α 2 contribution to the bulk stress for the TM mode, since the same formula as Eq. (6.21) applies for the TM mode as well. The result is only slightly different from that in Eq. (6.22): This leads to exactly the same logarithmic divergence in the plasma model as in Eq. (6.23). However, to get such terms for the other components of the stress tensor, we need second-order perturbative solutions for F and G, which we will deal with in the following section. As for the scattering contributions, it is evident that due to the softening produced by the plasma dispersion relation, the singular behavior in T s,E zz as the edge is approached from the inside goes away, consistent with the numerical results shown in Fig. 5 of Ref. [28]. (See further discussion in Sec. VIII.) The scattering part of T s,H zz in the plasma model has to be defined with an infrared cutoff, but certainly also does not diverge as the edge is approached.
So we have verified and extended the results of Ref. [24]: For a vacuum interface with a planar dielectric without dispersion, if the permittivity is continuous, but has a linear slope at the edge, the singularities in the normal-normal component of the stress tensor possess a universal 1/z 2 form, where z is the distance from the edge into the medium. If the permittivity is discontinuous, the normal-normal component of the stress tensor has a 1/z 3 singularity, and as shown in Appendix B, the singularity is reduced to logarithmic if the discontinuity is in the second derivative. As we will see in the next Section, the singularities in the energy density are one order higher for a linear discontinuity. Only the behavior of the potential at the edge of the dielectric is necessary to determine the singularities in form and magnitude, but this we have demonstrated through examples and a general perturbative analysis.

VII. OTHER STRESS TENSOR COMPONENTS
Let us now examine other components of the stress tensor, particularly in the continuous permittivity situation. The leading perturbative approximation yields the leading divergent structure, and the leading behavior near the edge. We will consider both the dispersive case with the plasma model, since that agrees, for the TE mode, with the scalar case, and is approximately realistic, as well as the situation when the permittivity is independent of frequency.

A. Leading-order contributions
Including the dispersive factor, the reduced TE energy density for the plasma model, where α = α 0 /ζ 2 , is for small α 0 z (exactly, for a linear potential) which agrees with the scalar energy density for a linear potential provided the conformal parameter ξ = 0 (or in the language of Ref. [28], β = −1/4), surprisingly, not the scalar conformal value of ξ = 1/6. (That is, the canonical stress tensor emerges, not the conformal one.) Thus we see that (setting ǫ = 1) Using the point-splitting methods of the Appendix, we find for the bulk contribution to the energy density, which coincides with the leading divergences found in Ref. [28]. Note that ∂ ∂δ δT ∆ 00 = T τ 00 (7.4) holds for the relation between the energy densities with the spatial and temporal cutoffs, as in Ref. [17]. And in the medium, just to the right of the edge, we find for the scattering contribution which when integrated over frequency and wavenumbers yields exactly the result as for the scalar case with β = −1/4 given by Eq. (6.7) of Ref. [28]. Had we assumed that α was independent of ζ, the sign of the potential term in Eq. (7.1) would have reversed, and we would have obtained instead for the bulk divergence (with temporal splitting) and for the edge singularity in the scattering part more singular than the behavior of T s,E zz in this non-dispersive model seen in Eq. (6.9). For the remaining diagonal components, from Eq. (3.8), for τ splitting, or ∆ splitting, respectively, in the plasma model, which exactly coincides with the leading scalar divergences found in Ref. [28] when we average over ρ x , ρ y there. It is easily checked that the trace identity (3.9) is satisfied: For the scattering part, which is exactly half the energy density found in Eq. (7.6) as required by the trace of the scattering part of the stress tensor being of O(α 2 0 ). For the TM mode in the plasma model the bulk part of the reduced energy density is which, upon integration, leads to the same result as Eq. (7.3). The transverse bulk parts of the reduced stress tensor are leading to the same result as Eq. (7.9), as required by the trace identity. For constant α the energy density divergence is the same as for the TE part, Eq. (7.7). The scattering part of the reduced energy density is which is twice the transverse reduced stress tensor components, as required by the trace identity. This possesses singularities, when ζ 2 = κ 2 cos 2 θ goes to zero, so the meaning of these seems somewhat obscure. However, if we adopt the nondispersive model and assume that α is constant, we can find the energy density singularity near the edge which is −9 times that from the TE mode, Eq. (7.8).

B. Second order perturbation theory
To proceed further, we need to work to the next order in perturbation theory. It is easy to work out the solutions to Eq. (6.1) to second order, assuming the potential is exactly linear. The two solutions are The expansion parameter is αγ. The Wronskian changes, but is still constant, The TM equation (4.1) is, assuming an exactly linear potential, which can also be straightforwardly solved to second order in α: Note that the terms of order αγ and of order α 2 γ 2 coincide with those of the TE solutions in Eq. (7.16). The Wronskian of these two solutions gives Now to get the order-α 2 corrections to the energy density, we have to use the second-order solutions, Eqs. (7.16) and (7.19). A straightforward calculation reveals, for the bulk contributions to the reduced energy density, Note that the O(α 0 ), O(α) and the O(α 2 γ 2 ), O(α 2 γ 3 ) terms are the same for the TE and TM contributions, which means that the divergences in the plasma model are the same, for example, in temporal point splitting for ǫ = 1 as defined in Appendix A, where the remainder is finite as δ → 0. This includes the results already found in Eq. (7.3), and coincides with the scalar divergences found in Ref. [28]. We can also straightforwardly find the next order corrections to the scattering part of the zz component of the reduced TE stress tensor, for example, with ǫ = 1, but the order α 3 correction means that the corresponding term in T s,E zz has one less power of z, so in the constant α situation, through this order, (In the plasma model, recall that there is no singularity in T s,E zz .) Dimensionally, since [α] = 1/L, the higher order corrections to the edge singularity must be subdominant.
Similarly, we can write for the TE part of the reduced energy density through order α 2 , t s,E 00 = −αγ which leads to, for constant α, the energy density through O(α 2 ) Again, the correction is necessarily subdominant.

VIII. EXACT LINEAR TE POTENTIAL
Of course, the linear TE problem is exactly solvable in terms of Airy functions, as seen in Refs. [25][26][27][28]. Independent solutions of Eq. (6.1) are (α = α 0 /ζ 2 )  By using the asymptotic expansion of the Airy functions for large argument, we straightforwardly obtain for the TE reduced scattering Green's function The above is valid if κ 3 /α 0 ≫ 1. If we now regard the potential as weak, we expand in powers of α 0 and obtain through second order This coincides exactly with the Green's function obtained from the perturbative solution (6.2), and leads, for example, to which follows from (6.7) when ǫ = 1 and α = α 0 /ζ 2 . But when one tries to integrate this over wavenumbers and frequency, one encounters an infrared divergence at κ = 0. Of course, such a divergence is not present in the exact solution, since the perturbative expansion is not valid for small κ. In fact, if the exact expression for t s,E zz is integrated the result is finite, but nonzero, at z = 0, as shown in Fig. 2, as earlier stated.
We can do the same type of calculation for the energy density. In this case the energy density does diverge as the edge is approached from within the medium, according to Eq. (7.6). In fact, the numerical integration of the exact formula fits this asymptotic formula quite well for small z, as shown in Fig. 3. T s,E xx has nearly identical behavior, except for the factor of 2 seen in Eq. (7.11).
The above figures were drawn with the assumption that the second solution G was exactly the second Airy function Bi. But, as noted in Sec. IV, the definition of the reflection coefficient is ambiguous, since the second solution may contain an arbitrary admixture of the first. The criteria given in Sec. IV do not apply, because both Ai and Bi behave as damped oscillatory functions for large negative z. The addition of the second solution is typically asymptotically exponentially subdominant, so this ambiguity does not appear in the asymptotic estimates. However, the ambiguity will affect the behavior away from the edge. We investigated this by substituting in the reflection coefficient Bi → Bi +λ Ai, where λ is a constant. (In fact, λ could be a function of κ.) In Fig. 4 we show how agreement with the estimate (7.6) is greatly improved by the choice of λ = 2/π. The reason for this particular value agreeing with the asymptotic estimate is, at present, mysterious.
The edge singularity is not altered when different constant values of λ are used as compared to the perturbative result because the leading asymptotic behavior of the Airy functions is so that when these are used for large κ and fixed z we see that the admixture parameter is related to the perturbation theory one by Here, the latter parameter is defined in the language of Sec. VI A by taking the second solution to be Thus, it is evident that the admixture of the first solution will be exponentially suppressed within the wavenumber integral.
The comparison between the perturbative value of the reflection coefficient and the exact one is shown in Fig. 5. Because the perturbative solutions are normalized such that F (0) = G(0) = 1, which is not the case for the Airy functions, an appropriate normalization factor must be supplied: What is plotted in the dotted curve in the figure is R PT = − π 8κ 3 e 4κ 3 /3 . These curves reveal that the validity of the perturbative solution depends on the inequality It will be noted from Figs. 2 and 3 that the stress tensor components rapidly go to zero as one goes deeper into the potential, as expected. To further explore this, we look at the Green's function, which represents the expectation  value of the product of the electric fields in the medium, for the case α 0 = 1, This is plotted, for z = z ′ , in Fig. 6. For even larger z than shown in the figure, the reflection coefficient may be replaced by its small-κ expansion, and then the resulting analytic form of the diagonal Green's function ultimately agrees with that found by numerical integration. (For a twenty-term expansion ofR E , the error of the analytic approximation is less than 1% for z > 11.)

IX. REFLECTED POTENTIALS
Of course, there is no net force on the semi-infinite slab we have been considering to this point. This is because T zz must vanish at infinity, and once the obvious bulk subtraction is made, T zz (0−) = 0, according to Eq. (4.9). So suppose we consider two bodies, constructed by placing the mirror image of our potential to the left of z = 0: that is, we assume ε(z) = ε(−z). These are two bodies in contact, not disjoint. Then for either the TE or TM mode, the Green's function may be constructed in terms of the fundamental solutions of the homogeneous equations,F andG, whereF → 0 as z → +∞, andG → 0 as z → −∞, g(z, z ′ ) = 1 AF (z > )G(z < ). (9.1) in terms of the effective Wronskian factor A. If we expand this out in terms of the solutions on the right for the semi-infinite slab, denoted as previously by F and G, we find for z, z ′ > 0 where α is the Wronskian term for the half-space. Here the reflection coefficient is Perturbatively, it is easy to check that to first order which are twice as big as the values found for the semi-infinite slab, in Eqs. (6.6) and (6.15), as would be expected, because the slope discontinuity is doubled.
In the case of the plasma model, T s zz is finite, and for an exact linear potential was solved explicitly in Sec. VIII-see Fig. 2. So in the case of two facing reflected linear potentials in contact, one might think that a finite force of one body upon the other could be determined, where we have restored the proper scaling with the coupling. Although this appears to be a finite attraction between the two slabs, the interpretation of this is suspect for the reasons stated in Sec. II, because the body is not immersed in a homogeneous medium. As there is no distance scale in the problem aside from the coupling, it is impossible to connect this to a change in the energy according to the principle of virtual work. Moreover, the ambiguity of separating bulk and scattering parts remains.

X. CONCLUSION
In this paper we have extended our previous calculations on the "soft wall" problem to the electromagnetic case. In the plasma dispersion model, the TE mode coincides with the scalar case considered in Ref. [28]. Without dispersion we recover the universal edge behavior found by Ref. [24]. We also reproduce the Weyl divergences found in the scalar case. We do this, first by considering explicitly solvable examples, and then by performing a generic perturbative analysis for small slopes in the dielectric potential.
Let us summarize the salient features. For the plasma model, where the potential may be defined by ε(z) − 1 = v(z)/ζ 2 we see universal Weyl singularities in the bulk stress tensor for both TE and TM polarizations: which coincide with the divergences found for a scalar field [28]. (The second derivative term is seen for the quadratic potential treated in Appendix B.) For the nondispersive model, with temporal splitting, These results are very similar to those seen for the quadratic potential treated in Appendix B, with the replacements α/z → −β, α 2 /z 2 → (β 2 /4) ln z. One might think one could remove the Weyl divergences by removing all terms with polynomial growth in z, for surely such growth deep within the material is unphysical. Unfortunately, the WKB analysis of Ref. [28] shows there must also be z 2 ln z terms in the linear plasma-model TE case, which is confirmed by numerical experiments, so such a procedure appears impossible.
Although we recover expected results, as well as some new features, our analysis remains incomplete. It hinges on a break-up between bulk and scattering contributions, which is not unique; however, it captures the essential asymptotic behavior for large wavenumbers. The suggestion that to achieve a finite stress one merely omits the bulk terms is plausible, but this is not a unique process. Moreover, there are finite, position-dependent contributions to the stress tensor contained in the bulk term that likely cannot be merely discarded.
while for the spatial cutoff (which, without loss of generality we can choose to be in the x direction), In these expressions we haven't written the remaining function ofκ within the integrals. The relation between the two cutoff factors is just that given in Ref. [17]: Appendix B: Quadratic potential Suppose the potential begins quadratically, that is, it is continuous, with a continuous first derivative, but a discontinuous second derivative at the edge, ε(z) = 1 + βz 2 . (B1) We then easily find the fundamental solution to first order in β: Here we again note that the terms proportional to ζ 2 are identical. The Wronskians of the solutions are First consider the bulk divergences. The identity (6.21) still holds with the potential αζ 2 ǫz here replaced by βζ 2 z 2 , so it is straightforward to compute in the plasma model, where βζ 2 = β 0 is a constant, which is just as expected from the WKB analysis of Ref. [28]. The divergent terms are again the same for the corresponding TM contributions. And for the energy density, with temporal splitting again as expected. Although for the TM part a singularity emerges in the ζ integration once again, the first two terms here are reproduced. For the nondispersive, constant β, case we obtain results precisely analogous to those in Eqs. (6.5a), (6.14) and Eq.
For the scattering parts, we need the reflection coefficients: Then, for the normal-normal stress tensor and the energy density we obtain terms which are less singular toward the edge than was the case for the linear potential for the non-dispersive case: (B8b) Notice that the ratios of the zz components are 43/3, while the energy densities are in the ratio −9, exactly as in the linear case, which reflects the fact that the angular integrations over cos 2 θ = ζ 2 /κ 2 are the same.