1 Introduction

The Casimir effect, in most general terms, may be defined as the quantum backreaction response to adiabatic changes in external macroscopic conditions. The effect was first predicted in the work of Hendrik Casimir from 1948 [3], in which the considered physical system was the quantum electromagnetic field, and the external conditions were created by the presence of parallel conducting plates; the backreaction effect is then a force exerted on the plates, tending to change their separation distance. More generally, one can consider any quantum system subject to external conditions characterized by a classical parameter; the tendency to change this parameter is then a generalized backreaction force.

Since that time the scope of the investigated physical situations in which such phenomena appear, both on the theoretical and on the experimental side, expanded enormously, see for example the book [4]. Regardless of this immense research many points remain unclear and many controversies are open. We mention only a few, that we judge as most relevant to the Casimir effect at large and to the content of this work.

With many experiments supporting theoretical predictions, there is little doubt that this effect does exist, but the relation of formulas most often appearing in theoretical works to specific experiments results is a much more nuanced question. There are many reasons for this situation, we mention only a few of them that appear in testing the Casimir force for two conducting plates. First, measurements are not conduced in absolute zero temperature, which on the theoretical side brings the need to model the thermal excitation of the conducting plates structure, see [12]. Second, real conducting plates are not ideally reflective for electromagnetic field, which means that the field does not obey strict boundary conditions. Consequently, this idealized setting should be replaced by a more realistic model of the field-plates interaction. Last, two parallel plates constitute a rather challenging setting for experiments, so they are often replaced by a configuration of a flat plate and part of a sphere. To adapt theoretical predictions to such setting one needs some approximation scheme (like e.g. the so called proximity force approximation). For more information and reference we direct the reader to [13].

We need to mention that on the purely theoretical side even the zero temperature Casimir effect is a subject of many problems and controversies, which at present day look far from being fully resolved. There are a variety of approaches to this phenomena, the most popular being the original Casimir zero-point energy method (see for example a review article [2]). Other techniques include the Green function method (see again [2]) and the scattering theory approach (see [6]). The relation between all these methods and resulting formulas is not quite transparent at this moment.

In this work we adopt the approach to the Casimir effect developed by Herdegen [10], motivated by algebraic methods in quantum field theory. We use this method to investigate a class of models describing two macroscopic spherically symmetric disjoint objects interacting with scalar field. We also consider appropriate limit in which both bodies are shrunk to points, resulting in two Dirac delta potentials. The use of the aforementioned framework allows us to make clear statements about the nature of our calculations, to derive rigorous results and make our assumptions about the class of the considered physical situations transparent. This becomes important when we compare our results to those attained by Scardicchio [17], Spreafico and Zerbini [18] and partially by Fermi and Pizzocchero [9] and Fermi [7]. We postpone this discussion to Sect. 9.

We also point to the fact, that within the algebraic field theory approach an investigation of the Casimir effect was also carried out by Dappiaggi, Nosari and Pianamonti, who considered the system of one or two parallel metal plates [5]. Their approach uses a different branch of algebraic quantum theory, the so-called functional formalism and as such is explicitly described by authors as “parallel and complementary” to that used here.

1.1 Sketch of the Framework and the Main Result

We briefly outline here physical and mathematical basis of the framework to be used in this article, for details see the original work [10]. We share the standard view that the generalized Casimir backreaction force should be obtained as a derivative of the global energy, which depends on the quantum state of the system and on some classical macroscopic parameters. But even at this point views on the nature of “Casimir energy” vary significantly. Here this will be the expectation value of one fixed operator H in a well defined ground state, which removes all possible ambiguity in its definition. Second, we suppose that macroscopic bodies are “heavy” enough, so a change in the classical parameters a, say the distance between the plates, is adiabatic, as is also the change of quantum ground state.

We now express statements of the last paragraph in a more precise way. Let the observables of the “free” system, say free quantum field, be represented in a Hilbert space \(\mathcal {H}\), and the dynamics be defined by the Hamiltonian H. Now one wants to introduce some external macroscopic bodies into the system, in our case two spherically symmetric bodies, whose presence will change the quantum dynamics. However, in order not to destroy the original quantum system the change cannot be to drastic, as to enforce a change of the representation of the basic quantum observables of the system. Thus we have to put severe restriction on an acceptable model of macroscopic perturbation: there is no change of Hilbert space \(\mathcal {H}\) in which our observables act. Once this is fulfilled, let the modified dynamics be given by the Hamiltonian \(H_{ \vec { a }}\), where \(\vec { a }\) is a parameter characterizing the macroscopic model. We use the vector notation \(\vec { a }\), since in our case this parameter is the vector connecting centers of the two considered bodies. We assume that the ground state of \(H_{ \vec { a }}\) is unique and we denote its vector \(\Omega _{ \vec { a }}\). Then the Casimir energy in this state is

$$\begin{aligned} \mathcal {E}_{ a } = ( \Omega _{ \vec { a }}, H \Omega _{ \vec { a }} ), \end{aligned}$$
(1)

and the Casimir generalized force is given by

$$\begin{aligned} \vec { \mathcal {F} }_{ \vec { a }} = -\frac{ \mathrm{d}^{ } { \mathcal {E}_{ a } } }{ \mathrm{d} { \vec { a } }^{ } }. \end{aligned}$$
(2)

In other words, we defined “Casimir energy” as the energy of the free field in the ground state of the system induced by macroscopic objects. Since due to symmetry the energy can depend only on \(a = |{ \, \vec { a }\, }|\), we omit the vector “arrow” in its subscript. The same convention applies to the number of particles to be introduced below.

Within this setting, our quantum system will be the massless scalar quantum field \(\phi ( t, \vec { x })\). The free field dynamics, and the dynamics modified by the presence of the centers, are governed, respectively, by the equations

$$\begin{aligned} \left( \frac{ \partial ^{ 2 } { } }{ \partial { t }^{ 2 } } + h^{ 2 } \right) \phi = 0, \qquad \left( \frac{ \partial ^{ 2 } { } }{ \partial { t }^{ 2 } } + h_{ \vec { a } }^{ \;\; 2 } \right) \phi = 0, \end{aligned}$$
(3)

where \(h^{ 2 } = -\Delta \) and \(h_{ \vec { a } }^{ \;\; 2 } = h^{ 2 } + V_{ \vec { a } }\), with \(V_{ \vec { a } }\) representing the external perturbation, such that \(h_{ \vec { a } }^{ \;\; 2 }\) is positive. To satisfy the algebraic demands mentioned above the square root operators h and \(h_{ \vec { a } }\) cannot differ to much, and the precise condition is that the operator \(h_{ \vec { a } }^{ \;\; -\frac{ 1 }{ 2 } }( h_{ \vec { a } } - h ) h^{ -\frac{ 1 }{ 2 } }\) has to extend to a Hilbert–Schmidt operator [10]. If this property is satisfied, the system is represented in the Fock space of the free field and as was shown in [10], the expectation values of the particle number \(\mathcal {N}_{ a }\) and of energy \(\mathcal {E}_{ a }\) are given by

$$\begin{aligned} \mathcal {E}_{ a }&= \mathrm {Tr}\left[ ( h_{ \vec { a }} - h ) h_{ \vec { a }}^{ -1 } ( h_{ \vec { a }} - h ) \right] , \end{aligned}$$
(4a)
$$\begin{aligned} \mathcal {N}_{ a }&= \mathrm {Tr}\left[ h^{ -1/2 } ( h_{ \vec { a }} - h ) h_{ \vec { a }}^{ -1 } ( h_{ \vec { a }} - h ) h^{ -1/2 } \right] . \end{aligned}$$
(4b)

The finiteness of \(\mathcal {N}_{ a }\) guaranties that the algebraic structure in the case of free field and of the field with macroscopic perturbation is the same, while finiteness of \(\mathcal {E}_{ a }\) is an additional physical requirement on a realistic system.

We decide to follow the work by Herdegen and Stopa [11] in modeling macroscopic bodies, so we choose \(V_{ \vec { a }}\) to be a finite rank operator, which we shall call a quasi-potential. This operator is nonlocal, but its nonlocality is under strict control, and the model is admissible within the adopted framework. We also investigate a family of rescaled quasi-potentials \(V_{ \vec { a },\, \uplambda }\), \(\uplambda \in ( 0, 1 ]\) and show that in the limit \(\uplambda \searrow 0\) our models tend to the model of a scalar field interacting with two Dirac deltas \(\delta ( \vec { x })\). This situation is in many ways similar to that considered in the article mentioned above, where the authors investigate canonical example of two parallel plates with Dirichlet and Neumann conditions. While strict two \(\delta ( \vec { x })\) are not admissible, since they violate both requirements (4a) and (4b), for all \(\uplambda \in ( 0, 1 ]\) the problem is well defined and we obtain the asymptotic expansion of the Casimir energy

$$\begin{aligned} \mathcal {E}( a, \uplambda )= & {} \mathcal {E}_{ \text {self} }( \uplambda ) + \frac{ 2 \alpha }{ \pi ^{ 3 } } \Bigg [ \frac{ \chi }{ \uplambda } \int \limits _{ 0 }^{ +\infty } \frac{ e^{ -2 l }\, dl }{ ( \gamma + l ) [ ( \gamma + l )^{ 2 } - e^{ -2 l } ] } \nonumber \\&+ \frac{ b_{ 1 } \chi }{ \gamma } \int \limits _{ 0 }^{ +\infty } \frac{ l^{ 2 } [ 3 ( \gamma + l )^{ 2 } e^{ -2 l } - e^{ -4 l } ] }{ ( \gamma + l )^{ 2 } [ ( \gamma + l )^{ 2 } - e^{ -2 l } ]^{ 2 } } \mathrm{d}l - \frac{ 2 }{ \gamma } \int \limits _{ 0 }^{ +\infty } \frac{ l e^{ -2 l } \, dl }{ ( \gamma + l ) [ ( \gamma + l )^{ 2 } - e^{ -2 l } ] } \nonumber \\&+ \frac{ 1 }{ \gamma } \int \limits _{ 0 }^{ +\infty } \frac{ ( 1 - l ) e^{ -2 l } }{ ( \gamma + l )^{ 2 } - e^{ -2 l } } \mathrm{d}l + R( a, \uplambda ) \Bigg ], \end{aligned}$$
(5)
$$\begin{aligned}&|{ R( a, \uplambda ) }| \le \text {const}( a )\, \uplambda , \end{aligned}$$
(6)

where

$$\begin{aligned} \gamma = \frac{ \alpha a }{ 2 \pi ^{ 2 } } > 1. \end{aligned}$$
(7)

The term \(\mathcal {E}_{ \text {self} }( \uplambda )\) does not depend on the separation \(\vec {a}\) and represents the self-energy of two independent interaction centers. The self-energy diverges in the limit, as should be expected, since in that situation it describes the self-energy of two \(\delta ( \vec { x })\). Because it does not contribute to the Casimir force, its divergence poses no problem for the force. The constants \(\alpha \), \(\chi \) and \(b_1\) characterize the initial model \(V_{ \vec { a } }\), which undergoes the scaling. However, their status is not equal. In the process of limiting the model tends to one of the \(\delta \)-models, which are parameterized only by \(\alpha \), the other two parameters die out in the resolvent limit. On the other hand, they survive in the expression for energy, which therefore is very strongly model-dependent.

As mentioned above, for the consistency of the model one has to demand that \(h_{ \vec { a } }^{ \;\; 2 }\) is positive. As it turns out, this is equivalent to the condition \(\gamma \ge 1\), which shows the presence in the model of a critical minimal value of the separation of interaction centers (the reader interested in details of this problem should go to chapter II.1.1 of the monograph [1]). For technical reasons, we adopt a marginally stronger form (7) of this condition.

