A mathematical analysis of Casimir interactions I: The scalar field

Starting from the construction of the free quantum scalar field of mass $m\geq 0$ we give mathematically precise and rigorous versions of three different approaches to computing the Casimir forces between compact obstacles. We then prove that they are equivalent.

the Casimir energy to a determinant computed from boundary layer operators. Such determinant formulae result in finite quantities that do not require further regularisation and have been obtained and justified in the physics literature [2, 21-24, 41, 42, 49, 51, 52]. This has not only resulted in more efficient numerical algorithms but also in various asymptotic formulae for the Casimir forces for large and small separations. Many of the justifications and derivations of these formulae are based on physics considerations of macroscopic properties of matter or of van der Waals forces. As such they often involve ill defined path integrals. From a mathematical point of view considerations that link the determinant formulae to the spectral approach initially taken by Casimir have been largely formal computations or involve ad-hoc cut-offs and regularisation procedures. We note that in the Appendix of [42] it is proved correctly that the Fredholm determinant in the final formula of the Casimir energy is well defined. Only recently a mathematical justification of these formulae was given in [34], relating it to the trace of a linear combination of powers of the Laplace operator. The determinant formulae are directly related to the multi-reflection expansion of Balian and Duplantier [1] that also yields a finite Casimir energy. We mention [43] in the mathematical literature where the Casimir energy of a piston configuration is expressed in terms of the zeta regularised Fredholm determinant of the Dirichlet to Neumann operator. Summarising, we list several ways to compute the Casimir force acting on a compact object that have been proposed and carried out: (1) Using a total energy obtained in some way by regularising the spectrally defined zeta function. This can be done either directly or by first considering a compact problem in a box and then taking the adiabatic limit. (2) By integrating the renormalised stress energy tensor around any surface enclosing the object. The list is non-exhaustive and other methods exist, such as for example the worldline approach (see e.g. [28] and references). Here we will restrict ourselves to the listed methods.
The aim of the present paper is to establish, in the case of finitely many compact objects, the precise mathematical meaning of each of the listed methods for the case of the scalar field and then prove that they give the same answer for the force (not necessarily for the energy). The main tool to achieve this will be the relative stress energy tensor. This tensor mimics the definition of the relative trace of [34] and seems to not have been defined or studied previously in the literature. Note that the renormalised stress tensor becomes unbounded and non-integrable [10,15] near the boundaries of objects and this makes it unsuitable to compute the total energy component from the energy density T 00 . In contrast, the relative stress energy tensor is smooth up to the obstacles and is much more regular when considering boundary variations. This relative stress energy tensor does not satisfy Dirichlet or Neumann boundary conditions and therefore integration by parts involves boundary contributions. The relative energy density can be defined entirely in terms of functional calculus of the Laplace operator. This relative energy density has been introduced in [34]. It was shown to be integrable and its integral can be interpreted as the trace of a certain operator. The main result of [34] states that this trace can be expressed as the determinant of an operator constructed from boundary layer operators, thus providing a rigorous justification of the method (3) linking it with method (1). To show the equivalence of methods (2) and (3) we must provide a formula for the variation of the relative energy when one of the objects is moved. To compute this variation we prove and use a special case of the Hadamard variation formula ( [27,31,50]) adapted to the non-compact setting. We then show that as a consequence of this formula that the variation of the total energy equals the surface integral of the spatial components of the relative stress energy tensor (see Theorem 5.9). This surface integral is also equal to the surface integral over the renormalised tensor (see Theorem 5.10). We will now give a more precise formulation of our main result.  For the free scalar field of mass m ≥ 0 let T ren be the renormalised stress energy tensor as defined in Section 2.2. In QFT terms this stress energy tensor is equivalent to the usual stress energy tensor obtained from canonical quantisation when normal ordering is used with respect to the free vacuum state. This is a smooth symmetric two-tensor away from ∂O, but it is singular at ∂O and the integral of T 00 over E does not converge. Let T rel be the relative stress energy tensor as given in Definition 3.3. The relative total energy E rel is defined as the integral of (T rel ) 00 (x) over R d which can be shown to exist and to be equal the trace of a certain combination of operators. In case m > 0 the regularised energy E reg is defined in Section 6, Definition 6.3 via zeta function regularisation. We would like to compute the force on an object O j due to the presence of the other objects. Approach (2) is to directly compute where T ren is the renormalised stress energy tensor given in Definition 2.3 and Σ is any smooth surface enclosing O j (i.e. homologous to ∂O j in E), n k is the exterior unit normal vector field, and dσ is the surface measure on Σ. If energy conservation holds, i.e. no energy is radiated off to infinity, one expects this force to be the directional derivative of the total energy when O j is moved rigidly. Thus, let Z be a constant vector field on R d . Let Y be a vector field that equals Z near O j and vanishes near O l when l = j. The vector field generates a flow Φ that, for near zero, moves the object O j rigidly and we end up with a configuration that depends on the parameter . The total energies E rel and E reg then also depend on in that way and become functions of in a small interval around zero. The change of these energies with respect to the flow can be interpreted as the change of energy needed to move object O j rigidly relative to the other objects in the direction of Z. Our main result may be stated as follows.
Theorem 1.1. The relative energy E rel ( ) is differentiable in at = 0 and its derivative equals F i Z i , where Moreover, where Q rel (ω) = Q iω (Q iω ) −1 is defined in detail in Section 3 and constructed out of boundary layer operators for the Laplacian. If m > 0 then E rel ( ) − E reg ( ) is constant near = 0.
We note that this mathematical theorem simply shows that all these proposed computational methods give the same Casimir interactions in the case of separated rigid bodies. The statement does not say anything about the actual origin of the Casimir force or its existence, which needs to be determined from experiments or physics considerations. There is however a strong argument for the expression to be a directly relevant physical quantity. Our point of view is that the stress energy tensor does not have an absolute meaning in this context, but rather is used to compare two vacuum states (normal ordering depends on a comparison state). If we would like to know the effect for the rigid objects O 2 , . . . O N on the rigid object O 1 the states to compare are not the ground state with Dirichlet conditions and the free ground state. It is rather the vacuum states obtained from the Laplacian with Dirichlet conditions imposed on O 1 alone and with Dirichlet conditions on all the objects. The comparison of these two states yields a stress energy tensor that is completely regular near O 1 , and the computation of the force based on this tensor leads directly to the above formula without regularisation. The paper is organised as follows. In Section 2 we review the rigorous construction of the free scalar field of mass m ≥ 0 in the presence of boundaries and show how this leads to a natural definition of the renormalised stress energy tensor, which is given in Section 2.2. We also review its most important properties and express it in terms of spectral quantities for the Laplace operator. Section 3 introduces the relative setting and gives the definition of the relative stress energy tensor and its basic properties. Some norm estimates on the relative resolvent are given in Section 4, which provides mathematical justifications for later proofs. In Section 5 we prove a Hadamard variation formula and compute the variation of the relative energy to establish the first part of the main theorem. In Section 6 we show that for m > 0 the renormalised version of the zeta function has a meromorphic continuation and can thus be used to define the regularised energy. This section also contains a proof that variations of the regularised energy and the relative energy coincide. To illustrate the method and relate it to the classical computations we treat the easier example of the one-dimensional Casimir energy explicitly in Appendix B. This example also illustrates that a divergence term for the time-component of the renormalised stress energy tensor that is normally neglected needs to be taken into account to obtained the correct result (see Remark 2.5).
In a follow-up paper we will establish a similar theorem for the electromagnetic field. We note here that the stress energy for the electromagnetic field is quite different from the scalar field and there are additional complications such as zero modes ( [60]) that are absent for the scalar field. Moreover, the boundary conditions for the electromagnetic field are slightly more complicated, and cannot be reduced to Dirichlet boundary conditions. We therefore decided to not attempt a unified treatment which would obscure the result by additional notations. Our approach is expected to carry over to other boundary conditions such a Neumann, mixed Dirichlet-Neumann, or Robin boundary conditions with the single layer operators replaced by the appropriate layer operators. As in the electromagnetic case additional technical problems need to be overcome in these cases due to the possible appearance of zero modes and singularities of the Dirichlet-to-Neumann map at zero. In this paper the Schwartz kernel of A will be denoted byȂ.