The \(\alpha \)-dependent limit models form the well-known family extensively investigated in the monograph on solvable models in quantum mechanics by Albeverio, et al. Solvable Models in Quantum Mechanics [1] (see chapter II.1).Footnote 1 However, we found it appropriate to briefly discuss our limiting procedure leading to the strict \(\delta \)-model, as it differs from that used by Albeverio et al. (chapter II.1.2 of the aforementioned work). Also, we stress the fact that our object of interest is not the limit model itself, but the approximating family of models.

2 The Model

In this article we investigate a class of models describing massless scalar field \(\phi ( t, \vec { x })\) interacting with two nonsingular spherically symmetric, nonoverlapping macroscopic bodies. This system is modeled by the following choice of the operator \(h_{ \vec { a }}^{ \;\; 2 }\):

$$\begin{aligned}&h_{ \vec { a }}^{ \;\; 2 } = -\Delta + V_{ \vec { a } }, \end{aligned}$$
(8a)
$$\begin{aligned}&\begin{aligned} V_{ \vec { a }}&= \sigma ( g )\left( | U_{ \vec { b }}\,g \rangle \langle U_{ \vec { b }}\,g | + | U_{ -\vec { b }}\,g \rangle \langle U_{ -\vec { b }}\,g | \right) \\&= \sigma ( g )( | U_{ \vec { b }}\,g \rangle \; | U_{ -\vec { b }}\,g \rangle ) \begin{pmatrix} \langle U_{ \vec { b }}\,g | \\ \langle U_{ -\vec { b }}\,g | \end{pmatrix} , \quad \vec { b }= \frac{ 1 }{ 2 } \vec { a }. \end{aligned} \end{aligned}$$
(8b)

Here \(g( \vec { x }) = g( |{\, \vec { x } \,}| )\) is a \(\mathcal {C}^{ \infty }( \mathbb {R}^{ 3 } )\), spherically-symmetric function with compact support and \(\sigma ( g )\) is a real g-dependent parameter to be defined below. Moreover, we assume for technical reasons that \(g( \vec { x })\) is real and positive. We also use notation \(U_{ \vec { b }}\,\) for the operator of translation by \(\vec { b }\).

In position representation the quasi-potential \(V_{ \vec { a }}\) is an integral operator with the kernel

$$\begin{aligned} V_{ \vec { a }}( \vec { x }, \vec { y }) = \sigma ( g )\left[ g( \vec { x }- \vec { b }\, )\, \overline{ g( \vec { y }- \vec { b }\, ) } + g( \vec { x }+ \vec { b }\, )\, \overline{ g( \vec { y }+ \vec { b }\, ) } \, \right] . \end{aligned}$$
(9)

Because we assume that the two bodies do not overlap, we set \(|{\, \vec { a } \,}| > L\), where \(\mathrm {supp}\,g \subset B( 0, L / 2 )\). The Fourier transform of \(g( \vec { x })\) is a real spherically-symmetric, analytic function

$$\begin{aligned} \begin{aligned} \widehat{g}( p ) = \widehat{g}( \vec { p }\, ) = \frac{ 1 }{ \sqrt{ 2 \pi } } \, \frac{ i }{ p } \int \limits _{ -\infty }^{ +\infty } \mathrm{d} x \, e^{ -i px } x g( x ) = \frac{ 1 }{ \sqrt{ 2 \pi } } \int \limits _{ -\infty }^{ +\infty } \mathrm{d} x \, x^{ 2 } \frac{ \sin px }{ px } g( x ), \end{aligned} \end{aligned}$$
(10)

where \(p = |{\, \vec { p } \,}|\). In this equation we use the extension of g(x) considered as the function of radius to the even function on \(\mathbb {R}\). This function is smooth, which is evident by the fact, that such extension is just the restriction of \(g(\vec {x})\) to one particular axis. In the same way we extend the function \(\widehat{g}( p )\). After such procedure

$$\begin{aligned} M_{ p }= \widehat{g}( p )^{ 2 } \end{aligned}$$
(11)

is an even, smooth function with the analytic extension to the complex plane.

The function \(\sigma ( g )\) is chosen as

$$\begin{aligned} \sigma ( g )^{ -1 } = -\alpha ( M_0 + \chi ( g ) ), \quad \alpha > 0, \end{aligned}$$
(12)

where \(\alpha \) is a parameter of the model and

$$\begin{aligned} \chi ( g ) = \frac{ 2 \pi }{ \alpha } \int \limits _{ -\infty }^{ +\infty } \mathrm{d} p \, M_{ p } = \sqrt{ 2 \pi }^{ 3 } \frac{ \check{M}_{ 0 }}{ \alpha }. \end{aligned}$$
(13)

Function \(\check{M}_{ x }\) is the one-dimensional inverse Fourier transform of \(M_{ p }\). For the sake of brevity we often suppress the argument of \(\chi \). This choice of \(\sigma (g)\) is motivated by the following perspective. We shall later consider a 1-parameter class of rescaled functions \(g_{ \uplambda }( \vec { x })\) (see Sect. 4) and a scaling limit \(\uplambda \searrow 0\); \(\sigma ( g )\) in the form above will ensure the existence of a well defined, physically interesting and non-trivial limit. The choice (12) makes the operator \(V_{ \vec { a }}\) invariant under the operation \(g( \vec { x }) \rightarrow C g( \vec { x })\), \(C> 0\). It is convenient to eliminate this freedom by demanding the following normalization:

$$\begin{aligned} \widehat{g}( 0 ) = \frac{ 1 }{\sqrt{ 2\pi }^{ 3 }} \int \limits _{ \mathbb {R}^{ 3 }} g( \vec { x })\, d^{ 3 }x = 1, \end{aligned}$$
(14)

so that

$$\begin{aligned} M_{ 0 } = \, \widehat{g}( 0 ) = 1, \end{aligned}$$
(15)

and

$$\begin{aligned} \sigma ( g )^{ -1 } = -\alpha ( 1 + \chi ( g ) ). \end{aligned}$$
(16)

The appearance of the parameter \(\alpha \) in the problem with two identical \(\delta ( \vec { x })\) is well established and known in literature. Such models, and their relation to less rigorous approach to \(\delta \)-potentials, are extensively discussed in “Introduction” and chapter II.1 of the monograph [1].

3 Resolvent and Positivity

In this section we investigate the resolvent operator of our model. We use the following notation

$$\begin{aligned} G_{ 0 }( w^{ 2 } ) = \left( w^{ 2 } - h^{ 2 } \right) ^{ -1 }, \qquad G( w^{ 2 } ) = \left( w^{ 2 } - h_{ \vec { a }}^{ \;\; 2 } \right) ^{ -1 }. \end{aligned}$$
(17)

We shall show that in the range of the separation of the interaction centers given by

$$\begin{aligned} a> \frac{ 2\pi ^{ 2 } }{ \alpha } \iff \gamma := \frac{ \alpha a }{ 2 \pi ^{ 2 } } > 1, \end{aligned}$$
(18)

the resolvent \(G( w^{ 2 } )\) is an analytic function for all \(w^{ 2 }\) outside the real, non-negative semi-axis, thus the operator \(h_{ \vec { a }}^{ \;\; 2 }\) is positive. Given \(w^{ 2 }\) we make w unique by demanding that \(\mathfrak {Im}\, w \ge 0\), and then the region of analyticity is covered by \(\mathfrak {Im}\, w>0\).

Because the quasi-potential \(V_{ \vec { a }}\) is a trace class operator, the Kato–Rosemblum theorem guaranties that standard scattering theory applies to this case (see [15], Thm. XI.8 and the discussion of the scattering theory at the beginning of chapter XI.3). For the purpose of this section we will need two equations from this theory

$$\begin{aligned} T( w^{ 2 } )&= V_{ \vec { a }}+ V_{ \vec { a }}\, G_{ 0 }( w^{ 2 } ) T( w^{ 2 } ), \end{aligned}$$
(19a)
$$\begin{aligned} G( w^{ 2 } )&= G_{ 0 }( w^{ 2 } ) + G_{ 0 }( w^{ 2 } ) T( w^{ 2 } )G_{ 0 }( w^{ 2 } ), \end{aligned}$$
(19b)

which follow from the definition of \(T( w^{ 2 } )\) and the Lippmann–Schwinger equations (see [19], chapters 8-a and 8-b).

To solve Eq. (19a) for \(T( w^{ 2 } )\) we again follow [11] and postulate that

$$\begin{aligned} T( w^{ 2 } )= ( | U_{ \vec { b }}\,g \rangle \; | U_{ -\vec { b }}\,g \rangle ) \mathcal {T}( w )\begin{pmatrix} \langle U_{ \vec { b }}\,g | \\ \langle U_{ -\vec { b }}\,g | \end{pmatrix} , \end{aligned}$$
(20)

where \(\mathcal {T}( w )\) is a \(2 \times 2\) complex matrix. Using this formula in Eq. (19a) we find

$$\begin{aligned} \begin{aligned} \mathcal {T}( w )&= \begin{pmatrix} -t( w ) &{} u( w ) \\ u( w ) &{} -t( w ) \end{pmatrix}^{ -1 }, \end{aligned} \end{aligned}$$
(21)

where

$$\begin{aligned} t( w )&= -\sigma ( g )^{ -1 }+ ( g,\, G_{ 0 }( w^{ 2 } ) g ) = \alpha + 2 \pi \int \limits _{ -\infty }^{ +\infty } \mathrm{d} p \, \frac{ w^{ 2 } }{ w^{ 2 } - p^{ 2 } } M_{ p }, \end{aligned}$$
(22a)
$$\begin{aligned} u( w )&= -( U_{ \vec { a }} g,\, G( w^{ 2 } ) g ) = 2 \pi ^{ 2 } M_{ w } \frac{ e^{ i w a } }{ a }; \end{aligned}$$
(22b)

we used here the explicit form of \(\sigma ( g )\) from (16) and for brevity we suppress the dependence on g. Also, to compute (22b) we used inequality (111) from “Appendix A.1”. It is easy to see that the operator norm of the matrix \(\mathcal {T}( w )\) is

$$\begin{aligned} \Vert \mathcal {T}( w ) \Vert = \max \left\{ \frac{ 1 }{ |{ t( w ) + u( w ) }| },\, \frac{ 1 }{ |{ t( w ) - u( w ) }| } \right\} . \end{aligned}$$
(23)

In “Appendix A.3” we show that

$$\begin{aligned} \mathfrak {Re}\, \int \limits _{ -\infty }^{ +\infty } \mathrm{d} u \, \frac{ w^{ 2 } }{ w^{ 2 } - u^{ 2 } } M_{ u } \ge 0, \quad \mathfrak {Im}\, w > 0, \end{aligned}$$
(24)

which implies that \(|{ \mathfrak {Re}\, t( w ) }| \ge \alpha \). Using estimate (111) from “Appendix A.1”, we also have the bound \(|{ u( w ) }| \le 2\pi ^{ 2 } / a\), which implies

$$\begin{aligned} \Vert \mathcal {T}( w ) \Vert \le \frac{ 1 }{ \alpha - 2\pi ^{ 2 } / a }, \quad \mathfrak {Im}\, w > 0, \end{aligned}$$
(25)

for all functions \(g( \vec { x })\) in the considered class. With the use of notation

$$\begin{aligned} \mathcal {F}_{ \vec { p } } = \widehat{g}( p ) ( e^{ -i \vec { b } \cdot \vec { p } } \quad e^{ i \vec { b } \cdot \vec { p } }\, ), \end{aligned}$$
(26)

we can express the integral kernel in the momentum space as

$$\begin{aligned} \begin{aligned} \langle \vec { p }\, | G( w^{ 2 } ) - G_{ 0 }( w^{ 2 } ) | \vec { q }\, \rangle = \frac{ \mathcal {F}_{ \vec { p } } }{ w^{ 2 } - \vec { p\, }^{ 2 } } \mathcal {T}( w )\frac{ \mathcal {F}_{ \vec { q } }^{ \; \dag } }{ w^{ 2 } - \vec { q\, }^{ 2 } }. \end{aligned} \end{aligned}$$
(27)

Positivity of the operator \(h_{ \vec { a }}^{ \;\; 2 }\) is equivalent to the existence of its resolvent on the real negative semi axis, i.e. for \(w = i \mu \), \(\mu > 0\). Because \(\Vert \mathcal {T}( i\mu ) \Vert \) is bounded by a constant, and

$$\begin{aligned} \int \frac{ \mathcal {F}_{ \vec { p }}^{ \; \dagger } \mathcal {F}_{ \vec { p }} }{ ( \mu ^{ 2 } + \vec { p }^{ \, 2 } )^{ 2 } }\, d^{ 3 }p \le \text {const} \frac{ 1 }{ \mu } \int \frac{ 1 }{ ( 1 + \vec { p }^{ \, 2 } )^{ 2 } }\, d^{ 3 }p, \end{aligned}$$
(28)

the difference \(G( -\mu ^{ 2 } ) - G_{ 0 }( -\mu ^{ 2 } )\) is a Hilbert–Schmidt operator and

$$\begin{aligned} \Vert G( -\mu ^{ 2 } ) - G_{ 0 }( -\mu ^{ 2 } ) \Vert \le \Vert G( -\mu ^{ 2 } ) - G_{ 0 }( -\mu ^{ 2 } ) \Vert _{ \text {HS} } \le \text {const} \frac{ 1 }{ \mu }. \end{aligned}$$
(29)

As \(\Vert G_{ 0 }( -\mu ^{ 2 } ) \Vert \le \text {const} / \mu ^{ 2 }\), we have

$$\begin{aligned} \Vert G( -\mu ^{ 2 } ) \Vert \le \text {const} \left( \frac{ 1 }{ \mu } + \frac{ 1 }{ \mu ^{ 2 } } \right) . \end{aligned}$$
(30)

This bound shows that \(h_{ \vec { a }}^{ \;\; 2 }\) is a non-negative operator. Also, for further use we note that it does not depend on the choice of function \(g( \vec { x })\).

4 Rescaled Model and Resolvent Limit

We now introduce a family of rescaled models with the quasi-potential \(V_{ \vec { a },\, \uplambda }\), with \(\uplambda \in ( 0, 1 ]\), by replacing function \(g( \vec { x })\) by \(g_{ \uplambda }( \vec { x })\), where

$$\begin{aligned} \begin{aligned} g_{ \uplambda }( \vec { x }) = \uplambda ^{ -3 } g( \tfrac{ \vec { x }}{ \uplambda } ) \iff \widehat{g}_{ \uplambda }( p ) = \widehat{g}( \uplambda p ), \end{aligned} \end{aligned}$$
(31)

and in consequence \(M_{ \uplambda ,\, p } = M_{ \uplambda p }\). The rescaled functions satisfy all our assumptions on functions \(g( \vec { x })\), in particular the scaling factor \(\uplambda ^{ -3 }\) ensures the validity of normalization condition (15).Footnote 2

We want to investigate the behavior of the system when \(\uplambda \searrow 0\). This procedure has a clear interpretation and is justified from physical point of view. Namely, we want to know what happens when the local strength of interaction of two point-like objects grows to infinity, and at the same time their size is shrinking to zero.

The basic quantities of the rescaled model are now

$$\begin{aligned}&\begin{aligned} \mathcal {T}_{ \uplambda }( w )&= \begin{pmatrix} -t_{ \uplambda }( w ) &{} u_{ \uplambda }( w ) \\ u_{ \uplambda }( w ) &{} -t_{ \uplambda }( w ) \\ \end{pmatrix}^{ -1 }, \end{aligned} \end{aligned}$$
(32)
$$\begin{aligned}&\mathcal {F}_{ \uplambda ,\, { \vec { p } } } = \widehat{g}( \uplambda p ) ( e^{ -i \vec { b } \cdot \vec { p } } \quad e^{ i \vec { b } \cdot \vec { p } }\, ), \end{aligned}$$
(33)
$$\begin{aligned}&\langle \vec { p }\, | G_{ \uplambda }( w^{ 2 } ) - G_{ 0 }( w^{ 2 } ) | \vec { q }\, \rangle = \frac{ \mathcal {F}_{ \uplambda ,\, { \vec { p } } } }{ w^{ 2 } - \vec { p }\,^{ 2 } } \mathcal {T}_{ \uplambda }( w )\frac{ \mathcal {F}_{ \uplambda ,\, { \vec { q } } }^{ \dagger } }{ w^{ 2 } - \vec { q }\,^{ 2 } }, \end{aligned}$$
(34)
$$\begin{aligned}&\begin{aligned} t_{ \uplambda }( w ) = -\sigma ( g_{ \uplambda })^{ -1 }+ ( g_{ \uplambda },\, G_{ 0 }( w^{ 2 } ) g_{ \uplambda }) = \alpha + 2 \pi \int \limits _{ -\infty }^{ +\infty } \mathrm{d}p \, \frac{ w^{ 2 } }{ w^{ 2 } - p^{ 2 } } M_{ \uplambda p }, \end{aligned}\end{aligned}$$
(35a)
$$\begin{aligned}&u_{ \uplambda }( w ) = -( U_{ \vec { a } }\, g_{ \uplambda }, G_{ 0 }( w^{ 2 } ) g_{ \uplambda }) = 2 \pi ^{ 2 } M_{ \uplambda w } \frac{ e^{ i w a } }{ a } . \end{aligned}$$
(35b)

The reader can check that thanks to our choice of \(\sigma ( g_{ \uplambda } )\), in Eq. (35a) two terms most singular in \(\uplambda \) cancel out. Also, it is important to note, that bound (25) for \(\Vert \mathcal {T}_{ \uplambda }( w ) \Vert \) is independent of the function \(g( \vec { x })\), which implies its independence from \(\uplambda \).

We now want to find the limit of operators \(h_{ \uplambda ,\, \vec { a }}^{ \; 2 }\) in the appropriate sense. We shall show below that there exists strong limit

$$\begin{aligned} \lim \limits _{ \uplambda \searrow 0 } G_{ \uplambda }( -\mu ^{ 2 } ) = {\widetilde{G}}( -\mu ^{ 2 } ), \qquad \mu > 0. \end{aligned}$$
(36)

We note that then \({\widetilde{G}}( -\mu ^{ 2 } )\) satisfies the bound (30). Therefore, arguing as in the proof of Trotter–Kato theorem (see Thm. VIII.22 in [16]), we can conclude, that \({\widetilde{G}}\) extends to an analytic operator valued function with the domain \(\mathbb {C} \setminus [ 0, +\infty )\). Moreover, from the explicit formula (37) below it is easy to see that \({\widetilde{G}}( -1 )\) is a bounded, strictly negative operator, so its range is dense in the Hilbert space. Therefore, \({\widetilde{G}}( w^{ 2 } )\) is the resolvent of a self-adjoint operator \({\widetilde{h}}_{ \vec { a }}^{ \;\; 2 }\), the so-called resolvent limit of \(h_{ \uplambda ,\, \vec { a }}^{ \; 2 }\).

To show (36) it is sufficient to prove that the integral kernel of the operator \(\langle \vec { p }\, | G_{ \uplambda }( -\mu ^{ 2 } ) - G_{ 0 }( -\mu ^{ 2 } ) | \vec { q }\, \rangle \) has the \(L^{ 2 }\)-limit as a function of \(\vec { p }\), \(\vec { q }\). This is an easy exercise which we omit. The limit function is easily extended to an analytic function:

$$\begin{aligned}&\langle \vec { p }\, | {\widetilde{G}}( w^{ 2 } ) - G_{ 0 }( w^{ 2 } ) | \vec { q }\, \rangle = - \frac{ 1 }{ ( \alpha - i 2 \pi ^{ 2 } w )^{ 2 } - ( 2 \pi ^{ 2 } e^{ i wa } / a )^{ 2 } } \nonumber \\&\quad \times \frac{ 1 }{ w^{ 2 } - \vec { p\, }^{ 2 } } ( e^{ -i \vec { b }\cdot \vec { p }} \quad e^{ i \vec { b }\cdot \vec { p }}\, ) \begin{pmatrix} \alpha - i 2 \pi ^{ 2 } w &{} 2 \pi ^{ 2 } \tfrac{ e^{ i wa }}{ a } \\ 2 \pi ^{ 2 } \tfrac{ e^{ i wa } }{ a } &{} \alpha - i 2 \pi ^{ 2 } w \\ \end{pmatrix} \begin{pmatrix} e^{ i \vec { b }\cdot \vec { q }} \\ e^{ -i \vec { b }\cdot \vec { q }} \end{pmatrix} \frac{ 1 }{ w^{ 2 } - \vec { q\, }^{ 2 } } .\nonumber \\ \end{aligned}$$
(37)

The resolvent \({\widetilde{G}}( w^{ 2 } )\) is known to describe the two \(\delta \)-model (see [1], section II.1.1.).

5 Spectral Analysis

The function \(g_{ \uplambda }( \vec { x })\) is in the class of all functions \(g( \vec { x })\) admitted in our model, so the following analysis applies both to the original, as well as to the rescaled versions of the model.

We already know that \(h_{ \vec { a }}^{ \;\; 2 }\) is a non-negative operator. Also, as \(h_{ \vec { a }}^{ \;\; 2 } - h^{ 2 }\) is a finite rank, thus compact operator, the essential spectrum of \(h_{ \vec { a }}^{ \;\; 2 }\) is \([ 0, \infty )\), as a result of Weyl Theorem (see [14], Thm. XIII.14 and Example 3 in section XIII.4). Using the form of \(G( -\mu ^{ 2 } )\) one easily shows that the resolvent has no pole singularity in \(\mu ^{ 2 } = 0\), so there is no zero eigenvalue. We shall show below, that for each \(0< a < b\) and each \(f \in \mathcal { S }( \mathbb {R}^{ 3 })\), the function

$$\begin{aligned} ( a, b ) \times i ( 0, 1 ] \ni w^{ 2 } \mapsto ( f,\, G( w^{ 2 } ) f ), \end{aligned}$$
(38)

is bounded. Then by standard theorems (see theorems XIII.19 and XIII.20 in [14]) the operator has no point or singular spectrum in (a, b).

To fill the missing step, we proceed as follow. Because

$$\begin{aligned} ( f,\, G( w^{ 2 } ) f ) = ( f,\, G_{ 0 }( w^{ 2 } ) f ) + ( f,\, [ G( w^{ 2 } ) - G_{ 0 }( w^{ 2 } ) ] f ) \end{aligned}$$
(39)

and \(\Vert \mathcal {T}( w ) \Vert \le \text {const}\) for all \(\mathfrak {Im}\, w > 0\) if \(\gamma = \alpha a / 2\pi ^{ 2 } > 1\), the problem is reduced to showing that for every function F in the Schwartz space \(\mathcal {S}( \mathbb {R}^{ 3 } )\) the term of the following form is bounded

$$\begin{aligned} \begin{aligned} \int \limits _{ \mathbb {R}^{ 3 }} d^{ 3 }p\, \frac{ F( \vec { p }\, ) }{ w^{ 2 } - \vec { p }\,^{ 2 } } = \int \limits _{ 0 }^{ +\infty } \mathrm{d} p \, \frac{ f( p ) }{ w^{ 2 } - p^{ 2 } }, \quad f( p ) = p^{ 2 } \int \limits _{ S^{ 2 } } \mathrm{d}\Omega ( \widehat{n})\, f( p\widehat{n}). \end{aligned} \end{aligned}$$
(40)