Scalar quantum field theory with Dirichlet boundary conditions
Let −∆ be the (positive) Dirichlet Laplacian imposed on the codimension one submanifold ∂O. By definition this is the unbounded self-adjoint operator defined on the Hilbert space L 2 (R d ) associated with the Dirichlet quadratic form q D (f, f ) = ∫ X ∇f 2 dx with form-domain being the Sobolev space H 1 0 (X). As a consequence of elliptic regularity ([59, Section 7.2]) we have for the domain of dom((−∆) s 2 ) equipped with its graph norm for any s ∈ 2N 0 the continuous inclusions and each subspace is an invariant subspace for −∆ in the sense that any bounded function of the operator as defined by spectral calculus will leave these subspaces invariant. The restriction of −∆ to L 2 (O j ) is the Dirichlet Laplacian on the interior of O j and therefore has compact resolvent. The restriction of −∆ to L 2 (E) has purely absolutely continuous spectrum [0, ∞). By comparison we also have the free Laplacian −∆ free on R d which corresponds to the case O = ∅. Throughout the paper we fix a mass parameter m ≥ 0.
Definition 2.1. The relativistic Hamiltonian H is defined to be the self-adjoint operator The space-time M we consider is the Lorentzian spacetime R×X with Minkowski metric. The forward and backward fundamental solutions G + − ∶ C ∞ 0 (M ) → C ∞ (M ) of the Klein-Gordon operator ◻ + m 2 with Dirichlet boundary conditions are given by where θ is the Heaviside step function. As usual in canonical quantisation one considers the difference G = G + − G − given by Here Here the inclusion of the domain in the Sobolev spaces follows for s ∈ 2N 0 from elliptic regularity up to the boundary ( [48,Theorem 4.18]) and for general s ≥ 0 by interpolation. In particular, this means that H −1 sin(H(t − t ′ )) has a distributional integral kernel. We can define a symplectic structure σ on This induces a symplectic structure on GC ∞ 0 (M ) that is well known to coincide with the standard symplectic structure on the space of solutions. Indeed, if we define u = Gf and v = Gg for f, g ∈ C ∞ 0 (M ) then u and v solve the Klein-Gordon equation with Dirichlet boundary conditions and In this equality the right hand side is independent of t.
2.1. Field algebra and the vacuum state. The field algebra of the real Klein-Gordon field is then the (unbounded) CCR * -algebra of this symplectic space. Instead of using the symplectic space W one can describe this algebra A directly as the complex unital * -algebra generated by symbols Physical states of this quantum system are states on this * -algebra. The construction and physical interpretations of such states usually relies on a Fock representation of A. This representation is chosen on physical grounds as a positive energy representation. We briefly explain this now. The group of time translations (T s f )(t, x) = f (t − s, x) commutes with ◻ + m 2 and G and therefore α t (φ(f )) ∶= φ(T t f ) defines a group of *automorphisms of A. If a state ω ∶ A → C is invariant then this group lifts to a group of unitary transformations U (t) on the GNS-Hilbert space which is uniquely determined by We say that π is a positive energy representation if this group is strongly continuous and its infinitesimal generator has non-negative spectrum. We will focus in this paper on the quasi-free ground state. This means that the state is completely determined by its two point distribution ) which is given explicitly as i.e. ω 2 is the integral kernel of the operator H −1 e −iH(t−t ′ ) . In case m = 0 the spectrum of H contains zero and this expression needs to be interpreted in the sense of quadratic forms. Namely, it follows from general resolvent expansions (for example [60, Theorem 1.5, 1.6 and 1.7]) that C ∞ 0 (X) is contained in the domain of the operator H − 1 2 . This follows from the formula In particular, the space of test functions C ∞ 0 (X) is contained in the form domain of H −1 and therefore, by the Schwartz kernel theorem, the operator H −1 has a distributional kernel in D ′ (X × X). This is of course also the case for general m > 0. We will denote the integral kernel of H −1 byH −1 , mildly abusing notation. One can check directly that ω 2 defines a positive energy representation. Instead of using −∆ we could also have used −∆ free , the free Laplace operator. This also defines a positive energy state on the free algebra of observables A free which we will denote by ω free , and similarly we use the notation H free = (−∆ free + m 2 ) 1 2 . There states can be compared by restricting them to certain subalgebras that are contained in both the algebra of observables and the free algebra of observables. For example if U is contained in a double cone in M then A(U ) which is generated by φ(f ), supp(f ) ⊂ U can be thought of as a subset of both A and A free and therefore both states can be restricted to this algebra.

2.2.
The renormalised stress-energy tensor. The classical stress energy tensor of the Klein-Gordon field for a smooth real-valued solution u is given by This is a symmetric 2-tensor on M and one can easily show, using the Klein-Gordon equation, that it is divergence-free. Here g is the Minkowski metric on M with signature (−1, 1, . . . , 1). The Euclidean metric on R d will be denoted by h. The components T cl ij of the stress energy tensor are the restrictions to the diagonal of the functions Q cl ij (x, y) defined on M × M by The quantum field theory counterpart of Q can be written in the field algebra as a field-algebra-valued bilinear form in the test functions f 1 , f 2 ∈ C ∞ (M ) as The expectation value of Q ij with respect to the state ω is then given in terms of the two point function ω 2 as Let K(t, x, y) be the distributional kernel of 1 2 H −1 e −iHt . Then the distribution ω(Q ij ) with respect to the ground state is given by The above expressions are formal and make sense only when paired with test functions. We will use such formal notation throughout the paper when there is danger of confusion. The expectation value of the stress energy tensor would correspond to the restriction of ω(Q ij ) to the diagonal as a distribution. Unfortunately, the distribution ω(Q ij ) is singular and cannot be restricted to the diagonal in a straightforward manner. If one is interested in relative quantities only then one can still define the renormalised expectation value of the stress energy tensor between the states. Both states ω and ω free are positive energy states and therefore satisfy the Hadamard condition (for example [61,Theorem 6.3]). By uniqueness of such Hadamard states the difference of the twopoint distributions is smooth near the diagonal in M × M . In the present case this can also be seen more directly as follows. Let K free (t, x, y) be the kernel of 1 2 H −1 free e −iH free t . We will consider the difference K(t, x, y) = K(t, x, y) − K free (t, x, y).
Proof. The distribution K(t, x, y) is a solution of the wave equation By [34,Lemma 5.1] integral kernel κ λ of the resolvent difference is smooth and satisfies on any compact subset U of X × X a C k -norm bound of the form for some ≥ 0 and δ ′ > 0. Therefore the integral representation Since the initial conditions are smooth the solution K is smooth where it is uniquely determined by the initial data. This is the case in a neighborhood of t = 0 in R × X × X.
Definition 2.3. The components (T ren ) ij of the renormalised stress energy tensor T ren are defined to be the restriction to the diagonal of the function ω(Q ij ) − ω free (Q ij ).
Theorem 2.4. The renormalised stress energy tensor is symmetric and is given by and Note that hereH −1 andH −1 free are the the integral kernels of H −1 and H −1 free respectively. Moreover, the expressionȂ ∆ means the restriction of the integral kernel, A, to the diagonal (See Section 1.1).
Remark 2.5. The terms of divergence forms in the renormalised stress energy tensor are commonly neglected in the literature, as one may naively think that they have zero contribution when integrating over the whole space. However, this is not the case. As it is not integrable due to the singular behaviour near the boundary, the divergence theorem does not apply in this case. See Appendix B for the simplest case. The problem disappears when we work with the relative stress energy tensor given in Definition 3.3. where is the kernel of a symmetric operator on X with respect to the real inner product (⋅, ⋅) and therefore satisfies This implies Using product rules, we have which gives Applying equation (8) to (10), we have From equations (8) and (9) we have In other words, (T ren ) 0i = (T ren ) i0 = 0 for 1 ≤ i ≤ d. Hence equations (6), (11) and (12) show that (T ren ) ij is symmetric tensor on M . Moreover, When i = 0, which yields the expressions for the renormalised stress energy tensor.
Theorem 2.6. The renormalised stress energy tensor is divergence-free and independent of t.
Proof. Let K(t, x, y) and K (t, t ′ , x, y) be the same as the in the previous theorem. Recall that the shorthand expression of (3) is given by Then we have Now we use product rules to get Hence, one has We use the same trick to show that Now use the symmetric properties of K in equations (7) and (8), we have The function K solves the Klein-Gordon equation, i.e.
which shows the divergence-free property of the renormalised stress energy tensor.