The function f(p) can be extended to a Schwartz function on the whole real line. After putting \(w^{ 2 } = k^{ 2 } + i \epsilon \), the integral takes the form

$$\begin{aligned} \begin{aligned} \int \limits _{ 0 }^{ +\infty } \mathrm{d} p \, f( p ) \frac{ 1 }{ k^{ 2 } + i \epsilon - p^{ 2 } }&= -\frac{ i \pi }{ 2 } \frac{ f( k ) }{ \sqrt{ k^{ 2 } + i \epsilon } } + \int \limits _{ 0 }^{ +\infty } \mathrm{d} p \, \frac{ f( p ) - f( k ) }{ k^{ 2 } - p^{ 2 } } \frac{ k^{ 2 } - p^{ 2 } }{ k^{ 2 } - p^{ 2 } + i \epsilon }. \end{aligned}\nonumber \\ \end{aligned}$$
(41)

Consider now \(k^{ 2 } \in ( a,\, b )\), \(0< a < b\). It is easy to see that the first term is bounded on this set. As for the second term, we first observe that for \(k>\sqrt{a}\), \(p\ge 0\), the function \(( f( p ) - f( k ) ) / ( k^{ 2 } - p^{ 2 } ) \) is smooth as a function of two variables, due to the well known fact that the difference quotient of a smooth function is also smooth. Because for \(p > b\) this function is bounded by \(\text {const} / p^{ 2 }\), we find that above expression is bounded on the set in question.

This means that the condition (18) is sufficient to ensure that \(h_{ \vec { a }}^{ \;\; 2 }\) has no singular or point spectrum and its absolutely continuous spectrum is equal to \([ 0, +\infty )\).

6 Casimir Energy and Number of Particles

We now turn to our main task in this article – calculation of the Casimir energy within the framework explained in Sect. 1. The key formulas for the energy and the number of particles are given in (4a) and (4b), respectively. Using the spectral properties of h in momentum space, we write the rhs of these formulas in the form

$$\begin{aligned} \mathrm {Tr}\, \mathcal {P}_{ \tau }( a ) = \int \limits _{ \mathbb {R}^{ 3 }} d^{ 3 }p\, p^{ -\tau } \langle \vec { p }\, | ( h_{ \vec { a }} - p ) h_{ \vec { a }}^{ \;\; -1 } ( h_{ \vec { a }} - p ) |\, \vec { p }\, \rangle , \end{aligned}$$
(42)

where \(\tau = 0\) for energy and \(\tau = 1\) for the number of particles. Since \(h_{ \vec { a }}\) has only a continuous spectrum equal to \([ 0, +\infty )\), we can use its spectral decomposition to transform the above formula into

$$\begin{aligned} \mathrm {Tr}\, \mathcal {P}_{ \tau }( a ) = \int \limits _{ \mathbb {R}^{ 6 } } \mathrm{d}^{ 3 } p \, \mathrm{d}^{ 3 } k \; \frac{ p^{ -\tau } }{ k } ( k - p )^{ 2 } | \langle \vec { p }\, | \, \vec { k }+ \rangle |^{ 2 }, \end{aligned}$$
(43)

where \(| \vec { k }+ \rangle \) denotes the standard eigenfunction of \(h_{ \vec { a }}\) (see [19], chapter 10-a). We now need another identity from scattering theory (see [19], chapter 8-c)

$$\begin{aligned} ( k - p ) \langle \vec { p }\, | \vec { k }+ \rangle = \frac{ \langle \vec { p }\, | \, T( k^{ 2 } + i 0 ) | \, \vec { k }\, \rangle }{ p + k }, \end{aligned}$$
(44)

where the Dirac delta contribution has vanished due to multiplication by \(k - p\). Using this, we obtain

$$\begin{aligned} \begin{aligned} \mathrm {Tr}\, \mathcal {P}_{ \tau }( a )&= \int \limits _{ \mathbb {R}^{ 6 } } \mathrm{d}^{ 3 } k \, \mathrm{d}^{ 3 } p \; \frac{ p^{ -\tau } }{ k ( p + k )^{ 2 } } M_{ p } M_{ k } \\&\quad \times \Bigl | \left( e^{ - i \vec { b }\cdot \vec { p }} \,\,\, e^{ i \vec { b }\cdot \vec { p }}\, \right) \mathcal {T}( k + i 0 ) \left( e^{ - i \vec { b }\cdot \vec { k }} \,\,\, e^{ i \vec { b }\cdot \vec { k }}\, \right) ^{ \dagger } \Bigr |^{ 2 }. \end{aligned} \end{aligned}$$
(45)

This can be simplified with the use of a relation which holds for any row matrices A, C and a square matrix B of the same dimension

$$\begin{aligned} \begin{aligned} |{ A B C^{ \dag } }|^{ 2 } = \mathrm {Tr}[ A^{ \dag } A B C^{ \dag } C B^{ \dag } ], \end{aligned} \end{aligned}$$
(46)

from which it follows

$$\begin{aligned}&\begin{aligned}&\bigl | ( e^{ -i \vec { b }\cdot \vec { p }} \,\,\, e^{ i \vec { b }\cdot \vec { p }}\, )\, \mathcal {T}( k + i 0 )\, ( e^{ -i \vec { b }\cdot \vec { k }} \,\,\, e^{ i \vec { b }\cdot \vec { k }}\, )^{ \dag } \bigr |^{ 2 } \\&\quad = \mathrm {Tr}\left[ \mathcal {M}( \vec { p }\, ) \mathcal {T}( k + i0 ) \mathcal {M}( \vec { k }\, ) \mathcal {T}^{ \dag }( k + i 0 ) \right] , \end{aligned} \end{aligned}$$
(47a)
$$\begin{aligned}&\mathcal {M}( \vec { k }\, ) = \left( e^{ -i \vec { b }\cdot \vec { k }} \,\,\, e^{ i \vec { b }\cdot \vec { k }}\, \right) ^{ \dagger } \left( e^{ -i \vec { b }\cdot \vec { k }} \,\,\, e^{ i \vec { b }\cdot \vec { k }}\, \right) = \begin{pmatrix} 1 &{} e^{ i \vec { k }\cdot \vec { a }} \\ e^{ -i \vec { k }\cdot \vec { a }} &{} 1 \\ \end{pmatrix} . \end{aligned}$$
(47b)

We now define functions t(k) and u(k) for a real argument k as the boundary values of the previously introduced functions:

$$\begin{aligned} t( k )&= \alpha + 2 \pi k^{ 2 } I_{ k } - i 2 \pi ^{ 2 } k M_{ k }, \end{aligned}$$
(48a)
$$\begin{aligned} u( k )&= 2 \pi ^{ 2 } M_{ k } \frac{ e^{ i ka } }{ a }, \end{aligned}$$
(48b)

where following [11] we denote

$$\begin{aligned} I_{ k } = \int \limits _{ -\infty }^{ +\infty } \mathrm{d} p \,{ \frac{ M_{ p } - M_{ k } }{ k^{ 2 } - p^{ 2 } } }. \end{aligned}$$
(49)

Computation of (48a) can be done using the method from Sect. 5. Then

$$\begin{aligned} \begin{aligned}&\mathcal {T}( k + i 0 ) = \begin{pmatrix} -t( k ) &{} u( k ) \\ u( k ) &{} -t( k ) \\ \end{pmatrix}^{ -1 } = -\frac{ 1 }{ t( k )^{ 2 } - u( k )^{ 2 } } \begin{pmatrix} t( k ) &{} u( k ) \\ u( k ) &{} t( k ) \\ \end{pmatrix} . \end{aligned} \end{aligned}$$
(50)

The matrix \(\mathcal {T}( k + i 0 )\) is well defined and bounded on the positive semi-axis, because the bound (25) remains true in the limit \(\mathfrak {Im}\, w \searrow 0\). In the following computations we will need the imaginary part of the matrix \(\mathcal {T}( k + i 0 )^{ -1 }\)

$$\begin{aligned} \mathfrak {Im}\, \mathcal {T}( k + i0 )^{ -1 } = 2 \pi ^{ 2 } k M_{ k } \begin{pmatrix} 1 &{}\quad \frac{ \sin ( ka ) }{ ka } \\ \frac{ \sin ( ka ) }{ ka } &{}\quad 1 \\ \end{pmatrix}. \end{aligned}$$
(51)

Taking insight from the work in [11], we introduce

$$\begin{aligned} \mathcal {L}( k ) = \int \limits _{ S^{ 2 }( 0,\, 1 ) } \mathrm{d}\Omega ( {\hat{k}} )\, \mathcal {M}( \vec { k }\, ) = 4 \pi \begin{pmatrix} 1 &{}\quad \frac{ \sin ( ka ) }{ ka } \\ \frac{ \sin ( ka ) }{ ka } &{}\quad 1 \\ \end{pmatrix}, \quad {\hat{k}} = \vec { k }/ | \vec { k }|, \end{aligned}$$
(52)

and point out that

$$\begin{aligned} \mathfrak {Im}\, \mathcal {T}( k + i0 )^{ -1 } = \frac{ \pi }{ 2 } k M_{ k } \mathcal {L}( k ). \end{aligned}$$
(53)

From this we have

$$\begin{aligned} \mathcal {L}( k ) = \frac{ i }{ \pi k M_{ k } } ( \mathcal {T}^{ \dag }( k + i0 )^{ -1 } - \mathcal {T}( k + i0 )^{ -1 } ). \end{aligned}$$
(54)

After integrating over the angles, the trace expressions take the form

$$\begin{aligned} \begin{aligned} \mathrm {Tr}\, \mathcal {P}_{ \tau }( a )&= \frac{ 2 }{ \pi } \int \limits _{ \mathbb {R}_{ + }^{ 2 } } \mathrm{d}^{ } k \, \mathrm{d}^{ } p \; \frac{ p^{ 2 - \tau } }{ ( p + k )^{ 2 } } M_{ p }\, \mathfrak {Re}\left\{ \mathrm {Tr}\left[ i\mathcal {L}( p ) \mathcal {T}\left( k + i 0 \right) \right] \right\} , \end{aligned} \end{aligned}$$
(55)

where the trace is bounded by a constant (see (25) and (51)) so these integrals are convergent. If we denote

$$\begin{aligned} h( x ) = \frac{ \sin ( x ) }{ x }, \end{aligned}$$
(56)

we can write this as

$$\begin{aligned} \mathrm {Tr}\mathcal {P}_{ \tau }( a ) = \mathfrak {Im}\int \limits _{ \mathbb {R}_{ + }^{ 2 } } \mathrm{d}^{ } k \, \mathrm{d}^{ } p \; f_{ \tau }( k, p ), \end{aligned}$$
(57)

where \(f_{ \tau }( k, p )\) is the boundary value of the function

$$\begin{aligned} f_{ \tau }( z, p ) = 16 \frac{ p^{ 2 - \tau } }{ ( p + z )^{ 2 } } M_{ p } \frac{ t( z ) + h( p a ) u( z ) }{ t( z )^{ 2 } - u( z )^{ 2 } }, \end{aligned}$$
(58)

which is analytical in the upper half-plane. The evaluation of these formulas can be greatly simplified when the integrals over k are rotated in the complex plane into the imaginary half-axis, using the contour of integration depicted in the Fig. 1. It is easily observed that for \(\mathfrak {Im}\, z \ge 0\) and fixed p

$$\begin{aligned} |{ f_{ \tau }( z, p ) }| \le \text {const} \, \frac{ 1 }{ \; |{ z + p }|^{ 2 } }, \end{aligned}$$
(59)

so the integral over \(\gamma _{ R }\) vanishes in the limit \(R \rightarrow +\infty \). Therefore, using properties of analytical functions we can change the integral over the path \(\gamma _{ 1 }\) to that over \(-\gamma _{ 2 }\). In particular the Casimir energy integral takes the form