The relative trace-formula and the Casimir energy
As mentioned in the introduction the renormalised stress-energy tensor T ren (x) becomes unbounded and non-integrable when x approaches the boundary of obstacles [4,10,15,25]. This prevents us from defining a renormalised total energy. One way to circumvent the problem is to introduce the relative framework of [34] which we now summarise. The main advantage of this construction is that it completely avoids ill-defined quantities and the need for regularisation.
More generally one defines the relative operator associated with a polynomially bounded Since all our operators are densely defined operators on the same Hilbert space L 2 (R d ) this combination makes sense. As a consequence of f being polynomially bounded the space C ∞ 0 (X) is in the domain of D f and therefore D f has a distribution kernel in D ′ (X × X).
To simplify our analysis later, we absorb the dependence of mass m in the functional f . We could write f m to emphasise the dependence on m, but the later analysis will not be affected by m. Therefore, we omit the m-dependence and have The main result of [34] is that for a large class of functions f , including the functions f (λ) = √ λ 2 and f (λ) = √ λ 2 + m 2 which are of interest in our context, the operator D f is trace-class and its trace can be computed by integrating the kernel on the diagonal. We now explain the precise statement of this result and its relation to the determinant of the boundary layer operator. In the following we will denote by λ and G free,λ =Ȓ free,λ are defined in an analogous way. By elliptic regularity these Green's distributions are smooth away from the diagonal x = y. Recall that for λ ∈ C we have the single layer potential operator be the Sobolev trace operator for s > 0. Properties of the Sobolev trace operator can be found, for instance, in [48]. One can write the above also as S λ = G free,λ ○ γ * . Restriction of S λ to the boundary defines an operator The operator Q λ is known to have the following properties. Since the boundary ∂O consists of N connected components ∂O j , we therefore have an orthogonal decomposition L 2 (∂O) = ⊕ N j=1 L 2 (∂O j ). Let p j ∶ L 2 (∂O) → L 2 (∂O j ) be the corresponding orthogonal projection. Now we havẽ In short,Q λ and T λ are respectively the diagonal and off diagonal part of Q λ . Now let D ν to be a sector in the upper half plane and it is given by where it suffices to consider 0 < < π 2 for our applications.
which is holomorphic in the upper half space and for some δ ′ > 0 satisfies the bound ).
See Theorem 1.6 of [34] for the above bound in the sector of the form {z Im(z) ≥ b Re(z) } for some δ > 0. Assume 0 < θ < π 2 and let S θ be the open sector Let P θ be the set of functions that are holomorphic and polynomially bounded in S θ .
Definition 3.1. The spaceP θ is defined to be the space of functions f such that f (λ) = g(λ 2 ) for some g ∈ P θ+ for some > 0 and there exists a > 0 such that We then have the following theorem.
In particular, choosing the function f (λ) = √ λ 2 + m 2 − m one obtains D f = H rel and therefore H rel is trace-class with trace equal to This follows immediately by deforming the contour integral using the exponential decay of Ξ in the upper half plane, considering the branch cut of Definition 3.3. The relative stress energy tensor is the renormalised stress energy tensor in the relative setting and it is defined as where (T ren ) O is the renormalised stress energy tensor for obstacle O and (T ren ) O i is the renormalised stress energy tensor for obstacle O i , which are defined at the beginning of this section.
Remark 3.4. One can also consider other versions of a relative stress energy tensor.
dropping the connectedness requirement of obstacles and work with The corresponding energy encodes the amount of work needed to separate the two obstacle configurations A and B. This quantity can also be expressed in terms of T rel for A and B. It is easy to see that this equals and therefore working with T rel only does not result in a loss of generality.
Theorem 3.5. T rel is smooth on X and extends smoothly to E as well as to O. The function (T rel ) 00 is integrable on R d .
Proof. By Theorem 5 we have The theorem was shown in [34] for the part 1 2H rel ∆ and the same method of proof can also be applied to the second term. We provide the full details here for the sake of completeness. We use two estimates proved in [34] which we now explain. Recall that the relative resolvent is given by For the integral kernelȒ rel,λ we write G rel,λ . As shown in [34] in the proof of Theorem 1.5 the integral kernel G rel,λ of R rel,λ is smooth up to the boundary on E as well as to O and its C k -norms on compact subsets K ⊂ E × E satisfy the bound for some ≥ 0 for λ in the sector containing the imaginary axis. We consider the operator  (21) and (22)]. Let λ ∈ D ν (See (15), i.e. a sector in the upper half plane) and k O,λ (x, y) denote the integral of R O,λ − R free,λ , then we have This shows , for some positive C 1 and C 2 . By Corollary 2.8 of [34] we have ).
Now we can conclude that and Let (17) and (18) to get the decay rate of h rel (x) by integrating over λ. That is, for ρ(x) > 1, h rel (x) has a decay of ρ −k (x)e −Cmρ(x) with both k and C being positive. This warrants the integrability of ∆h rel (x) for d ≥ 2, m ≠ 0 and ρ(x) > 1. Therefore, we will now focus on the case m = 0. By integrability of the integrand we can interchange differentiation and integration and therefore get Again, by integrating over λ, we have Let Ω ⊂ R d be an open set with dist(Ω, ∂O) > 1 and ϕ ∈ C ∞ 0 (R d ) satisfy 0 ≤ ϕ ≤ 1 and ϕ = 1 in a neighbourhood of R d Ω. Then we have the decomposition The integrability of ∆(ϕh rel ) in equation (20) follows from the smoothness property of the kernel of (H −1 ) rel at the diagonal as shown above. Thereby the integrability of ∆h rel on R d is equivalent to the one of ∆[(1 − ϕ)h rel ] on supp(1 − ϕ). This follows immediately from equation (19). Therefore, we have shown the integrability of ∆((H −1 ) rel ∆ ) on R d . Finally, the integrability of T rel follows from Theorem 2.4 and the definition of relative stress energy tensor (i.e. Definition 3.3).
Definition 3.6. The relative energy is defined as Theorem 3.7. We have the equality Proof. We have where B r is the ball of radius r centred at the origin. As ∆((H −1 ) rel ∆ ) has only jumptype discontinuity across ∂X, we can apply the divergence theorem to the integral on B r for sufficiently large r. That is , where (⋅) + and (⋅) − are the exterior and interior limits respectively. From equation (19), we also have which implies the contribution of the integral over ∂B r vanished as r → ∞ and therefore . From equation (16), we then have This shows We start by showing that for the restrictions to ∂O q we have the following identity To see this we temporarily denote by k(x, y) the integral kernel of H −1 O − H −1 Oq . This kernel vanishes O q × O q and the interior normal derivative therefore vanishes trivially.
We therefore only need to concern ourself with the exterior normal derivative. As shown in the proof of Theorem 3.5 the kernel (H −1 ) rel is smooth. One concludes from this, using Theorem 2.2, that The kernel k satisfies Dirichlet boundary conditions in both variables in the sense that k(x, y) = 0 if y ∈ ∂O q or if x ∈ ∂O q . By the chain rule y k(x, y)) y=x , which therefore vanishes on ∂O q . We then have Hence we have ∫ R d ∆((H −1 ) rel ∆ )dx = 0, which verifies the first equation in (21). The representation of the trace of H rel in terms of Ξ follows via Ξ ′ (λ) = −2λTr R rel,λ from Theorem 4.2 in [34], Theorem 3.2 and equation (16).

Estimates on the relative resolvent
In preparation for the proof of the variational formula we will need some additional estimates on the relative resolvent R rel,λ , which we collect in this section. We have the following well known layer potential representation (see for example [34,Equ. (19)]) which gives Also, let ρ be the function defined as For i, j, k ∈ Z, let ρ i,j;k be the functions defined as The following proposition partially follows from [ is a Hilbert-Schmidt operator whose Hilbert-Schmidt norm is bounded by Proof. Recall that we could also write the single layer potential operator as The second bound follows from the spectral representation of −∆ free . Therefore, we conclude the proof for part (1). Part (3) follows immediately from the bound on the Dirichlet-to-Neumann operator in [34]. For part (2) In particular, 0 < δ < dist(V, ∂O). Let k ∈ N. Now we have from [34] the estimate (1 + log (Im λ)) e −2δ Im λ for d = 4 (Im λ) d+2k−4 1 + Im λ 4−d e −2δ Im λ for d ≠ 2, 4 where ∆ = ∆ x + ∆ y is the Laplace operator on V × V ′ . In other words, we have Now by taking Sobolev trace, we have Statement (2) follows, using the properties of the Hilbert-Schmidt norm (for example Section A.3.1 in [59]).
where E s = ( √ −∆ ∂O + 1) s . To prove the last property of S λ , we are left with proving the bound for Note that the explicit kernel of φG free,λ φ, denoted by K λ (x, y), is given by is the Hankel function. By the Schur test and estimates on the free resolvent (see Appendix A for details), we have φG free,λ φ L 2 →H 2 ≤ C r 0 ,φ ρ 1,0;0 (Im λ) for d = 2, 1 for d ≠ 2.
Taking the adjoint, we get the same bound for H −2 → L 2 . Finally, using the estimate (24) with χ = 1 − φ, we obtain Remark 4.2. In the case of odd dimensions an easier argument can be used to provide a weaker estimate that is still sufficient for the purposes of this paper. Using the strong Huygens principle one deduces ([17, Section 3.1]) which implies This gives φ S λ H s →L 2 ≤ C φ,ν (1 + Im λ) for − 3 2 ≤ s ≤ − 1 2 . Applying inequality (24) for χ = 1 − φ and using −2s+1 (2) and (3) In this section, we will show that a version of Hadamard variation of the renormalised stress energy tensor (T ren ) ij defined in Section 2.2 is related to the Hadamard variation formula for the resolvent. We will follow the methods developed in [27,32,33,50] and derive a Hadamard variation formula for the relative resolvent, then apply it to the relative stress energy tensor. Short proofs of Hadamard variation formula can be found in [46,54] for the case of bounded domains. Since we are dealing with an unbounded domain we extend theory to non-compact setting for the special case of boundary translation flows (see Definition 5.1) in Theorem 5.3, which are sufficient for our purposes. Following Peetre's derivation of Hadamard variation formula, we define the following variational derivative.

Hadamard variation formula and equivalence of approaches
Definition 5.2. Let u be a (weak- * ) C 1 curve of functions in D ′ (U ). The variational derivative at = 0 is given by where θ Y u = lim →0 ϕ * u −u 0 is the variational derivative defined by Garabedian-Schiffer's in [27]. Here the derivative θ Y u is understood in the weak- * -sense and the action of the vector-field Y is understood in the sense of distribution.
Note that θ Y u is different from the standard (conventional) Lie derivative, as u may have an additional dependence on . In fact, the last term, Y u in equation (26), should be understood as the conventional Lie derivative of u 0 . The derivation of Hadamard variational formula for the resolvent associated with the Dirichlet Laplace operator usually starts with the energy quadratic form (see [27,33,50]). The energy quadratic form associated with the Dirichlet Laplacian −∆ U on U is given by where λ ∈ C and Im λ > 0. Using the diffeomorphism flow ϕ , one can pull-back the quadratic form from U to U 0 = U , which gives a one-parameter family of quadratic forms on U , i.e.
Note that the energy form (27) and the induced energy forms (28) are related by The operator associated with the energy form (27) is the Dirichlet Laplace operator, whereasẼ defines a one-parameter families of elliptic operators on U for sufficiently small . Let G λ, = G U ,λ be kernel of the resolvent for the Dirichlet Laplacian on U . Then from equations (27), (28) and (29), we have in the sense of distributions where y is in the interior of U and by elliptic regularity G λ, (x, y) is then smooth at the boundary and therefore it makes sense to define its boundary value.
As we would like to study the variation of resolvents, it is convenient to consider the variation as distributions on U ×U . In other words, for R ∈ D ′ (U ×U ) and u, v ∈ C ∞ 0 (U ), we have, from the Schwartz kernel theorem, where the first two brackets correspond to the pairing between distributions and test functions while the third and the forth brackets are the L 2 inner products on U and U respectively. It is not hard to see that the existence of the variational derivative of R in the sense of (26) is implied by the existence of θ Y R. From equation (31), the existence of θ Y R = lim →0 ϕ * R −R 0 in the weak- * -sense is equivalent to the existence of the standard derivative of r( ) = ⟨u, ϕ * R ϕ * − v⟩ with respect to for all u, v ∈ C ∞ 0 (U ). The following theorem is the well known Hadamard variation formula for the resolvent extended to case of unbounded domains in our setting. Theorem 5.3. Let ϕ be a boundary translation flow, then the variational derivative of R λ, exists in the weak- * topology. Let G λ, (x, y) be the kernel of R λ, , then its variational derivative is given by Proof. Firstly, we prove the existence of θ Y R λ in weak- * -sense. We know that, from equation (23), ϕ * R λ, ϕ * − − R free,λ = −ϕ * S λ, Q −1 λ, S t λ, ϕ * − . We will therefore establish differentiability of for any fixed test functions g, f ∈ C ∞ 0 (U ) and compute its derivative. In the last term of equation (33), the operator S t λ, is the transpose operator to S λ, obtained from the real inner product, i.e. S t λ, f = S t λ, f . Since the free resolvent is smooth off the diagonal and the kernel of S λ, is smooth on U × ∂U . Here f ∈ C ∞ (∂U ), g ∈ C ∞ 0 (U ) and G free,λ is the kernel of the free resolvent. To establish differentiability by the product rule it is sufficient to prove the existence of d d ϕ * S λ, ϕ * − in the C ∞ -topology of integral kernels on U × ∂U , and the existence of d d ϕ * Q −1 λ, ϕ * − at = 0 in the weak- * -sense. The free resolvent kernel is smooth off the diagonal and therefore the above formula for ϕ * S λ, ϕ * − shows differentiability of the smooth kernel in the parameter at = 0 and the classical sense. We are thus left with proving the existence of d d ϕ * Q −1 λ, ϕ * − at = 0 in the weak- * -sense. Note that ϕ * Q −1 λ, ϕ * − is a one-parameter family of maps from C ∞ (∂U ) to C ∞ (∂U ), i.e. the spaces do not depend on . Similar to equation (14) we have the splittingQ where p j are the orthogonal projections L 2 (∂U ) → L 2 (∂U j ) and ∂U j , j = 1, . . . , N are the connected components of ∂U . Define Q λ, = ϕ * Q λ, ϕ * − ,Q λ, = ϕ * Q λ, ϕ * − and T λ, = ϕ * T λ, ϕ * − . By the definition ofQ λ, , we havẽ where f ∈ C ∞ (∂U ) and f i = p i f . As ϕ is a boundary translation flow, we obtain the following relationships.Q λ, =Q λ,0 =Q λ andQ −1 λ, =Q −1 λ . Now, from the decomposition of Q λ in equation (14), one obtains The family T λ, is a differentiable family of smoothing operators for sufficiently small (i.e. no obstacles are overlapping) and its derivative in therefore, by Taylor's remainder estimate, exits as a family of smooth kernels. Hence, Q λ, is differentiable in at = 0 as a family of operators from H s (∂U ) to H s+1 (∂U ) for any s ∈ R. We now use that Q λ,0 is invertible and the inverse Q −1 λ,0 is a pseudodifferential operator of order one, and maps H s (∂U ) to H s−1 (∂U ) . Since the space of invertible operators is open the inverses Q −1 λ, exist near = 0 as maps from H s (∂U ) to H s−1 (∂U ). Hence, Q −1 λ, is differentiable in at = 0 as a family of operators from H s (∂U ) to H s−1 (∂U ). In particular the derivatives d d Q λ, and d d Q −1 λ, exist in the weak- * -sense. Hence the variational derivative of R λ, exists in the weak- * sense and it is given by where Y ′ means the action of Y on the second variable. It remains to compute the derivative. To do this we consider the inhomogeneous problem Let ν be the exterior unit normal of U , y be an interior point of U such that ϕ (y) =ỹ and e(u, v) = ∇u ⋅ ∇v − λ 2 uv. By applying to equation (35), taking derivative in of equation (35) and using equation (30) and Peetre's computations [50], one has Let Using the symmetric property of G λ, (x, y) = G λ, (y, x) and equations (34), (36) and (37), we obtain Using the boundary conditions G λ, (ϕ (x), y) = 0 for x ∈ ∂U , we arrive at the Hadamard variation formula for the Dirichlet resolvent.

Application of the Hadamard variation formula to the relative resolvent.
We now apply the Hadamard variation formula to our setting with finitely many obstacles, combining the above formulae for U = E and U = O. We have from equation where (⋅) + means taking limits from E to the boundary ∂E and (⋅) − means taking limits from O to the boundary ∂O. We will now use the variational formula for the relative resolvent to prove the following theorem. Hence, we define O = ϕ (O) and similarly O j, = ϕ (O j ). In this way we can define the relative resolvent and its integral kernel G rel,λ, depending on the parameter .
Theorem 5.4. Let ϕ be a boundary translation flow, d ≥ 2 and λ ∈ D ν . Then R rel,λ, is for each λ ∈ D ν a C 1 trace-class operator valued function of near the point = 0. Its derivativeṘ rel,λ,0 equals δ Y R rel,λ and there exists δ > 0 such that L 2 -trace-norm oḟ R rel,λ,0 is bounded by Its kernel,Ġ rel,λ,0 =Ȓ rel,λ,0 , is given bẏ Proof. We start by showing that the family is Fréchet differentiable in the Banach space of trace-class operators with continuous derivative, i.e. the function is C 1 as a trace-class operator valued function on I, where is a fixed sufficiently small compact interval. As in the previous section we have We can decompose its variation form as in the proof of Theorem 5.3 Then we split the last term into a product of three terms, i.e. S λ, ϕ * − , ϕ * (Q −1 λ, −Q −1 λ, )ϕ * − , and ϕ * S t λ, . The first operator S λ, ϕ * − is given by Since ∂O is a disjoint union of the components ∂O j the operators S λ, ϕ * − splits into a sum Here we used the fact that Z is constant and equal to Z j near ∂O j and that the free Green's function is translation invariant. Since T j ( ) is C 1 as a family of maps H 1 (R d ) → L 2 (R d ) this shows that S λ, ϕ * − and its adjoint are C 1 as families of bounded operators L 2 (∂O) → L 2 (R d ). As shown in the proof of Theorem 5.3, the operator ϕ * Q λ, ϕ * − is independent of and therefore ϕ * Q λ, ϕ * − =Q λ + ϕ * T λ, ϕ * − =Q λ + T λ, . The map T λ, has smooth integral kernel that depends smoothly on for sufficiently small . This family is therefore C 1 as a family of trace-class operators. We temporarily denote ϕ * Q λ, ϕ * − by J , As G free,λ is translation invariant, ϕ * Q λ, ϕ * − is independent of and hence we denote it byJ. Then, by the above J = J 0 + r and J 0 −J is trace-class. Moreover, the remainder term r is of the form r =J 0 + ρ( ), where ρ( ) 1 = o( ) as → 0. By the Neumann series we have We have used here that trace-class operators form an ideal in the algebra of bounded operators, and the norm estimate AB 1 ≤ A B 1 holds. We conclude that ϕ * (Q −1 λ, −Q −1 λ, )ϕ * − is C 1 as a family of trace-class operators. We now compute the derivative in of all the three terms. We obtaiṅ For the derivative of the second term, we have as in Theorem 5.3 Therefore, the variation of the relative resolvent is given bẏ To estimate the trace-norm ofṘ rel,λ,0 , note that the first term in the above equation can be estimated by The third term in equation (41) can be bounded the same as the first term. For the second term one can estimate it by Combining equations (42) and (43), one has Ṙ rel,λ,0 1 ≤ C(1 + (Im λ) 2 ) 4 ρ(Im λ) E 1 T λ 1 + E 1Ṫλ,0 1 .
As in Theorem 5.3, T λ, is a smoothing operator that depends smoothly on , as long as we have dist(∂O j, , ∂O k, ) > 0 for all pairs of obstacles. Since the obstacles are compact, T λ, is a smoothing operator on compact domains and hence it is also a trace-class operator from H s (∂O) → H s (∂O) for all s ∈ R. This proves the first part of the theorem. Also, from Theorem 5.3, we know that δ Y G rel,λ,0 exists in the weak- * sense. By equations (33) and (34), we know that the kernel ofṘ rel,λ,0 coincides with δ Y G rel,λ,0 . In other words, the variational derivative exists in a stronger sense. We can therefore apply the variation formula (38) to the relative resolvent, which giveṡ Since the interior parts of G O,λ, and G O j ,λ, are the same, the interior contributions of δ Y G O,λ,0 in the expression (38) cancel out with the ones of δ Y G O j ,λ,0 . Therefore, we are left with only the exterior contributions as shown in equation (40). Our next proposition gives a relationship between the trace of the variation of the relative resolvent with a local trace on the boundary.
Proposition 5.6. The trace of the variation of the relative resolvent on L 2 (X) is also given by Proof. Let γ E,ν ∶ H 3 2 (E) → L 2 (∂O) be the Sobolev trace after taking the exterior normal derivative (ν is the exterior normal vector field) and To see that this map is well defined we note that where S O,λ and Q O,λ are the same as S λ and Q λ defined in (22), but with emphasis on the dependence on O. Then R free,λ γ * E,ν maps L 2 (∂Ω) continuously to L 2 (R d ). The operator S t λ γ * E,ν is the double layer operator on the boundary ∂O and maps continuously L 2 (∂Ω) → L 2 (∂Ω) (see for example [12]). Since Q −1 λ is a pseudodifferential operator of order one, and S λ continuously maps H −1 (∂O) to H where Q O i ,λ is the operator Q λ in equation (22) when O is replaced by O i in the definition. As in (14) we have used the decomposition Q λ =Q λ + T λ . Since T λ is smoothing, Similarly, as the free Green's function is smooth off the diagonal the operator p k S t λ γ * E,ν p i has smooth integral kernel for k = i. In particular these operators are trace-class as maps from L 2 (∂O) to L 2 (∂O). Since S λ as well as S λQ −1 λ is bounded from L 2 (∂O) to L 2 (R d ) this shows that for every i ∈ {1, . . . , N } the operator where Y i = ⟨Y, ν E i ⟩ ∂O i is viewed as a multiplication operator acting on L 2 (∂O i ) ⊂ L 2 (∂O). Taking the trace, we have Here the cyclic permutation under the trace is justified because of the nuclearity of We then obtain By Theorem 3.2 one can define trace-class operators g rel = D f = i π ∫Γ λf (λ)R rel,λ dλ for f (λ) = g(λ 2 ) and g ∈ P θ . We now have the following.
Proposition 5.7. Let g ∈ P θ and ϕ be a boundary translation flow. Then δ Y g rel is a C 1 trace-class operator valued function of near the point = 0. Its derivativeġ rel = δ Y g rel satisfies where g ′ (z) = dg dz (z). Proof. By Theorem 5.4 the operator R rel (z) is, for fixed z ∉ S θ , in the Banach space C 1 (I, L 1 ) of trace-class operator valued C 1 -functions on a compact interval I containing zero. Differentiation defines a closed operator from on C(I, L 1 ) with domain C 1 (I, L 1 ). By Theorem 5.4 the derivative of R rel,λ is integrable in λ. An application of Hille's theorem to the Bochner integral defining g rel in the Banach space of trace-class operators shows that differentiation in commutes with integration. We therefore know that g rel is differentiable and Let f (λ) = g(λ 2 ). Using Proposition 5.6 and integration by parts in λ, we have In the special case g(z) = √ z + m 2 Proposition 5.7 shows differentiability of H rel with respect to in the space of trace-class operators at = 0. Using Theorem 3.7 and differentiating under the trace then gives the following theorem.
Theorem 5.8. The variation of the relative energy is given by δ Y E rel = 1 2 Tr(δ Y H rel ). We will now use the Hadamard variation formula to compute this variation.

Variation of the Klein-Gordon energy tensor.
Theorem 5.9. Let Y be a smooth boundary translation vector field. The variation of the Klein-Gordon energy generated by Y is equal to the boundary integral of its spatial tensor contracted with Y . That is where the integration on the right-hand-side is at the exterior boundary and ν is the exterior normal for E.
Proof. From equations (5), we have We know that • ((T ren ) Op ) ij is smooth on a neighbourhood of O q for p ≠ q (by Theorem 2.6), • ((T ren ) Op ) ij is divergence-free (by Theorem 2.6), • Y is constant on a neighbourhood of O q (by assumptions) for all q ∈ {1, . . . N }. Therefore, F i = ((T ren ) Op ) ij Y j is smooth and divergence-free on a neighbourhood of O q . In other words, we have As ((T ren ) O − (T ren ) Op ) ij is vanishing at the boundary ∂O p , we have Altogether, we have The second term in the above equation can be expressed as where we used the properties (11) and (30) Equations (45) and (46) imply Since ∫ ∂O (T rel ) ij ν i Y j dσ is integrating at the exterior boundary, we have Applying Proposition 5.7 to g(z) = √ z + m 2 , we have Op ⟨Y, ν⟩] − = 0 as Y is a boundary translation vector field. Therefore, we have That is Finally, equations (47), (48) and Theorem 5.8 complete the proof.
An application of an analogue of Theorem 5.9 to calculate the Casimir force in dimension one can be found in Appendix B. Finally, we have the following theorem.
Theorem 5.10. Let T ren be the renormalised stress energy tensor in Theorem 2.3 and let Y be a boundary translation flow as in Theorem 5.9. We assume further that Y is constant near O p for some p ∈ {1, . . . , N } and vanishes near O q if q = p. Let S be any smooth hypersurface in E that is homologous to ∂O p in E and let ν be the exterior normal vector field of S. Then the variation of the relative energy generated by Y is equal to where Z is the unique constant vector field on R d whose restriction to ∂O q equals Y .
Proof. As in the proof of Theorem 5.9, we have We also know that ∂O p is homologous to S R , a sphere with sufficiently large radius To get a decay property of (T ren ) Op at infinity, we first recall that For m ≠ 0, we would have an exponential decay of e −Cmρ(x) for H O −H free (x, x), as explained in the proof of Theorem 3.5. For m = 0, one could use the estimates of (17) and (18) to obtain that, for d ≥ 2, Moreover, estimate (19) also implies Applying the above estimates to Theorem 2.4, one concludes that This implies lim R→∞ ∫ S R ((T ren ) Op ) ij ν i ν j ⟨Y, ν⟩dσ = 0. As (T ren ) O = T ren (see Definition 2.3 and the first two paragraphs of Section 3), this completes the proof.
6. The zeta regularised energy and the equivalence of (1) and (3) In this section we assume that m > 0. In that case it is well known that for Re(s) > (d + 1) 2 the operator is trace-class (see for example [5]) and we can therefore define the renormalised zeta functions ζ O (s) as where ξ(λ) is the spectral shift function of the problem. The Birman-Krein formula applies to this setting and we have where S(λ) is the stationary scattering matrix of the problem. One can also use the Mellin tranform to write The following Lemma should be well known but we could not find a reference for the precise statement. It is a simple consequence of heat kernel expansions [3,26,30] and Kac's principle of not feeling the boundary [39,45]. We also refer to [35,40] for more details on obstacle scattering theory and the Birman-Krein formula.
Lemma 6.1. The function h O (t) is exponentially decaying as t → ∞ and has a full asymptotic expansion as t ↘ 0 of the form where the infinite sum is understood in the sense of asymptotic summation. The coefficients a k are integrals over ∂O of locally determined quantities expressed in terms of the extrinsic and intrinsic curvature of the boundary and its derivatives. In particular, Proof. The exponential decay follows immediately from the representation by means of the spectral shift function and m > 0. Now h O (t) is the trace of the difference of the two heat operators e −tH 2 O and e −tH 2 free with integral kernels K O (t, x, y) and K free (t, x, y), respectively. Since the difference is trace-class and the integral kernel is smooth we have where B R is a ball of radius R, i.e. integration is over a large ball of radius R with the obstacles removed. The heat kernel difference satisfies not feeling the boundary estimates. For example a general finite propagation speed estimate ( [45]) gives Let U ⊂ E be an open neighbourhood that contains ∂E, i.e. ∂E ⊂ U ⊂ E. Then we have for t < 1 that for any N > 0. This computation is therefore purely local.
where the b 2k are integrals of locally defined quantities over M , and the b 2k+1 are integrals of locally defined quantities over ∂M which are determined by the boundary conditions. When considering differences heat kernels of Laplace operators with different boundary conditions the b 2k terms cancel and only the terms containing b 2k+1 remain.
It follows as usual (for example [29,Section 1.12]) that ζ O has a meromorphic continuation to the complex plane. If d is odd then there are finitely many poles at . . , 1} with residue at d−1 2 − k determined by the coefficients a k . In this case the values at non-positive integers are also expressible in terms of a k . In case d is even poles may be located at the points { d−1 2 − k k ∈ N 0 }. Definition 6.3. The regularised energy E reg is then defined where FP z=a f (z) denotes the finite part of the meromorphic function f at the point a, i.e. the constant term in the Laurent expansion of f about the point a.
In particular, in case d is odd we have E reg = 1 2 ζ O (− 1 2 ). We can also define a zeta regularized energy E j reg for every object O j . Obviously, E j reg does not depend on the position of O j in R d and is also invariant under active rotations of the object. Since the heat coefficients are local quantities the relative zeta function is an entire function. Since the relative quantities are trace-class for s < 0 we also have that E rel = 1 2 ζ rel (− 1 2 ) and therefore Thus, E rel − E reg does not change if the individual objects are translated or rotated.