$$\begin{aligned} \begin{aligned} \mathcal {E}( a )&= \mathfrak {Re}\int \limits _{ \mathbb {R}_{ + }^{ 2 } } \mathrm{d}^{ } p \, \mathrm{d}^{ } k \; f_{ 0 }( ik, p ). \end{aligned} \end{aligned}$$
(60)
Fig. 1
figure 1

Contour of integration of function f(z)

Full size image

We note that \(|{ h( \mathrm{pa} ) }|\) and \(|{ u( z ) }|\) are uniformly bounded in a and tend to zero for \(a \rightarrow +\infty \) [for u(z) use bound (111)]. Therefore

$$\begin{aligned} \mathcal {E}_{ 0 } = \lim _{ a \rightarrow +\infty } \mathcal {E}( a ) = 16\, \mathfrak {Re}\int \limits _{ \mathbb {R}_{ + }^{ 2 } } \mathrm{d}p\, \mathrm{d}k\, \frac{ p^{ 2 } }{ ( p + i k )^{ 2 } } \, \frac{ M_{ p } }{ t( ik ) }. \end{aligned}$$
(61)

The physical sense of this formula is as follows: it describes the energy of two centers being “infinitely” apart, so it is equal to twice the self-energy of one center in empty space. Subtracting it from the total energy we obtain the interaction energy

$$\begin{aligned} \mathcal {E}_{ \text {int} }( a ) = 16\, \mathfrak {Re}\int \limits _{ \mathbb {R}_{ + }^{ 2 } } \mathrm{d}^{ } p \, \mathrm{d}^{ } k \; \frac{ p^{ 2 } }{ ( p + i k )^{ 2 } } \, M_{ p } \, \frac{ u( ik )^{ 2 } + h( pa ) t( i k ) u( i k ) }{ t( ik ) [ t( ik )^{ 2 } - u( i k )^{ 2 } ] }. \end{aligned}$$
(62)

7 Casimir Energy of The Rescaled Model

We recall that our model with a scaled \(g_{ \uplambda }( \vec { x })\) replacing \(g( \vec { x })\) approximates for \(\uplambda \searrow 0\) the system of two point-like interaction centers. Their separation a has to satisfy two restrictions:

  1. (i)

    \(a > \frac{ 2\pi ^{ 2 } }{ \alpha }\) (see Sect. 3, Eq. (18)), and

  2. (ii)

    \(a > \uplambda L\), the diameter of support of \(g_{ \uplambda }( \vec { x })\) (see Sects. 2 and 4).

Since the parameter L is scaled with \(\uplambda \), its initial value is of no particular importance, and without restricting generality we can put

$$\begin{aligned} L = \frac{ 2 \pi ^{ 2 } }{ \alpha }, \end{aligned}$$
(63)

which leaves one dimensional parameter in the model. With this identification the set of admitted scaling parameters is \(\uplambda \in ( 0, 1 ]\). However, for the sake of simplification of the following discussion we restrict the set to \(\uplambda \in ( 0, 1 / 2 )\).

The replacement of \(g( \vec { x })\) by \(g_{ \uplambda }( \vec { x })\) in formulas for \(\mathcal {E}_{ 0 }\) and \(\mathcal {E}_{ \text {int} }( a )\) amounts to the use of functions \(M_{ \uplambda p }\), \(t_{ \uplambda }( ik )\) and \(u_{ \uplambda }( i k )\) (see Sect. 3) in place of \(M_{ p }\), t(ik) and u(ik) in (61) and (62). While the a-independent term \(\mathcal {E}_{ 0 }( \uplambda )\) does not contribute to the Casimir force, and for that reason we omit it from consideration, it is worth noting that it is divergent in the limit. This should be expected, as it describes twice the self-energy of \(\delta ( \vec { x })\). It will prove convenient to introduce dimensionless integration variables \(q = \mathrm{pa}\), \(l =\mathrm{ka} \) and to denote

$$\begin{aligned} r_{ \uplambda }( l )&:= \frac{ \gamma }{ \alpha } t_{ \uplambda }( \tfrac{ il }{ a } ) = \gamma + \frac{ 1 }{ \pi } \int \limits _{ -\infty }^{ +\infty } \mathrm{d} u \, \frac{ l^{ 2 } }{ l^{ 2 } + u^{ 2 } } M_{ \uplambda u / a }, \end{aligned}$$
(64a)
$$\begin{aligned} n_{ \uplambda }( l )&:= \frac{ \gamma }{ \alpha } u_{ \uplambda }( \tfrac{ il }{ a } ) = M_{ i \uplambda l / a } e^{ -l }. \end{aligned}$$
(64b)

The interaction energy takes now the form

$$\begin{aligned} \begin{aligned} \mathcal {E}_{ \text {int} }( a, \uplambda )&= \frac{ 4 \alpha }{ \pi ^{ 4 } \gamma }\, \mathfrak {Re}\int \limits _{ \mathbb {R}_{ + }^{ 2 } } \mathrm{d}^{ } q \, \mathrm{d}^{ } l \; \frac{ q^{ 2 } }{ ( q + il )^{ 2 } } M_{ \uplambda q / a } \left( Q_{ \uplambda }^{ ( 1 ) }( a, l ) + h( q ) Q_{ \uplambda }^{ ( 2 ) }( a, l ) \right) . \end{aligned}\nonumber \\ \end{aligned}$$
(65)

We introduce here useful notation

$$\begin{aligned} Q_{ \uplambda }^{ ( 1 ) }( a, l ) = \frac{ n_{ \uplambda }( l )^{ 2 } }{ r_{ \uplambda }( l ) [ r_{ \uplambda }( l )^{ 2 } - n_{ \uplambda }( l )^{ 2 } ] }, \qquad Q_{ \uplambda }^{ ( 2 ) }( a, l ) = \frac{ n_{ \uplambda }( l ) }{ r_{ \uplambda }( l )^{ 2 } - n_{ \uplambda }( l )^{ 2 } }. \end{aligned}$$
(66)

We now want to find the asymptotic expansion of this energy for \(\uplambda \searrow 0\). We shall prove below that the result is

$$\begin{aligned}&\mathcal {E}_{ \text {int} }( a, \uplambda )\nonumber \\&\quad = \frac{ 2 \alpha }{ \pi ^{ 3 } } \Bigg [ \frac{ \chi }{ \uplambda } \int \limits _{ 0 }^{ +\infty } \frac{ e^{ -2 l }\, dl }{ ( \gamma + l ) [ ( \gamma + l )^{ 2 } - e^{ -2 l } ] } + \frac{ b_{ 1 } \chi }{ \gamma } \int \limits _{ 0 }^{ +\infty } \frac{ l^{ 2 } [ 3 ( \gamma + l )^{ 2 } e^{ -2 l } - e^{ -4 l } ] }{ ( \gamma + l )^{ 2 } [ ( \gamma + l )^{ 2 } - e^{ -2 l } ]^{ 2 } } \, \mathrm{d}l \nonumber \\&\qquad - \frac{ 2 }{ \gamma } \int \limits _{ 0 }^{ +\infty } \frac{ l e^{ -2 l }\, dl }{ ( \gamma + l ) [ ( \gamma + l )^{ 2 } - e^{ -2 l } ] } + \frac{ 1 }{ \gamma } \int \limits _{ 0 }^{ +\infty } \frac{ ( 1 - l ) e^{ -2 l } }{ ( \gamma + l )^{ 2 } - e^{ -2 l } }\, \mathrm{d}l + R( a, \uplambda ) \Bigg ],\nonumber \\ \end{aligned}$$
(67)

where the remainder fulfills the bound

$$\begin{aligned} \begin{aligned} |{ R( a, \uplambda ) }| \le \mathrm {const}( a )\, \uplambda . \end{aligned} \end{aligned}$$
(68)

Constant \(b_{ 1 } > 0\) depends on the function \(g( \vec { x })\), its precise meaning is explained in “Appendix B.1”, Eq. (128). The strategy for the proof of this formula is as follows.

  1. (a)

    Obtain the leading terms of the asymptotic expansion in \(\uplambda \) of the integrals over q in formula (65); estimate the rest.

  2. (b)

    Find leading terms of the Taylor expansion in \(\uplambda \) of \(Q_{ \uplambda }^{ ( i ) }( a, l )\); estimate the rest.

  3. (c)

    Combine these results to obtain result (67) and (68).

Technical details of (a) and (b) are contained in “Appendices C and D”.

The result for (a) is

$$\begin{aligned}&\mathcal {I}_{ 1 } = \mathfrak {Re}\int \limits _{ 0 }^{ +\infty } \mathrm{d} q \, \frac{ q^{ 2 } }{ ( q + il )^{ 2 } } M_{ \uplambda q / a } = \frac{ 1 }{ \uplambda } \, \frac{ \pi \gamma \chi }{ 2 } - \pi l + \mathcal {R}_{ 1 }( a, l, \uplambda ), \end{aligned}$$
(69a)
$$\begin{aligned}&|{ \mathcal {R}_{ 1 }( a, l, \uplambda ) }| \le \uplambda \frac{ 3 \pi }{ 2 \gamma } l^{ 2 }, \end{aligned}$$
(69b)
$$\begin{aligned}&\mathcal {I}_{ 2 } = \mathfrak {Re}\, \int \limits _{ 0 }^{ +\infty } \mathrm{d} q \, \frac{ q \sin ( q ) }{ ( q + il )^{ 2 } } M_{ \uplambda q / a } = \frac{ \pi }{ 2 } ( 1 - l ) e^{ -l } + \mathcal {R}_{ 2 }( a, l, \uplambda ), \end{aligned}$$
(70a)
$$\begin{aligned}&|{ \mathcal {R}_{ 2 }( a, l, \uplambda ) }| \le \uplambda \frac{ \pi }{ 2 \gamma } l ( 2 + l ) e^{ -l / 2 }. \end{aligned}$$
(70b)

For (b) one obtains

$$\begin{aligned}&Q_{ \uplambda }^{ ( 1 ) }( a, l ) = \frac{ e^{ -2 l } }{ ( \gamma + l ) [ ( \gamma + l )^{ 2 } - e^{ -2 l } ] } + \mathcal {R}_{ 3 }( a, l, \uplambda ), \end{aligned}$$
(71)
$$\begin{aligned}&\begin{aligned} Q_{ \uplambda }^{ ( 1 ) }( a, l )&= \frac{ e^{ -2 l } }{ ( \gamma + l ) [ ( \gamma + l )^{ 2 } - e^{ -2 l } ] }\\&\quad + \uplambda \, \frac{ b_{ 1 } }{ \gamma } \, \frac{ l^{ 2 } [ 3 ( \gamma + l )^{ 2 } e^{ -2 l } - e^{ -4 l } ] }{ ( \gamma + l )^{ 2 } [ ( \gamma + l )^{ 2 } - e^{ -2 l } ]^{ 2 } } + \mathcal {R}_{ 4 }( a, l, \uplambda ), \end{aligned} \end{aligned}$$
(72)
$$\begin{aligned}&Q^{ ( 2 ) }_{ \uplambda }( a, l ) = \frac{ e^{ -l } }{ ( \gamma + l )^{ 2 } - e^{ -2 l } } + \mathcal {R}_{ 5 }( a, l, \uplambda ), \end{aligned}$$
(73)
$$\begin{aligned}&|{ \mathcal {R}_{ 3 }( a, l, \uplambda ) }| \le \uplambda \frac{ l e^{ -l } d }{ \gamma ^{ 2 } } \left\{ 1 + \frac{ l }{ \gamma } + d \left[ 2l ( \gamma + l ) + e^{ -l } \right] \right\} , \end{aligned}$$
(74)
$$\begin{aligned}&|{ \mathcal {R}_{ 4 }( a, l, \uplambda ) }| \nonumber \\&\quad \le \uplambda ^{ 2 } \frac{ l e^{ -l } d }{ \gamma ^{ 3 } } \bigg ( l d \bigg \{ l ( \gamma + l ) + l^{ 2 } + \frac{ 13 }{ 24 } e^{ -l } + 4 d \left[ l ( \gamma + l ) + \tfrac{ 1 }{ 2 } e^{ -l } \right] ^{ 2 } \bigg \} \nonumber \\&\qquad + \frac{ ( \gamma + l )^{ 2 } }{ \gamma } \left\{ d \left[ 2 l ( \gamma + l ) + e^{ -l } \right] + \frac{ l }{ \gamma ^{ 2 } } \right\} + \frac{ l }{ 24 \gamma } ( 13 \gamma + 25l ) \bigg ), \end{aligned}$$
(75)
$$\begin{aligned}&|{ \mathcal {R}_{ 5 }( a, l, \uplambda ) }| \le \uplambda \frac{ l e^{ -l / 2 } d }{ \gamma } \left\{ \frac{ 1 }{ 2 } + d \left[ 2 l ( \gamma + l ) + e^{ -l } \right] \right\} . \end{aligned}$$
(76)