Proof of main theorem
In this section, we will prove our main theorem (Theorem 1.1) by combining the results we obtained in the previous sections.
For the operator χG free,λ φ, we obtain from equation (52) that which gives the verification of the estimate (24).
Appendix B. The method illustrated for the 1-D Casimir effect In this appendix, we illustrate Theorem 5.9 in its simplest form i.e. for the case of the 1-dimensional Casimir effect with m = 0. This will also illustrate the advantages of the relative framework. Let a 1 < b 1 < a 2 < b 2 , where O 1 = (a 1 , b 1 ) and O 2 = (a 2 , b 2 ) are the obstacles. Then we have x, y < a 2 G (a 2 ,b 2 ) a 2 < x, y < b 2 G (b 2 ,+∞) x, y > b 2 .
In particular, The same calculation yields for x, y < b. That is When restricting to the diagonal, we have Now for a < x, y < b, we have When restricting to the diagonal, we have Equation (53) gives From equations (54) and (55), we have This equation shows that H rel (x) is continuous, which is consistent with the claim in the proof of Theorem 5.8. Integrating over R, we have Similarly, one has the renormalised counterpart of H rel (x), which is given by (H O ) ren = (H O − H free ) ∆ . Note that this only corresponds to the first term in T 00 of (5), i.e.
It is easy to see that H ren is not integrable. Therefore, some regularisation schemes would be needed at this point. One way is by heat-kernel regularisation (see, for instance, [25]). However, this only resolves the non-integrability problem of the first term of T 00 . We also need to integrate the term 1 8 ∆[(H −1 − H −1 free ) ∆ ] in equation (5) over R. This is also ill-defined, as it is not integrable. We will see that these problems disappear when we work in the relative setting. Restricting equation (53) to the diagonal and then taking the action of Laplacian, we have Integrating spectral variable k alongΓ and then over the space variable x, we have hence E rel = 1 2 Tr R (H rel ) + .
Note that using heat-kernel regularisation, one would also obtain E rel , see [25]. Equations (56) and (58) agree with Theorem 3.7. Note that equation (57) also shows that where all the three terms are ill-defined as they are not integrable. For instance, (a 1 −x) 2 singularity when approaching x = a 1 . This justifies Remark 2.5. Now let X be a smooth vector field that generates a movement of (right) obstacle 2 to the right with a constant speed v. Moreover, X is zero around (left) obstacle 1. In other words, we move the obstacle 2 to the right by v and keep the obstacle 1 stationary. Now the variation of the relative energy is given by δ X E = v ⋅ ∂ a 2 E. The left hand side of equation in Theorem 5.8 becomes Now the identity (48) used in the proof of Theorem 5.9 says It becomes Note that Combining equations (58), (60) and (61), we have verified the identity (48) in one dimensional cases. Moreover, equations (59), (60) and (61) are consistent with Theorem 5.8 and Theorem 5.9.