We use here for the sake of brevity function \(d( a, l, \uplambda )\) defined in “Appendix D” and we also suppress it arguments.

Step (c) in the proof of the expansion is now easily achieved, with the details of estimation of the rest given by the formulas:

$$\begin{aligned} R( a, \uplambda ) = R_{ 1 }( a, \uplambda ) + R_{ 2 }( a, \uplambda ) + R_{ 3 }( a, \uplambda ) + R_{ 4 }( a, \uplambda ) + R_{ 5 }( a, \uplambda ), \end{aligned}$$
(77)

with

$$\begin{aligned} R_{ 1 }( a, \uplambda )= & {} \frac{ 4 }{ \pi ^{ 4 } \gamma } \int \limits _{ 0 }^{ +\infty } \mathrm{d} l \, \mathcal {R}_{ 1 }( a, l, \uplambda ) Q^{ ( 1 ) }_{ \uplambda }( a, l ), \end{aligned}$$
(78)
$$\begin{aligned} |{ R_{ 1 }( a, \uplambda ) }|\le & {} \uplambda \frac{ 6 }{ \pi ^{ 3 } \gamma ^{ 2 } } \int \limits _{ 0 }^{ +\infty } \mathrm{d} l \, l^{ 2 } e^{ -l } d, \end{aligned}$$
(79)
$$\begin{aligned} R_{ 2 }( a, \uplambda )= & {} \frac{ 1 }{ \uplambda } \, \frac{ 2 \chi }{ \pi ^{ 3 } } \int \limits _{ 0 }^{ +\infty } \mathrm{d} l \, \mathcal {R}_{ 4 }( a, l, \uplambda ), \end{aligned}$$
(80)
$$\begin{aligned} |{ R_{ 2 }( a, \uplambda ) }|\le & {} \uplambda \frac{ 2 \chi }{ \pi ^{ 3 } \gamma ^{ 3 } } \int \limits _{ 0 }^{ +\infty } \mathrm{d} l \, l e^{ -l } \mathrm{d} \Bigg ( l d \Bigg \{ l ( \gamma + l ) + l^{ 2 } + \frac{ 13 }{ 24 } e^{ -l } \nonumber \\&+ 4 d \left[ l ( \gamma + l ) + \frac{ 1 }{ 2 } e^{ -l } \right] ^{ 2 } \Bigg \} \nonumber \\&+ \frac{ ( \gamma + l )^{ 2 } }{ \gamma } \left\{ \mathrm{d} \left[ 2 l ( \gamma + l ) + e^{ -l } \right] + \frac{ l }{ \gamma ^{ 2 } } \right\} + \frac{ 1 }{ 24 } \, \frac{ l }{ \gamma } ( 13 \gamma + 25 l ) \Bigg ), \end{aligned}$$
(81)
$$\begin{aligned} R_{ 3 }( a, \uplambda )= & {} -\frac{ 4 }{ \pi ^{ 3 } \gamma } \int \limits _{ 0 }^{ +\infty } \mathrm{d} l \, l \mathcal {R}_{ 3 }( a, l, \uplambda ), \end{aligned}$$
(82)
$$\begin{aligned} |{ R_{ 3 }( a, \uplambda ) }|\le & {} \uplambda \frac{ 4 }{ \pi ^{ 3 } \gamma ^{ 3 } } \int \limits _{ 0 }^{ +\infty } \mathrm{d} l \, l^{ 2 } e^{ -l } \mathrm{d} \left\{ 1 + \frac{ l }{ \gamma } + \mathrm{d} \left[ 2 l ( \gamma + l ) + e^{ -l } \right] \right\} , \end{aligned}$$
(83)
$$\begin{aligned} R_{ 4 }( a, \uplambda )= & {} \frac{ 4 }{ \pi ^{ 3 } \gamma } \int \limits _{ 0 }^{ +\infty } \mathrm{d} l \, \mathcal {R}_{ 2 }( a, l, \uplambda ) Q_{ \uplambda }^{ ( 2 ) }( a, l ), \end{aligned}$$
(84)
$$\begin{aligned} |{ R_{ 4 }( a, \uplambda ) }|\le & {} \uplambda \frac{ 2 }{ \pi ^{ 3 } \gamma ^{ 2 } } \int \limits _{ 0 }^{ +\infty } \mathrm{d} l \, l ( 2 + l ) e^{ -l } d, \end{aligned}$$
(85)
$$\begin{aligned} R_{ 5 }( a, \uplambda )= & {} \frac{ 2 }{ \pi ^{ 3 } \gamma } \int \limits _{ 0 }^{ +\infty } \mathrm{d} l \, \mathcal {R}_{ 5 }( a, l, \uplambda ) ( 1 - l ) e^{ -l }, \end{aligned}$$
(86)
$$\begin{aligned} |{ R_{ 5 }( a, l, \uplambda ) }|\le & {} \uplambda \frac{ 2 }{ \pi ^{ 3 } \gamma ^{ 2 } } \int \limits _{ 0 }^{ +\infty } \mathrm{d} l \, e^{ -3 l / 2 } l ( 1 + l ) \mathrm{d} \left\{ \frac{ 1 }{ 2 } + \mathrm{d} \left[ 2 l ( \gamma + l ) + e^{ -l } \right] \right\} . \end{aligned}$$
(87)

To attain more information about the remainder, one should use a bound for d, e.g. as given by the formula (147a). We shall not pursue this point further in details, we only note that for large \(\gamma \) we can use inequality (147b), obtaining

$$\begin{aligned} |{ R( a, \uplambda ) }| \le \uplambda \sum _{ i = 1 }^{ N } C_{ i } \frac{ 1 }{ \gamma ^{ m_{ i } } } \, \frac{ 1 }{ ( \gamma ^{ 2 } - 1 )^{ n_{ i } } }, \quad m_{ i }, n_{ i } \in \mathbb {N}_{ + }, \end{aligned}$$
(88)

for some constants N and \(C_{ i }\). This result shows that the remainder \(R( a, \uplambda )\) vanishes with \(\uplambda \), when \(\gamma \) is fixed, which is equivalent to fixed distance a, and that with \(\uplambda \) fixed it vanishes when \(a \rightarrow +\infty \).

8 Local Energy Density

In this section we want to discuss the local energy density in the state \(\Omega _{ \vec { a }}\) in our model. For the class of models described in Sect. 1, the point-splitting procedure gives the general formula for this quantity:

$$\begin{aligned} \mathcal {E}( \vec { a }, \vec { x }) = T_{ \vec { a }}( \vec { x }, \vec { x }), \end{aligned}$$
(89)

where \(T_{ \vec { a }}( \vec { x }, \vec { y })\) is the kernel of the distribution \(T_{ \vec { a }}( \varphi , \psi )\) defined by

$$\begin{aligned} T_{ \vec { a }}( \varphi , \psi ) = \frac{ 1 }{ 4 } ( \varphi , ( h_{ \vec { a }} - h ) \psi ) + \frac{ 1 }{ 4 } ( \nabla \varphi , ( h_{ \vec { a }}^{ -1 } - h^{ -1 } ) \nabla \psi ), \end{aligned}$$
(90)

where the second term is understood as the scalar product of two gradients. Functions \(\varphi ( \vec { x })\) and \(\psi ( \vec { x })\) are chosen as real with compact support. The above formulas, in a somewhat inexplicit form, were derived in [10] (see formulas (6.8)–(6.10) with (5.13) of this reference). In the form stated above they were given in [11] (see (73), (74) in this reference).

To find \(T_{ \vec { a }}( \vec { x }, \vec { y })\) we follow approach from [11]. First, we note two elementary identities

$$\begin{aligned} \frac{ 1 }{ k } = \frac{ 2 }{ \pi } \int \limits _{ 0 }^{ +\infty } \mathrm{d} r \, \frac{ 1 }{ r^{ 2 } + k^{ 2 } }, \quad \; k - p = -\frac{ 2 }{ \pi } \int \limits _{ 0 }^{ +\infty } \mathrm{d} r \, r^{ 2 } \left( \frac{ 1 }{ r^{ 2 } + k^{ 2 } } - \frac{ 1 }{ r^{ 2 } + p^{ 2 } } \right) . \end{aligned}$$
(91)

By spectral representation of \(h_{ \vec { a }}\) and h, they allow us to write the scalar products from (90) as

$$\begin{aligned} ( \varphi , ( h_{ \vec { a }} - h ) \psi )= & {} -\frac{ 2 }{ \pi } \int \limits _{ 0 }^{ +\infty } \mathrm{d} r \, r^{ 2 } ( \varphi , [ G( -r^{ 2 } ) - G_{ 0 }( -r^{ 2 } ) ] \psi ), \end{aligned}$$
(92a)
$$\begin{aligned} (\nabla \varphi , ( h_{ \vec { a }}^{ -1 } - h^{ -1 } ) \nabla \psi )= & {} \frac{ 2 }{ \pi } \int \limits _{ 0 }^{ +\infty } \mathrm{d} r \, ( \nabla \varphi , [ G( -r^{ 2 } ) - G_{ 0 }( -r^{ 2 } ) ] \nabla \psi ), \end{aligned}$$
(92b)

where the integral kernel in scalar products is given by (see (27)):

$$\begin{aligned}&{[} G( -r^{ 2 } ) - G_{ 0 }( -r^{ 2 } ) ]( \vec { x }, \vec { y }) \nonumber \\&\quad = \frac{ 1 }{ \sqrt{ 2 \pi }^{ 6 } } \int \limits _{ \mathbb {R}^{ 3 }} d^{ 3 } p\, e^{ i \vec { p }\cdot \vec { x }} \frac{ \mathcal {F}_{ \vec { p }} }{ \vec { p }^{ \; 2 } + r^{ 2 } }\, \mathcal {T}( i r ) \int _{ \mathbb {R}^{ 3 }} d^{ 3 } q\, e^{ -i \vec { y }\cdot \vec { q }} \frac{ \mathcal {F}_{ \vec { q }}^{ \; \dag } }{ \vec { q }^{ \; 2 } + r^{ 2 } }. \end{aligned}$$
(93)

This can be much simplified if we use notation

$$\begin{aligned} J( r, \vec { x })= & {} \left( H( r, |{\, \vec { x }- \vec { b } \,}| ) \qquad H( r, |{\, \vec { x }+ \vec { b } \,}| ) \right) , \end{aligned}$$
(94a)
$$\begin{aligned} H( r, |{\, \vec { x } \,}| )= & {} \frac{ 1 }{ \sqrt{ 2\pi }^{ 3 } } \int \limits _{ \mathbb {R}^{ 3 }} d^{ 3 }p\; \widehat{g}( p ) \frac{ e^{ i \vec { p }\cdot \vec { x }} }{ \vec { p }^{ \; 2 } + r^{ 2 } }. \end{aligned}$$
(94b)

The function \(H( r, |{ \vec { x } }| )\) and its derivative in the second argument \(H'( r, |\, \vec { x }\, | )\) have obvious bounds

$$\begin{aligned} |{ H( r, |{\, \vec { x } \,}| ) }| \le \text {const} \, \frac{ 1 }{ 1 + r^{ 2 } }, \qquad |{ H'( r, |{\, \vec { x } \,}| ) }| \le \text {const} \, \frac{ 1 }{ 1 + r^{ 2 } }. \end{aligned}$$
(95)

Using (94a) we write

$$\begin{aligned}{}[ G( -r^{ 2 } ) - G_{ 0 }( -r^{ 2 } ) ]( \vec { x }, \vec { y }) = J( r, \vec { x }) \mathcal {T}( i r ) J^{ T }( r, \vec { y }). \end{aligned}$$
(96)

The integral \(T_{ \vec { a }}( \vec { x }, \vec { y })\) can be now easily found by integration by parts:

$$\begin{aligned} T_{ \vec { a }}( \vec { x }, \vec { y })= & {} -\frac{ 1 }{ 2\pi } \int \limits _{ 0 }^{ +\infty } \mathrm{d} r \, r^{ 2 } J( r, \vec { x }) \mathcal {T}( i r ) J^{ T }( r, \vec { y }) \nonumber \\&+ \frac{ 1 }{ 2\pi } \int \limits _{ 0 }^{ +\infty } \mathrm{d} r \, \sum \limits _{ j = 1 }^{ 3 } \nabla _{ j } J( r, \vec { x }) \mathcal {T}( i r ) \nabla _{ j } J^{ T }( r, \vec { y }), \end{aligned}$$
(97)

where the convergence of integrals follows from (95) and the bound on matrix \(\mathcal {T}( w )\) given in Eq. (25).

After putting \(\vec { x }= \vec { y }\) we arrive at the expression for the energy density

$$\begin{aligned} \mathcal {E}( \vec { a }, \vec { x })= & {} -\frac{ 1 }{ 2 \pi } \int \limits _{ 0 }^{ +\infty } \mathrm{d} r \, \frac{ t( ir ) }{ t( ir )^{ 2 } - u( ir )^{ 2 } } [ H'( r, |\, \vec { x }+ \vec { b }\, | )^{ 2 } - r^{ 2 } H( r, |\, \vec { x }+ \vec { b }\, | )^{ 2 } \nonumber \\&+ H'( r, |\, \vec { x }- \vec { b }\, | )^{ 2 } - r^{ 2 } H( r, |\, \vec { x }- \vec { b }\, | )^{ 2 } ] \nonumber \\&+\frac{ 1 }{ \pi } \int \limits _{ 0 }^{ +\infty } \mathrm{d} r \, \frac{ u( ir ) }{ t( ir )^{ 2 } - u( ir )^{ 2 } } \bigg [ r^{ 2 } H( r, |\, \vec { x }+ \vec { b }\, | ) H( r, |\, \vec { x }- \vec { b }\, | ) \nonumber \\&- \frac{ ( \vec { x }+ \vec { b }\, ) \cdot ( \vec { x }- \vec { b }\, ) }{ |\, \vec { x }+ \vec { b }\, | |\, \vec { x }- \vec { b }\, | } H'( r, |\, \vec { x }+ \vec { b }\, | ) H'( r, |\, \vec { x }- \vec { b }\, | ) \bigg ], \end{aligned}$$
(98)

a quantity well-defined for all \(\vec { x }\).

Similarly as in the global energy case, we now want to exclude the self-energy of the two interaction centers. For this end we put \(\vec { x }- \vec { b }= \vec { y }\) (this shifts one of the centers to the origin) and then take the limit \(a \rightarrow \infty \). In this way we find

$$\begin{aligned} \mathcal {E}_{ \text {single} }( \vec { y }) = \lim _{ a \rightarrow \infty } \mathcal {E}( \vec { a }, \vec { b }+ \vec { y }) = -\frac{ 1 }{ 2 \pi } \int \limits _{ 0 }^{ +\infty } \mathrm{d} r \, \frac{ 1 }{ t( ir ) } [ H'( r, |\, \vec { y }\, | )^{ 2 } - r^{ 2 } H( r, |\,\vec { y }\, | )^{ 2 } ].\nonumber \\ \end{aligned}$$
(99)

This formula describes energy density around single interaction center in the whole space. We now subtract from the total energy density the self-energy of the two centers

$$\begin{aligned} \begin{aligned} \mathcal {E}_{ \text {self} }( \vec { a }, \vec { x })&= \mathcal {E}_{ \text {single} }( \vec { x }+ \vec { b }\, ) + \mathcal {E}_{ \text { single } }( \vec { x }- \vec { b }\, ), \end{aligned} \end{aligned}$$
(100)

which results in interaction energy

$$\begin{aligned} \mathcal {E}_{ \text {int} }( \vec { a }, \vec { x })= & {} -\frac{ 1 }{ 2 \pi } \int \limits _{ 0 }^{ +\infty } \mathrm{d} r \, \frac{ u( ir )^{ 2 } }{ t( ir ) [ t( ir )^{ 2 } - u( ir )^{ 2 } ] } [ H'( r, |\, \vec { x }+ \vec { b }\, | )^{ 2 } - r^{ 2 } H( r, |\, \vec { x }+ \vec { b }\, | )^{ 2 } \nonumber \\&+ H'( r, |\, \vec { x }- \vec { b }\, | )^{ 2 } - r^{ 2 } H( r, |\, \vec { x }- \vec { b }\, | )^{ 2 } ] \nonumber \\&+ \frac{ 1 }{ \pi } \int \limits _{ 0 }^{ +\infty } \mathrm{d} r \, \frac{ u( ir ) }{ t( ir )^{ 2 } - u( ir )^{ 2 } } \bigg [ r^{ 2 } H( r, |\, \vec { x }+ \vec { b }\, | ) H( r, |\, \vec { x }- \vec { b }\, | ) \nonumber \\&- \frac{ ( \vec { x }+ \vec { b }\, ) \cdot ( \vec { x }- \vec { b }\, ) }{ |\, \vec { x }+ \vec { b }\, | |\, \vec { x }- \vec { b }\, | } H'( r, |\, \vec { x }+ \vec { b }\, | ) H'( r, |\, \vec { x }- \vec { b }\, | ) \bigg ]. \end{aligned}$$
(101)

8.1 Local Energy Density in the Scaling Limit

We now turn to the investigation of the energy density when the function \(g( \vec { x })\) is replaced by \(\widehat{g}_{ \uplambda }( \vec { x })\) and the limit \(\uplambda \searrow 0\) is performed. More precisely, we consider \(\mathcal {E}( \vec { a }, \vec { x }, \uplambda )\) with fixed a and \(\vec { x }\) and compute this quantity in aforementioned limit. This means that in (100) and (101) we need to replace t(ir), u(ir) by \(t_{ \uplambda }( ir )\), \(u_{ \uplambda }( ir )\) (see Sect. 4) and \(H( r, |\, \vec { x }\, | )\) by

$$\begin{aligned} H_{ \uplambda }( r, |\, \vec { x }\, | ) = \frac{ 1 }{ \sqrt{ 2\pi }^{ 3 } } \int \limits _{ \mathbb {R}^{ 3 }} d^{ 3 }p\; \widehat{g}( \uplambda p ) \frac{ e^{ i \vec { p }\cdot \vec { x }} }{ \vec { p }^{ \, 2 } + r^{ 2 } }. \end{aligned}$$
(102)

Result of this is the following

$$\begin{aligned}&\mathcal {E}_{ \text {self} }( \vec { a }, \vec { x }, \uplambda ) = -\frac{ 1 }{ 2 \pi } \int \limits _{ 0 }^{ +\infty } \mathrm{d} r \, \frac{ 1 }{ t_{ \uplambda }( ir ) } [ H_{ \uplambda }'( r, |\, \vec { x }+ \vec { b }\, | )^{ 2 } - r^{ 2 } H_{ \uplambda }( r, |\, \vec { x }+ \vec { b }\, | )^{ 2 } ] \nonumber \\&\quad - \frac{ 1 }{ 2 \pi } \int \limits _{ 0 }^{ +\infty } \mathrm{d} r \, \frac{ 1 }{ t_{ \uplambda }( ir ) } [ H_{ \uplambda }'( r, |\, \vec { x }- \vec { b }\, | )^{ 2 } - r^{ 2 } H_{ \uplambda }( r, |\, \vec { x }- \vec { b }\, | )^{ 2 } ]. \end{aligned}$$
(103)
$$\begin{aligned}&\mathcal {E}_{ \text {int} }( \vec { a }, \vec { x }, \uplambda ) = -\frac{ 1 }{ 2 \pi } \int \limits _{ 0 }^{ +\infty } \mathrm{d} r \, \frac{ u_{ \uplambda }( ir )^{ 2 } }{ t_{ \uplambda }( ir ) [ t_{ \uplambda }( ir )^{ 2 } - u_{ \uplambda }( ir )^{ 2 } ] } \nonumber \\&\quad \times \big [ H_{ \uplambda }'( r, |\, \vec { x }+ \vec { b }\, | )^{ 2 } - r^{ 2 } H_{ \uplambda }( r, |\, \vec { x }+ \vec { b }\, | )^{ 2 } + H_{ \uplambda }'( r, |\, \vec { x }- \vec { b }\, | )^{ 2 } - r^{ 2 } H_{ \uplambda }( r, |\, \vec { x }- \vec { b }\, | )^{ 2 } \big ] \nonumber \\&\quad + \frac{ 1 }{ \pi } \int \limits _{ 0 }^{ +\infty } \mathrm{d} r \, \frac{ u_{ \uplambda }( ir ) }{ t_{ \uplambda }( ir )^{ 2 } - u_{ \uplambda }( ir )^{ 2 } } \bigg [ r^{ 2 } H_{ \uplambda }( r, |\, \vec { x }+ \vec { b }\, | ) H_{ \uplambda }( r, |\, \vec { x }- \vec { b }\, | ) \nonumber \\&\quad - \frac{ ( \vec { x }+ \vec { b }\, ) \cdot ( \vec { x }- \vec { b }\, ) }{ |\, \vec { x }+ \vec { b }\, | |\, \vec { x }- \vec { b }\, | } H_{ \uplambda }'( r, |\, \vec { x }+ \vec { b }\, | ) H_{ \uplambda }'( r, |\, \vec { x }- \vec { b }\, | ) \bigg ]. \end{aligned}$$
(104)

We now want to investigate the limit of this expressions when \(\uplambda \searrow 0\). Inequality (111) shows that for any fixed \(\vec { x }\), if \(\uplambda \in \big ( 0, \min ( 1/2, \alpha |\, \vec { x }\, | / \pi ^{ 2 } ) \big )\), the integral defining \(H( r, |\, \vec { x }\, | )\) can be computed using the residuum technique and the result is

$$\begin{aligned} H_{ \uplambda }( r, |{\, \vec { x } \,}| ) = \sqrt{ \frac{ \pi }{ 2 } } \widehat{g}( i \uplambda r ) \frac{ e^{ -r |{\, \vec { x } \,}| } }{ |{\, \vec { x } \,}| }, \quad H_{ \uplambda }'( r, |\, \vec { x }\, | ) = -\sqrt{ \frac{ \pi }{ 2 } } \widehat{g}( i \uplambda r ) ( 1 + r |{\, \vec { x } \,}| ) \frac{ e^{ -r |{\, \vec { x } \,}| } }{ |{\, \vec { x } \,}|^{ 2 } }.\nonumber \\ \end{aligned}$$
(105)

This also shows that we can use the dominated convergence theorem to obtain the limit \(\uplambda \searrow 0\):

$$\begin{aligned}&\mathcal {E}_{ \text {single} }( \vec { x }, 0 ) = -\frac{ 1 }{ 2 } \int \limits _{ 0 }^{ +\infty } \mathrm{d} r \, \frac{ 1 }{ \alpha + 2 \pi ^{ 2 } r } \frac{ e^{ -2 r |{\, \vec { x } \,}| } }{ |\, \vec { x }\,|^{ 4 } } \left( 1 + 2 r |{\, \vec { x } \,}| \right) , \end{aligned}$$
(106)
$$\begin{aligned}&\mathcal {E}_{ \text {self} }( \vec { a }, \vec { x }, 0 ) = \mathcal {E}_{ \text {single} }( \vec { x }+ \vec { b }, 0 ) + \mathcal {E}_{ \text {single} }( \vec { x }- \vec { b }, 0 ), \end{aligned}$$
(107)
$$\begin{aligned}&\mathcal {E}_{ \text {int} }( \vec { a }, \vec { x }, 0 ) = -\frac{ 1 }{ 8\pi ^{ 2 } } \int \limits _{ 0 }^{ +\infty } \mathrm{d} l \,\frac{ e^{ -2 l } }{ ( \gamma + l ) [ ( \gamma + l )^{ 2 } - e^{ -2 l } ] } \nonumber \\&\times \Bigg [ \frac{ e^{ -2 l |{\, \vec { x }+ \vec { b } \,}| / a } }{ |\, \vec { x }+ \vec { b }\, |^{ 4 } } \left( 1 + 2 l \frac{ |\, \vec { x }+ \vec { b }\, | }{ a } \right) + \frac{ e^{ -2 l |{\, \vec { x }- \vec { b } \,}| / a } }{ |\, \vec { x }- \vec { b }\, |^{ 4 } } \left( 1 + 2 l \frac{ |{\, \vec { x }- \vec { b } \,}| }{ a } \right) \Bigg ] \nonumber \\&+ \frac{ 1 }{ 4 \pi ^{ 2 } } \int \limits _{ 0 }^{ +\infty } \mathrm{d} l \, \frac{ e^{ -l } }{ ( \gamma + l )^{ 2 } - e^{ -2 l } } \frac{ e^{ -l ( |\, \vec { x }+ \vec { b }\, | + |\, \vec { x }- \vec { b }\, | ) / a } }{ |\, \vec { x }+ \vec { b }\, | |\, \vec { x }- \vec { b }\, | } \Bigg [ \frac{ l^{ 2 } }{ a^{ 2 } } \Bigg ( 1 - \frac{ ( \vec { x }+ \vec { a }) \cdot ( \vec { x }- \vec { a }) }{ |\, \vec { x }+ \vec { b }\, | |\, \vec { x }- \vec { b }\, | } \Bigg )\nonumber \\&- \frac{ ( \vec { x }+ \vec { b }\, ) \cdot ( \vec { x }- \vec { b }\, ) }{ |{\, \vec { x }+ \vec { b } \,}| |{\, \vec { x }- \vec { b } \,}| } \frac{ 1 + l ( |{\, \vec { x }+ \vec { b } \,}| + |{\, \vec { x }- \vec { b } \,}| ) / a }{ |{\, \vec { x }+ \vec { b } \,}| |{\, \vec { x }- \vec { b } \,}| } \Bigg ]. \end{aligned}$$
(108)

For the sake of clarity we have written the last integral in dimensionless integration variable \(l = ra\). We also point out that due to the nature of \(\delta \)’s interactions, the above expressions for local energy density are meaningful only if the conditions \(\vec { x }\ne \pm \vec { b }\) are satisfied.

9 Discussion of Results

In this work, we have investigated the Casimir effect in a well-defined model of the scalar field \(\phi ( \vec { x })\) with two nonsingular interaction centers in three space dimensions. With appropriate assumptions the evolution generator is positive and has purely absolutely continuous spectrum equal to \([ 0, +\infty )\). We have shown that the Casimir energy and the number of particles given by (4a) and (4b) are finite. We have also considered a \(\uplambda \)-rescaled model which for \(\uplambda \searrow 0\) tends to the model of two Dirac \(\delta \)-s in the sense of the resolvent limit of operators. In this limit the Casimir energy diverges, but we derived its asymptotic expansion (67) in \(\uplambda \). Analogous result was achieved for the local density of energy in vacuum state \(\Omega _{ \vec { a }}\).

The asymptotic expansion of the global energy (67) is the central result of the article. As our model is well-posed only when \(\chi > 0\) and \(b_{ 1 }>0\), we conclude that there is no universal (independent of the approximating model) expression for the Casimir force in the setting considered here. The first term in this expression is obviously positive and decreasing in \(\gamma \). Numerical analysis shows that the same is true for the second and the sum of the third and fourth terms. Therefore, the force is repulsive in all considered models.

In consequence, our analysis applies to systems which may very closely approximate strict point interaction, but shows that the limit model is not a correct idealization for the consideration of the backreaction force. Precise prediction for the force depends on particular values of \(\chi \) and \(b_{ 1 }\), but it is reasonable to think that in vicinity of the limit model the term proportional to \(1 / \uplambda \) will be dominating all others.

These conclusions are in full concord with the results obtained earlier, within the approach used here, for a different system: scalar and electromagnetic fields in interaction with two parallel plates. Admissible regularization of this system in momentum space was considered by Herdegen [10], and the configuration space regularization was used later by Herdegen and Stopa [11]; our analysis follows the latter technique. Each of these two regularization methods resulted in the same conclusion: for the Neumann boundary conditions in addition to the universal term in the backreaction energy per area \(-\pi ^{ 2 } / 1440 a^{ 3 }\) predicted by Casimir, there appears a term of the form \(c/\uplambda a^2\), which diverges in the limit of strict boundary conditions.

In the light of that development, our work further confirms the hypothesis that on theoretical level the Casimir force containing only simple universal terms is the result of unjustified regularization. However, as already mentioned in Sect. 1, experiments are involved in many complications of physical and technical nature. In our view, results of these experiments do not exclude existence of non-universal terms in Casimir force. On the other hand, the approach used here must be developed further, to admit more realistic settings and finite temperature. Such development, at this moment, remains an open task.

Beside the global energy, we have also considered the local energy density. For \(\vec { x }\ne \pm \vec { b }\) (outside the position of the limit \(\delta \)-centers) this density has a well-defined limit for \(\uplambda \searrow 0\), given by (108). This energy density is totally independent of the choice of a particular model satisfying our assumptions, and as a result one cannot recover the formula for the global Casimir energy integrating it over whole space. This means that significant part of the global energy is concentrated at the position of interaction centers. Moreover, the integral is in fact not well-defined due to the terms of order \(1 / |{\, \vec { x }\pm \vec { b } \,}|^{ 4 }\) in formula (108).

We need to stress that although such integration of energy density is mathematically dubious, Herdegen and Stopa found that in the case of two plates it yields for the scalar field the canonical expression \(-\pi ^{ 2 } / 1440 a^{ 3 }\). This difference may be due to greater singularity (three dimensional) of the present model, and in this way shows that the finite result mentioned above is rather exceptional.

We now want to compare our results with those already found in literature. Scardicchio [17] uses a variant of the zero-point energy approach to consider a class of systems, among them the scalar field interacting with the two-deltas system. His predictions for such system are radically different from ours, giving a universal (regularization-independent), attractive force. This again shows, that the system considered in the present article supplies a stronger test for the comparison of various approaches. General remarks on the comparison of zero-point idea with the present approach may be found in [10].

In their work Spreafico and Zerbini [18] compute the partition function for the scalar field interacting with one or two deltas in finite temperature. Using this, one can derive, in appropriate limit, a formula for their prediction of the Casimir energy of this system in zero temperature. This expression is very different both from our result, as well as that of Scardicchio. The authors do not analyze (nor even write explicitly) their formula, which turns out to be difficult to interpret and would need a numerical analysis.

It is worth pointing out, that there do exists some similarities of the form of integrals appearing in all three predictions of the Casimir energy mentioned above. However, these similarities are quite superficial, and their source is easy do understand. All three works make use, although in quite different ways, of the resolvent of the operator “\(-\Delta + \mu \delta ( \vec { x }+ \vec { b }\, ) + \mu \delta ( \vec { x }- \vec { b }\, )\)”. The differences in predictions should not be a surprise, as all three approaches are based on different ideas, and understand “Casimir energy” differently.

We are not aware of an alternative local energy formula for the system of two \(\delta \)s. However, Fermi and Pizzocchero [9] have a formula for the density in case of a single \(\delta \), which should be juxtaposed with our expression (106). Their expression depends on a parameter \(\xi \) responsible for a conformal addition to the energy–momentum tensor (see authors’ discusion of this topic in appendix A in [8]). In the flat spacetime context this is a nonstandard modification; in the present article we work with the standard energy density for \(\xi = 0\) (see section 6 in [10] for more details).

Our formula for energy density is very different from that of [9] with \(\xi = 0\). In fact, for any value of \(\xi \) these two formulas do not coincide. We believe that the reason for this is as follows. In our work we use a quasi-potential, whose support is compact. Outside that support the energy density is well defined without any regularization. The limit \(\uplambda \searrow 0\) shrinks the support to a point, outside which the density is uniquely defined. On the other hand, Fermi and Pizzacchero start by regularizing the field using a highly nonlocal operator dependent on parameter u, and apply the (seemingly) local prescription for energy density to this nonlocal field (in effect, this is not a local density). After some process of regularization and renormalization they proceed to the limit \(u \rightarrow 0\), which is supposed to restore the local field. We are not aware of any arguments why this circuitous procedure should yield back the correct local expression for energy.

When Fermi and Pizzacchero try to find the global Casimir energy by integrating their \(u = 0\) limit energy density over \(\mathbb {R}^{ 3 } \setminus \{ 0 \}\), they find their formula diverges. A cure of this deficiency is attempted by Fermi in [7], where he reverts the order of operations. He first integrates the regularized, u-depended energy density, and only then attempts the limit \(u \rightarrow 0\). As it turns out, this does not solve the problem and the energy is still infinite. The author concludes that this confirms the expectation that the delta model is a bad idealization of real physical phenomena at small scale.

This conclusion is in concord with the results of the algebraic analysis of singular potentials, as discussed in articles by Herdegen and our present work. Also, one can easily check that the self energy (61) diverges in the scaling limit. However, we believe that our analysis helps to better understand the source of difficulties.

We would like to end this discussion by stressing once more the difference of the approach adopted here from most of other treatments present in literature, including the articles mentioned above. These works take for consideration models with singular interaction (this singularity may be precisely identified in algebraic language, as mentioned in Introduction). Secondly, popular prescriptions for “Casimir energy”, like e.g. zero mode summation, are also by their very nature singular. These singularities need some regularization/renormalization to yield some finite result. The nature of these operations, although often mathematically ingenious, is physically not transparent. On the other hand, the present approach takes into consideration models well defined in one common representation for changing parameters, and identifies the Casimir energy as a well defined expectation value for all values of those parameters. In the present work the model was that of non-singular quasi-potentials. It is then a mathematically legitimate, and physically wholly justified problem to investigate the behaviour of the model and energy in singular limit of point interactions. The limit itself is not a justified idealization, but the asymptotic expansion gives predictions for acceptable systems “not far” from idealized physical settings.

Numerical computations in this work were done using numpy v. 1.16.1 package in Python 3.6.8. The source code for this computation can be acquired by writing to the authors on email address: kziemianfvt@gmail.com.