Casimir interactions between two short-range Coulomb sources

Annals of Physics 404 (2019) 132–157

Contents lists available at ScienceDirect

Annals of Physics journal homepage: www.elsevier.com/locate/aop

Casimir interactions between two short-range Coulomb sources ∗

Yu. Voronina a , I. Komissarov b , , K. Sveshnikov a a

Department of Physics and Institute of Theoretical Problems of MicroWorld, Moscow University, 119991, Leninsky Gory, Moscow, Russia b Department of Physics, Columbia University, New York, NY 10027, USA

article

info

Article history: Received 18 November 2018 Accepted 21 February 2019 Available online 1 March 2019 Keywords: One dimensional Dirac-Coulomb system Vacuum polarization Non-perturbative quantum electrodynamics Casimir interaction

a b s t r a c t We calculate the Casimir interaction between two short range extended charge sources, embedded in a background of one dimensional massive Dirac fermions. We explore both the induced polarization density ρv ac (x) and the corresponding integral charge Qv ac , as well as the Casimir energy E v ac and its contribution to the interaction between sources. For sources with the same charge, by means of the novel ln [Wronskian] contour integration techniques we find that the interaction energy between two sources can exceed sufficiently large negative values and simultaneously reveal the features of a long-range force in spite of non-zero fermion mass, which could significantly influence the properties of such quasi-one-dimensional systems. For large distances between sources we recover that their mutual interaction is governed first of all by the structure of the discrete spectrum of a single source, in dependence on which it can be tuned to give an attractive, a repulsive, or an (almost) compensated Casimir force with various rates of the exponential fall-down, quite different from the standard exp(−2ms) law. © 2019 Elsevier Inc. All rights reserved.

1. Introduction There is now a lot of interest to the study of various non-perturbative vacuum polarization effects in sufficiently strong static or adiabatically slowly varying Coulomb fields, generated by ∗ Corresponding author. E-mail addresses: [email protected] (Y. Voronina), [email protected] (I. Komissarov), [email protected] (K. Sveshnikov). https://doi.org/10.1016/j.aop.2019.02.014 0003-4916/© 2019 Elsevier Inc. All rights reserved.

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

133

localized extended charge sources (see e.g., [1–15] and refs. therein). Such effects have attracted a considerable amount of theoretical and experimental activity in 3 + 1 D heavy ions collisions, where for Z > Zcr ,1 ≃ 170 a non-perturbative reconstruction of the vacuum state is predicted, which should be accompanied by a number of nontrivial effects including the vacuum positron emission ([1–3,5,7,9–11] and refs. therein). Similar in essence effects are expected to take place both in 2 + 1 D (planar hetero-structures like graphene [12–14]) and in 1 + 1 D (one dimensional ‘‘hydrogen ion’’ [9–11,15]). In this work we explore another essentially non-perturbative effect, namely the vacuum polarization in quasi-one-dimensional QED caused by one or two short-range sources in the form of potential wells. The main question of interest is how the non-perturbative vacuum effects, including the effects of super-criticality, manifest themselves in Casimir forces between such sources. Physical examples of such systems, where the contribution from E int v ac should be quite significant, fill the range from a cluster of heavy ions with large Z [3,5,7], up to charged impurities in low-dimensional nanostructures and cold atoms systems [12–14,16–18]. It should be mentioned that the Casimir interaction caused by the fermion vacuum polarization has been intensively studied from different points of view and in different geometries [19]. In particular, in [20] the main focus was made on the study of the energy density and interaction between two ‘‘Dirac spikes’’ as a function of a single ‘‘spike’’ parameters and the distance between them. In this model each ‘‘spike’’ is represented by a square barrier, which enters the fermion dynamics as an additional mass term, considered in the δ -limit via transfer-matrix, that allows to imitate the boundary conditions for a fermionic bag. In [21,22] the Casimir interaction between two square potential barriers (‘‘scatterers’’), mediated by the massless fermions, has been considered. The Casimir force between the scatterers was found for both the case of finite width and strength of the barriers and in the δ -limit. In the case of identical scatterers, separated by a large distance, the interaction force between them reveals the conventional attractive asymptotics. At the same time, for a more general case of inequivalent scatterers it was shown that the magnitude and sign of the force depends on relative spinor polarizations of potential barriers. Unlike [20–22], in this paper we consider the Casimir interaction of two short-range Coulomb sources, which enter the fermion dynamics as the electrostatic potentials. In the case of the scalar coupling, considered in [20–22], the scatterers affect equally the positive- and negative-frequency fermionic modes. In the case of vector coupling the behavior of electronic and positronic components is principally different and leads to a number of new effects, the most significant of which is the discrete levels diving into the lower continuum and related non-perturbative effects of vacuum reconstruction, when the sources attain the overcritical region. In particular, the Casimir interaction energy between the sources with charges of the same sign can acquire large negative values and so significantly influence the dynamics of such quasi-one-dimensional systems. Moreover, in such systems the Casimir force shows up a nontrivial behavior with increasing distance between sources, when there become possible as separate vertical jumps caused by the effect of creation–annihilation of discrete levels at the lower threshold, as well as different exponent rates and signs in the decreasing law and also the features of a long-range force despite the non-zero fermion mass. The single short-range Coulomb source is described as a square well of width 2a and depth V0 V (x) = −V0 θ (a − |x|) .

(1)

Actually the potential (1) is nothing else but a projection of the Coulomb potential of a spherical shell of radius R0 and charge Z onto x-axis U(x) = −Z α (θ (R0 − |x|)/R0 + θ (|x| − R0 )/|x|) ,

(2)

strongly screened for |x| > R0 . In this work we consider the system of two such sources, separated by the distance s. The component of the vacuum polarization energy E int v ac , responsible for their interaction, is defined as E int v ac (s) = E v ac ,2 (s) − E v ac ,2 (s → ∞),

(3)

where E v ac ,2 (s) is the total vacuum polarization energy for the system containing two potentials like (1), located at the distance s from each other, while E v ac ,2 (s → ∞) = 2 E v ac with the latter being the vacuum energy of a single source.

134

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

Fig. 1. The ‘‘critical charges’’ Zcr ,1 , Zcr ,2 , Zcr ,3 for three lowest levels in the single well (1).

It would be worth to note that since we consider here the sources with several parameters (for a single well these are the depth V0 and the half-width a), the subcritical and overcritical regions for a concrete level are defined by a set of pairs (V0 , a), rather than by a single quantity Zcr , as it occurs whenever a concrete relation between the size and charge of the Coulomb source is implied. In the case of a single source (1) in the diagram (V0 , a) the subcritical and overcritical regions are separated by a curve (see Fig. 1). Therefore under the notion of a ‘‘critical charge’’ Zcr ,i for the ith discrete level we will imply the whole set of the source parameters, which separate the sub- and overcritical regions from each other, rather than one definite quantity. In particular, in Fig. 1 there are shown three first ‘‘critical charges’’ Zcr ,1 , Zcr ,2 , Zcr ,3 for the potential (1). As in other works on vacuum polarization in the strong Coulomb fields [23–26], the radiative corrections from virtual photons are neglected. Henceforth, if it is not stipulated separately, the units h ¯ = m = c = 1 are used. Thence the coupling constant α = e2 is also dimensionless, and the numerical calculations, illustrating the general picture, are performed for α = 1/137. However, it would be worthwhile to note that for the effective electron–hole vacuum in the quasi-onedimensional systems like nanotubes, as in graphene, the actual value of α , and hence, the magnitude of the Casimir effects could be quite different from the ones in the pure vacuum QED. 2. Evaluation of the induced density for a single short-range Coulomb source The non-perturbative approach to evaluation of the induced vacuum density ρv ac (x) and, in turn, of the vacuum energy E v ac is based on the Wichmann–Kroll (WK) method ([23–27] and refs. therein). Its application to quasi-one-dimensional Dirac–Coulomb systems has been considered in detail in Refs. [9–11]. So here we will discuss only those points, which are crucial for calculation of Casimir interaction between two sources like (1). The WK method reduces the calculation of the vacuum density to integration of the trace of the Green function TrG for the Dirac equation (DE) along the contours P(R) and E(R) on the first sheet of the Riemann energy surface (Fig. 2)

ρvac (x) = −

|e| 2

( lim

R→∞

1 2π i

∫

dϵ TrG(x, x; ϵ ) + P(R)

1 2π i

)

∫

dϵ TrG(x, x; ϵ ) .

(4)

E(R)

At the same time, TrG(x, x; ϵ ) can be represented as TrG(x, x; ϵ ) =

1 J(ϵ )

ψLT (x)ψR (x),

(5)

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

135

Fig. 2. Special contours in the complex energy plane, used for the representation of the vacuum charge density via contour integrals.

with ψL (x) and ψR (x) being regular at x = ±∞ solutions of the Dirac spectral problem [α p + β + V (x)] ψ (x) = ϵ ψ (x) ,

(6)

in which α = σ2 , β = σ3 are used, while J(ϵ ) is their Wronskian J(ϵ ) = ψL,2 (x)ψR,1 (x) − ψL,1 (x)ψR,2 (x) .

(7)

Since DE with the potential (1) is even under reflection, the functions ψL (x) and ψR (x) are connected via parity transformation ψL (x) = βψR (−x) and so can be written in the following form

ψR (x) = θ (x − a)Y (x) + θ (a − |x|)(AΦ (x) + BΨ (x)) + θ (−x − a)(CY (x) + DW (x)) , ψL (x) = θ (x − a)(CW (x) + DY (x)) + θ (a − |x|)(−AΦ (x) + BΨ (x)) + θ (−x − a)W (x) ,

(8)

where the coefficients A, B, C , D are found from the continuity conditions for ψL (x), ψR (x) at the points x = ±a, while Φ (x), Ψ (x), Y (x), W (x) represent linearly independent solutions of (6) in the corresponding regions of continuity of the potential (1), namely

) (√ V0 + ϵ + 1 sin [j(ϵ, V0 )x] √ , Φ (x) = V0 + ϵ − 1 cos [j(ϵ, V0 )x] (√ ) √V0 + ϵ + 1 cos[j(ϵ, V0 )x] , Ψ (x) = − V0 + ϵ − 1 sin[j(ϵ, V0 )x] (√ ) (√ ) 1 + ϵ e−γ (ϵ )x 1 + ϵ eγ (ϵ )x √ √ Y (x) = , W (x) = , 1 − ϵ eγ (ϵ )x − 1 − ϵ e−γ (ϵ )x

(9)

where j(ϵ, V0 ) =

√

(V0 + ϵ )2 − 1 ,

γ (ϵ ) =

√

1 − ϵ2 .

(10)

136

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

From the continuity conditions for ψL (x), ψR (x) one finds A=

D=

[Ψ (x), Y (x)]a

,

B=

[Y (x), Φ (x)]a

,

[Ψ (x), Φ (x)]a [Ψ (x), Φ (x)]a [Ψ (x), Y (x)]a [Φ (x), Y (x)]a , C =2 [Ψ (x), Φ (x)]a [Y (x), W (x)]a [Ψ (x), Y (x)]a [W (x), Φ (x)]a + [Y (x), Φ (x)]a [Ψ (x), W (x)]a [Ψ (x), Φ (x)]a [Y (x), W (x)]a

(11)

,

where [f (x), g(x)]a = f2 (a)g1 (a) − g2 (a)f1 (a) is the Wronskian of the functions f (x) and g(x) at x = a. Proceeding further, from the explicit form of ψL (x) and ψR (x) one finds the following expression for their Wronskian [Ψ (x), Y (x)]a [Φ (x), Y (x)]a J(ϵ ) = −2 . (12) [Ψ (x), Φ (x)]a Compiling the answers (5), (8)–(12), one finds the following expression for TrG(x, x; ϵ ) of DE with the potential (1) TrG(x, x; ϵ ) =

⎧ ( ) W T (x)W (x) [Φ (x), W (x)]a [Ψ (x), W (x)]a W T (x)Y (x) ⎪ ⎪ − + x < −a , ⎪ ⎪ ⎪ [Ψ (x), Y (x)]a ) ⎪ [W (x), Y (x)]a (2[W (x), Y (x)]a [Φ (x), Y (x)]a ⎨ 1 [Ψ (x), Y (x)]a T [Φ (x), Y (x)]a T Φ (x)Φ (x) − Ψ (x)Ψ (x) , |x| ⩽ a , ⎪ 2[Ψ (x), Φ (x)]a [Φ (x), Y (x)]a ( [Ψ (x), Y (x)]a ⎪ ) ⎪ T T ⎪ [Φ (x), W (x)]a Y (x)Y (x) [Ψ (x), W (x)]a W (x)Y (x) ⎪ ⎪ ⎩ , x>a. − + [W (x), Y (x)]a 2[W (x), Y (x)]a [Φ (x), Y (x)]a [Ψ (x), Y (x)]a

(13)

Substitution of the explicit form of solutions (9) into (13) yields the final expression for TrG in the region |x| ⩾ a TrG(x, x; ϵ ) = √

ϵ

1 − ϵ2 e2γ (ϵ )(a−|x|)

V0 sin[2j(ϵ, V0 )a] , (14) γ (ϵ ) j(ϵ, V0 )γ (ϵ ) cos[2j(ϵ, V0 )a] + [1 − ϵ (V0 + ϵ )] sin[2j(ϵ, V0 )a] whereas in the region |x| < a one obtains V0 cos[2j(ϵ, V0 )x] − (V0 + ϵ )[1 − ϵ (V0 + ϵ )] cos[2j(ϵ, V0 )a] 1 TrG(x, x; ϵ ) = j(ϵ, V0 ) j(ϵ, V0 )γ (ϵ ) cos[2j(ϵ, V0 )a] + [1 − ϵ (V0 + ϵ )] sin[2j(ϵ, V0 )a] (15) (V0 + ϵ )γ (ϵ ) sin[2j(ϵ, V0 )a] + . j(ϵ, V0 )γ (ϵ ) cos[2j(ϵ, V0 )a] + [1 − ϵ (V0 + ϵ )] sin[2j(ϵ, V0 )a]

+

In the next step we consider the asymptotics of TrG in the energy plane on the segments of the large circle C1 (R), C2 (R), C3 (R), C4 (R) (see Fig. 2). For the segments C1 (R) and C2 (R) (|ϵ| → ∞, 0 < Arg(ϵ ) < π ) in the upper half-plane one obtains

⎧ ⎪ ⎪ ⎨i + i + O(|ϵ|−4 ) , |x| > a , 2ϵ 2 TrG(x, x; ϵ ) → ⎪ ⎪ ⎩i + i − iV0 + O(|ϵ|−4 ) , |x| ⩽ a . 2ϵ 2 ϵ3

(16)

⎧ ⎪ ⎪ ⎨−i − i + O(|ϵ|−4 ) , |x| > a , 2ϵ 2 TrG(x, x; ϵ ) → ⎪ ⎪ i iV0 ⎩ −i − 2 + 3 + O(|ϵ|−4 ) , |x| ⩽ a . 2ϵ ϵ

(17)

In turn, in the lower half-plane for the segments C3 (R) and C4 (R) (|ϵ| → ∞, π < Arg(ϵ ) < 2π ) there follows

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

137

From (16) and (17) one finds that the integration contour in (4) can be reduced to the imaginary axis, whence there follows the final expression for the vacuum density in the non-perturbative approach

∫ +i∞ |e| dϵ TrG(x, x; ϵ ) 2π i −i∞ ∫ +∞ |e| dy TrG(x, x; iy) . = 2π −∞

ρvac (x) =

(18)

In the presence of negative discrete levels −1 ⩽ ϵn < 0 the expression (18) should be replaced by [9–11,24]

[ ψn† (x)ψn (x)

∑

ρvac (x) = |e|

1

+

2π

−1⩽ϵn <0

]

+∞

∫

dy TrG(x, x; iy)

.

(19)

−∞

Proceeding further, it is worth to note the following general properties of TrG(x, x; ϵ ) under the change of the sign of the external field V0 → −V0 and by the complex conjugation TrGV0 (x, x; ϵ ) = −TrG−V0 (x, x; −ϵ ) ,

(20)

TrGV0 (x, x; ϵ )∗ = TrGV0 (x, x; ϵ ∗ ) . Combining the properties (20), one finds TrGV0 (x, x, iy)∗ = −TrG−V0 (x, x; iy).

(21)

There follows from (20), (21) that Re[TrGV0 (x, x; iy)] is odd in V0 and even in y, whereas Im[TrGV0 (x, x; iy)] is even in V0 and odd in y. Hence, ρv ac (x) is determined only via Re[TrGV0 (x, x; iy)] and so is an odd in V0 real function. However, the expressions for the vacuum density (18) and (19) do not provide the condition of vanishing integral vacuum charge Qv ac in the subcritical region Z < Zcr . At the same time, in contrast with the vacuum charge, induced by the unscreened Coulomb potentials [9–11], in the present case the induced charge Qv ac is finite for the whole range of external source parameters a and V0 . The reason is that Re[TrG(x, x; iy)] converges exponentially in the argument x:

√ 1+y2 (a−|x|)

Re[TrG(x, x; iy)] ∼ e2

,

|x| → ∞ .

(22)

Moreover, the real part of the Green function Re [TrG(x, x; iy)] decreases exponentially in y for |y| → ∞ provided |x| > a and polynomially as O(1/y3 ) for |x| ⩽ a. According with [9–11,24], matching between ρv ac (x) and vanishing integral vacuum charge Qv ac is achieved via renormalization of the linear in the external potential V0 component of the induced density. Namely, the term in ρv ac (x) proportional to V0 should be replaced by the corresponding (1) renormalized expression ρv ac (x), obtained from the first-order perturbation theory in a quite similar way with the corresponding expressions found for the unscreened case in [9–11]. In the case under consideration this expression reads

√

ρ

(1) v ac (x)

=

2 α V0

π2

+∞

∫

dq 0

cos qx sin qa q

(

arcsinh(q/2)

1−2 √ q 1 + (q/2)2

) .

(23) (1)

Note that since the external potential Aext 0 (x) has finite jumps at x = ±a, ρv ac (x) is also discontinuous at the same points (see Fig. 3). (3+) Following Ref. [24], for these purposes one should first introduce the density ρv ac (x) as

[ ρ

(3+) v ac (x)

= |e|

∑ −1⩽ϵn <0

ψn† (x)ψn (x)

+

1 2π

∫

]

+∞

dy (TrG(x, x; iy) − TrG (x, x; iy)) (1)

−∞

,

(24)

138

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

(1)

Fig. 3. One-loop vacuum density ρv ac (x) for a = 1, V0 = 1.

where TrG(1) (x, x; ϵ ) is defined by the first Born approximation for the Green function in (6) and can be easily found analytically TrG(1) (x, x; iy) =

V0

[

(1 + y2 )3/2

√ 1+y2

e−2|x|

√

sinh(2a 1 + y2 ) θ (|x| − a)

√ ( ) ] √ 2 + 1 − e−2a 1+y cosh(2x 1 + y2 ) θ (a − |x|) .

(25)

(3+)

By construction, the density ρv ac (x) does not contain the linear in Z terms, since the subtracted term (3+) TrG(1) (x, x; ϵ ) is fully responsible for them. As a result, for Z < Zcr ,1 ρv ac (x) does not contribute to the integral vacuum charge. In absence of negative discrete levels, which prevent the analytic continuation, this statement can be easily proved via transition∫ to the complex x-plane, where ρv(3ac+) (x), represented as a uniformly converging integral (|e|/2π ) dy[TrG(x, x; iy) −∫TrG(1) (x, x; iy)], (3+) + turns out to be an analytic function of x. In this case the induced charge Qv ac = dx ρv3ac (x) can be expressed as a contour integral along the half-circle in the upper half-plane (Im x > 0), which (3+) vanishes exactly, since ρv ac (x) decreases exponentially for |x| → ∞. Thus, the renormalized density ρvRac (x) should be defined as (3+) ρvRac (x) = ρv(1) ac (x) + ρv ac (x) ,

(26)

(1)

where ρv ac (x) is the perturbative renormalized density (23). This expression for ρvRac (x) provides vanishing of the vacuum charge in the subcritical region Z < Zcr ,1 . In presence of negative discrete levels vanishing of Qv ac for Z < Zcr ,1 follows from the model-independent arguments, based on the initial expression for the vacuum density [9–11]. ren A more detailed picture of changes in ρVP (x) is given in Refs. [1–3,5,7], based on the U. Fano approach [28]. The main result is that, when the discrete level ψn (x) dives into the lower continuum, the change in the vacuum density equals to

⏐ ∆ρVP (x) = −|e| ψn (x)† ψn (x)⏐ϵ

n =−1

.

(27)

However, it should be noted that this approach contains a number of approximations, and so the expression (27) is exact only in the vicinity of the corresponding Zcr , what is explicitly shown in Refs. [9–11]. The reliable evaluation of ρvRac (x) for the whole range of the source parameters should be based on expressions (24) and (26) with subsequent check of the integer value for Qv ac . The induced charge densities, calculated in this way for the potential well (1) with the half-width a = 1, are shown in Figs. 4, 6 and 7. In Fig. 4 the renormalized charge density is shown for the purely perturbative case corresponding to V0 = 1, thereon for V0 = 2.85, when the first critical potential for the given half-width V0,cr1 ≃ 2.8621 is not reached yet and further for V0 = 2.88, when the

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

139

Fig. 4. ρvRac (x) for the source (1) with a = 1 and V0 = 1, 2.86, 2.87, 4.29, 4.3 (color online).

first (even) discrete level has already dived into the lower continuum. Thereafter there are shown the densities for V0 = 4.28, when V0,cr2 ≃ 4.2969 is not reached yet, and, finally, for V0 = 4.31, just after the second (odd) level has dived into the lower continuum. The direct numerical check confirms that the integral vacuum charge for V0 = 1 and 2.85 vanishes exactly, Qv ac = −|e| for V0 = 2.88 and 4.28, and Qv ac = −2|e| for V0 = 4.31. The corresponding ‘‘critical charges’’ for both even and odd levels are defined from the equation

√

sin(2a (V0 − 1)2 − 1) = 0 .

(28)

Here it should be noted that, as it follows from Fig. 4, the difference in the induced densities ρvRac (x) for the cases V0 = 2.85, 2.88 and V0 = 4.28, 4.31 looks quite small, although their integral values in both cases differ by |e|. The origin of this circumstance is that the contributions from negative † discrete levels |e| ψn (x)ψn (x), which approach the lower threshold, are highly dispersed along the x-axis and tend to the uniform distribution for ϵ → −1. In Fig. 5 the evolution of the charge density of the first (even) discrete level by approaching the first critical value a = 1 , V0,cr1 ≃ 2.8621 is shown. Such behavior of the induced density differs significantly from the corresponding behavior of the electronic WF at the lower threshold for the unscreened case, when the charge distribution of the discrete level by approaching Zcr is well localized. The reason is that in contrast to the purely Coulomb problem, in DE with short-range source like (1) the state with ϵF = −1 is the starting one of the continuous spectrum with vanishing wave number. It can be assumed that such a regime of discrete levels diving into the lower continuum (namely, absence of overcritical metastable states) could lead to absence of Fano resonances (that means, appearance of sharp jumps by π in the positron scattering phase). These resonances has been much discussed in Refs. [1–3,5,7] in connection with possible positronic peaks in heavy ion collisions, which should signal the achievement of the supercritical region and non-perturbative QED vacuum reconstruction. Actually, the calculations confirm that the positronic scattering phase does not reveal the corresponding resonance behavior, and the poles originating from the discrete levels dived into the lower continuum on the second sheet of the Riemann energy surface in this case are absent. Nevertheless, this circumstance does not pose any serious problems for the self-consistent treatment of the considered picture of the vacuum shells formation due to dived discrete levels. As it was shown in Ref. [2,29], the resonances, which originate from such levels, show up now in a different fashion and can be easily detected via indirect methods (e.g., through the analysis of the positronic part of the transition coefficient T+ ). As it follows from Fig. 6, the main changes in the profile of the function ρvRac (x) with increasing V0 occurs when one of the discrete levels approaches the lower threshold (the vacuum shell, created after diving of the first even level into the lower continuum, is localized within the potential well).

140

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

Fig. 5. The evolution of the charge density of the first (even) discrete level by approaching the well’s half-width to the first critical value a = 1 , V0,cr1 ≃ 2.8621.

Fig. 6. The evolution of ρvRac (x) with increasing depth of the potential well V0 in the range V0,cr1 < V0 < V0,cr2 for a = 1 (color online).

In Fig. 7 the evolution of the induced density is shown for increasing well’s width with constant depth V0 = 5. The solid line corresponds to the purely perturbative case with a = 0.4, the dashed one to a = 0.65, when the first critical half-width of the well for the given potential acr ,1 ≃ 0.4055 is already crossed and the lowest (even) level has reached the lower continuum, while the last (dotted) one corresponds to a = 0.9, when the second critical half-width acr ,2 = 2acr ,1 ≃ 0.8111 is exceeded and the next lowest (odd) level has dived into the lower continuum. The direct numerical calculation confirms that the corresponding induced charges equal to Qv ac = 0, −|e|, −2|e|, as expected.

3. Evaluation of Casimir energy via ln[Wronskian] contour integration Non-perturbative polarization effects, which in the overcritical region show up in ρv ac (x) via localized vacuum shells formation, lead also to corresponding changes in the vacuum energy E v ac .

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

141

Fig. 7. ρvRac (x) V0 = 5 a = 0.4, 0.65, 0.9.

Quite similar to emergence of the vacuum shells, these changes turn out to be essentially nonperturbative and so cannot be restored directly from ρv ac (x) [10]. It should be noted that the shells formation rate and increase of their contribution to E v ac depend strongly on the number of spatial dimensions. In 1 + 1 D this contribution is minimal, thence the behavior of E v ac in the non-perturbative region depends merely on the renormalization term and on the decrease of the contribution from the perturbative term in E v ac by increasing V0 and/or a, rather than on the number of levels dived into the lower continuum. The starting point for the non-perturbative evaluation of the vacuum energy is the Schwinger vacuum average [1–3,5,7]

E v ac =

1

(

2

) ∑

ϵn −

ϵn <ϵF

∑

ϵn

ϵn ⩾ϵF

−

1

(

) ∑

2

ϵn −

ϵn <ϵF

A

∑ ϵn ⩾ϵF

,

ϵn

(29)

0

where ϵF = −1, the label A indicates the external field, while the label 0 corresponds to the free case. For the subsequent analysis it is convenient to separate in (29) the contributions from the discrete spectrum and both continua and apply to the latter the well-known tool, which replaces it by the scattering phase δ (ϵ ) integration (see e.g., [20,30,31] and refs. therein). After a number of almost evident steps one obtains [10]

E v ac =

1 2π

+∞

∫

δtot (ϵ ) dϵ + 1

1 2

∑

(1 − ϵn ) ,

(30)

−1⩽ϵn <1

where δtot (ϵ ) is the total phase shift for the given |ϵ|, including the contributions from scattering states in both continua, while in the discrete spectrum the (effective) electron rest mass is subtracted from each level in order to retain in E v ac the interaction effects only. Such approach to calculation of E VP turns out to be quite effective, since the total phase shift δtot (ϵ ) behaves in both (IR and UV) limits much better, than each of the elastic phases separately, and is automatically an even function of the external potential. More concretely, in the Coulomb potentials with non-vanishing source size δtot (ϵ ) is finite for ϵ → 1 and decreases ∼ 1/ϵ 3 for ϵ → ∞, what provides the convergence of the phase integral in (30) [9–11,13,14]. The sum over bound energies of discrete levels 1 − ϵn in the case of short-range sources like (1) turns out to

142

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

be finite from the very beginning, since such potentials allow for a finite number of discrete levels only. As a result, the expression (30) is finite without any additional renormalization. However, the convergence of E v ac in the form (30) is completely due to the specifics of 1 + 1 D and in no way means no need for a renormalization. The need for renormalization through the fermionic loop is based on the analysis of ρv ac (x), from which there follows that without such UV-renormalization the integral induced charge will not have the value that follows from obvious physical arguments [1–3,5,7,9–11,13,14,27]. Another obvious condition on E v ac is that for V0 → 0 (1) the value of E v ac should coincide with the expression for E v ac , obtained within the first order PT quite similar to the unscreened case considered in [9–11] E (1) v ac

=

2V02

+∞

∫

π2

dq 0

sin2 qa

(

arcsinh(q/2)

1−2 √ q 1 + (q/2)2

q2

) .

(31)

It is easy to verify that the non-renormalized vacuum energy (30) does not satisfy the latter requirement, since the introduced below renormalization coefficient (33), which provides also the (1) correspondence between E Rv ac and E v ac for V0 → 0, in the general case turns out to be non-zero. For these reasons, in complete analogy with ρv ac (x) we should pass from E v ac to the renormalized vacuum energy E Rv ac . Actually, this procedure is equivalent to a finite renormalization and normalization conditions as known from perturbative QED (see, e.g., [5,32]). In the practically useful form E Rv ac should be represented as follows E Rv ac = E v ac + λV02 ,

(32)

where (1)

λ = lim

E v ac (V0 ) − E v ac (V0 )

V02

V0 →0

.

(33)

It should be also noted that the renormalization coefficient λ depends solely on the shape of the external potential and in the general case is a sign-alternating function with zeros (see below). The evaluation of E Rv ac via the sum of the phase integral and discrete levels is considered in detail in Refs. [9–11,13,14] for the unscreened Coulomb potential and in the present case will differ only by certain minor technical details. At the same time, there exists an alternative approach for E v ac , similar to calculation of vacuum density ρv ac (x) via integration of specially constructed function along the WK contours (Fig. 2). It is easy to see that the function F (ϵ, V0 ) =

ϵ (dJ(ϵ )/dϵ) , J(ϵ )

(34)

where J(ϵ ) is the Wronskian for the spectral problem (6), reveals all the pole properties, which are required for the representation of the expression (29) via integrals along the WK contours, since actually J(ϵ ) is nothing else, but the Jost function of the spectral problem (6) [10]. The real zeros of J(ϵ ) lie in the interval −1 ⩽ ϵn < 1 and coincide with discrete energy levels, while the complex ones reside on the √ second sheet of the Riemann energy surface with negative imaginary part of the wavenumber k = ϵ 2 − 1 and for Re k > 0 correspond to the elastic resonances. Moreover, both J(ϵ ) and TrG as functions of the wavenumber k reveal the same reflection symmetry f ∗ (k) = f (−k∗ ) of the Jost function. To represent E v ac via integration along the WK contours, it suffices to pass to the difference E (ϵ, V0 ) = F (ϵ, V0 ) − F (ϵ, 0), normalized on the free case. As a result, the induced vacuum energy can be represented as E v ac = −

(∫

1 4π i

dϵ E (ϵ, V0 ) +

lim

R→∞

P(R)

∫

dϵ E (ϵ, V0 )

)

.

(35)

E(R)

The asymptotics of the function E (ϵ, V0 ) on the upper half of the large circle (|ϵ| → ∞, 0 < Arg(ϵ ) < π ) reads E (ϵ, V0 ) → −

V0

ϵ

+

V0 (2ia + V0 )

ϵ2

+ O(1/|ϵ|3 ) ,

(36)

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

143

while on the lower one (|ϵ| → ∞, π < Arg(ϵ ) < 2π ) one finds E (ϵ, V0 ) → −

V0

V0 (−2ia + V0 )

+

ϵ

ϵ2

+ O(1/|ϵ|3 ).

(37)

The first terms in the asymptotics (36) and (37) compensate each other due to the opposite direction of integration on the left and right half-circles, while the others disappear for R → ∞. Upon taking into account the (possible) negative discrete levels, in complete analogy with the result for the vacuum density [9–11,24] one finds the final expression for E v ac in the following form E v ac =

+∞

∫

1 2π

dy E (iy, V0 ) − −∞

∑

ϵn ,

(38)

−1⩽ϵn <0

where for the potential (1) the integrand in (38) takes the form E (iy, V0 ) =

iV0 y (V0 [V0 + 2iy]γ (iy) + 2aj2 (iy, V0 )γ 2 (iy) sin[2aj(iy, V0 )])

, V0 )γ 3 (iy)(j(iy, V0 )γ (iy) cos[2aj(iy, V0 )] + (1 − iV0 y + y2 ) sin[2aj(iy, V0 )]) 2ia V0 y j(iy, V0 ) γ 3 (iy) cos[2aj(iy, V0 )] −2 . j (iy, V0 )γ 3 (iy) (j(iy, V0 )γ (iy) cos[2aj(iy, V0 )] + (1 − iV0 y + y2 ) sin[2aj(iy, V0 )]) j2 (iy

(39) The direct calculation shows that both approaches to the vacuum energy (30) and (38) lead with a high precision to the same results. It would be worthwhile noticing here that the transformation of the initial WK-contours or, equivalently, the Feynman one CF , specially shown in Fig. 2, to the imaginary axis, employed above, is quite similar to the Wick rotation, which is a well-known tool in the theory of relativistic bound states (see, e.g.,[32]). The behavior of the renormalization term in (32) is determined by the coefficient λ(a), which can be represented as

λ(a) = λ1 (a) − λ2 (a) ,

(40)

where λ1 emerges from the perturbative contribution

λ1 (a) =

a

− I1 (a) , ∫ ∞

π 4

I1 (a) =

dq

π2

sin2 (qa)

(

arcsinh(q/2)

to the vacuum energy

) ,

1−2 √ q 1 + (q/2)2

q2

0

(1) E v ac

(41)

while λ2 corresponds to the first (quadratic) in V0 term in E v ac , which is found from the Born series (see [9–11,13,14,24])

√ λ2 (a) =

a

π

−

1 16

+ I2 (a) ,

I2 (a) =

1 4π

∞

∫

dy 0

e−4a

1+y2

(1 + y2 )2

.

(42)

By means of the relation λ1 (a) + λ2 (a) = a/π , which derivation requires for some additional considerations and so will be presented separately [33], the integral I1 (a) in (41) can be expressed via I2 (a), what gives

λ(a) =

a

π

− 2λ2 (a) =

1 8

−

a

π

− 2I2 (a).

(43)

Let us consider now the asymptotics of λ(a) for a ≪ 1 and a ≫ 1. For a ≪ 1 the main contribution to I2 (a) comes from the region y ∈ [−1/(2a), 1/(2a)], whereas the integral over the remaining region turns out to be O(a3 ). Expanding the exponent in the vicinity of a = 0 and calculating the integrals, one obtains I2 (a → 0) =

1 16

−

a

π

+ a2 + O(a3 ) .

(44)

144

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

Fig. 8. The dependence of the renormalization coefficient λ on the half-width of the potential well a.

Inserting the latter in (43), one finds the asymptotics of λ(a) for a ≪ 1 a λ(a → 0) = − 2a2 + O(a3 ) .

π

(45)

To find the asymptotics for a ≫ 1, one should firstly perform in the integral I2 (a) the change of variables y = sinh t and thereafter take the third derivative with respect to a. Proceeding further, by means of the integral representation of the MacDonald function it is possible to restore I2 (a) from its third derivative, what gives I2 (a) =

24

π

∞

∫

∫

′

∞

da

da a′

a

′′

∫

∞

a′′

da′′′ K0 (4a′′′ ).

(46)

Inserting the asymptotics of the MacDonald function for the large argument

√ K0 (z → ∞) =

π 2z

e−z

(

( )) 1+O

1 z

(47)

in (46) and performing a triple integration over a, one obtains the function with exponential falldown for growing a. So I2 (a ≫ 1) vanishes exponentially, what allows to write the asymptotics of the coefficient λ(a) for large a neglecting the exponentially small corrections as

λ(a → ∞) =

1 8

−

a

π

.

(48)

As a result, for small a the coefficient λ(a) grows linearly, while for large a it decreases with the same modulus slope 1/π . Hence, the renormalization coefficient λ(a) should vanish at certain a = acr . In the considered case of the well (1) one has acr ≃ 0.297. The shape of λ(a) on the most significant interval 0 < a < 1, where it turns out to be sign-alternating, is given in Fig. 8. Now let us consider the role of λ(a) in the behavior of E Rv ac (V0 ). For a < acr the renormalization coefficient yields the positive growing contribution ∼ V02 . In the case a > acr the renormalization coefficient is negative, hence the corresponding term λV02 turns out to be rapidly decreasing and shifts E Rv ac (V0 ) into negative region. When a ≃ acr , the value of E Rv ac (V0 ) is determined by competition between the sum over discrete levels and the phase integral. The performed calculations of E Rv ac (V0 ) for various values of the half-width a, when the renormalization coefficient (33) acquires all the three most representative values — positive, (almost) vanishing and negative, confirm the above mentioned general arguments. In particular, in Fig. 9a– d the curves E Rv ac (V0 ) for a = 0.1, acr , 1, 10 are shown. The number of vacuum shells for the given values of a is approximated as 0.06 V01.01 , 0.16 V01.04 , 0.39 V01.15 , 0.17 V03.92 , correspondingly. For a = 0.1 the renormalization coefficient λ ≃ 0.017, hence E Rv ac increases quadratically with V0 . For

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

145

Fig. 9. The dependence E Rv ac (V0 ) for several most representative values of the half-width of the well a: (a) a = 0.1 ; (b) a = acr ≃ 0.297; (c) a = 1; (d) a = 10.

a = acr the contribution of the renormalization term vanishes, and so the behavior of the vacuum energy with growing V0 in the overcritical region is determined merely by the contribution of the discrete spectrum, which increases ∼ 0.18 V01.01 . For a = 1, 10 the renormalization coefficients are equal to λ ≃ −0.195 and λ ≃ −3.06, respectively, while the non-renormalized vacuum energy in these cases grows as ∼ 0.73 V00.98 and ∼ 3.8 V01.64 . Hence, now with growing V0 the dominating contribution to E Rv ac (V0 ) will be determined by the renormalization term λV02 . As a result, E Rv ac (V0 ) becomes a rapidly decreasing function, which with growing V0 can acquire large negative values. Meanwhile, in the region of sufficiently small V0 the renormalized vacuum energy retains still the quadratic growth due to perturbative effects (see Fig. 9c,d). Let us consider now separately the contributions from the phase integral, discrete spectrum and the renormalization term to E Rv ac (V0 ) for the most representative cases a = acr and a = 10. For a = acr the behavior of separate components of E Rv ac is shown in Fig. 10. In Fig. 10a the behavior of the phase integral is presented. Fig. 10b demonstrates the behavior of the nonrenormalized vacuum energy E v ac (solid line), the total bound energy of discrete levels (the dashed one), the phase integral (denoted by dots) and, finally, the number of vacuum shells. In Fig. 10c the solid line shows the behavior of the renormalized vacuum energy E Rv ac (V0 ), while the dashed one corresponds to its approximation with account for the jumps by (−1), without which it is described by the function ∼ 0.22 V00.97 . Fig. 11 provides the same information for the case a = 10. The approximation function in Fig. 11c is ∼ −0.0008 V08.2 without jumps taken into account when plotting the curve. It should be noted also that the number of vacuum shells for a = 10 grows more rapidly, than for a = acr .

146

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

Fig. 10. (a): The phase integral; (b): the non-renormalized E v ac (V0 ), the sum of discrete levels bound energy, phase integral and the number of vacuum shells; (c): E Rv ac (V0 ) and its approximation for a = acr .

4. Casimir energy and vacuum density for two short-range Coulomb sources

Now – having dealt with the induced vacuum properties of the single source this way – let us turn to the configuration of such identical short-range Coulomb sources of the form V2 (x) = −V0 θ (|x| − d) θ (d + a − |x|) .

(49)

Let us note that now the separate sources have the total width a, that provides the restoration of the initial potential well (1) with the width 2a for d → 0. The perturbative induced density and polarization energy in this case are equal to

√

ρ

(1) v ac ,2

=

2 α V0

π2

+∞

∫

dq

cos qx(sin[q(a + d)] − sin[qd]) q

0

(

arcsinh(q/2)

1−2 √ q 1 + (q/2)2

) (50)

and (1) E v ac ,2

=

2V02

π2

+∞

∫

dq 0

[sin(q(a + d)) − sin (qd)]2 q2

(

arcsinh(q/2)

1−2 √ q 1 + (q/2)2

) ,

(51)

correspondingly. It is easy to see that in the limit d → 0 (50) coincides with (23), while (51) — with (31). The trace of the Green function is constructed through the expression (5), where the regular for x → −∞ and x → +∞ solutions of DE (6) for the potential (49) take the form

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

147

Fig. 11. (a): The phase integral; (b): the non-renormalized E v ac (V0 ), the sum of discrete levels bound energy, phase integral and the number of vacuum shells; (c): E Rv ac (V0 ) and its approximation for a = 10.

ψR (x) = θ (x − d − a) W (x) + θ (d ⩽ x < d + a) (−AΦ (x) + BΨ (x)) + θ (d − |x|) (CW (x) + DY (x)) + θ (−d − a ⩽ x < −d) (−E Φ (x) + F Ψ (x)) + θ (−x − d − a) (GW (x) + HY (x)) , ψL (x) = θ (x − d − a) (HW (x) + GY (x))

(52)

+ θ (d ⩽ x < d + a) (E Φ (x) + F Ψ (x)) + θ (d − |x|) (CY (x) + DW (x)) + θ (−d − a ⩽ x < −d) (AΦ (x) + BΨ (x)) + θ (−x − d − a) Y (x) , the functions Ψ (x), Φ (x), W (x), Y (x) are defined as in (9), while the coefficients A, B, C , D, E, F , G, H are found from the continuity conditions for ψL (x) and ψR (x) at the points of discontinuity of the potential (49). Further procedure of calculation and renormalization of the vacuum density and energy repeats completely the case of the single source and so does not need any special discussion besides two points. The first is that the explicit form of the integrand E (iy, V0 ) in the analog of (38) for the two-well configuration takes more than one page and so will be omitted. The second one is the structure of the renormalization term in E Rv ac ,2 . As in the case of one potential well, the calculation

148

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

of E Rv ac ,2 requires the renormalization in the second order with respect to the external potential E Rv ac ,2 = E v ac ,2 + Λ(a, d)V02 ,

(53)

where

Λ(a, d) = Λ1 (a, d) − Λ2 (a, d) .

(54)

The components of the renormalization coefficient Λ1 (a, d) and Λ2 (a, d), which play the same role as λ1 (a) and λ2 (a) for the single source, are expressed via the latter in the following form

Λ1 (a, d) = λ1 (a + d) + λ1 (d) + 2λ1 (a/2) − 2λ1 (d + a/2) , Λ2 (a, d) = λ2 (a + d) + λ2 (d) + 2λ2 (a/2) − 2λ2 (d + a/2) .

(55)

Combining both parts of (55) and using the identity λ1 (a) + λ2 (a) = a/π , one finds that the components of the renormalization coefficient (55) are subject of the same identity

Λ1 (a, d) + Λ2 (a, d) = a/π .

(56)

As a result, the renormalization coefficient for the two-well problem (49) can be represented as

Λ(a, d) = a/π − 2Λ2 (a, d)

= a/π − 2λ2 (a + d) − 2λ2 (d) − 4λ2 (a/2) + 4λ2 (d + a/2) .

(57)

Now let us list the results for found this way ρvRac ,2 and E Rv ac ,2 , which are necessary for the further analysis of the Casimir interaction between separate sources. The most significant here is the dependence of these quantities on the parameter d , 0 ≤ d ≤ ∞, which defines the distance between the sources in such a way that the distance s between the centers of the wells is given by s = a + 2d .

(58)

At first let us explore the features of the discrete spectrum of DE with the potential (49). For d → ∞ the wells become independent, while the spectrum of DE — twice degenerate. More concretely, with growing d the even levels increase, while the odd ones, in contrast, decrease, and so for d → ∞ the neighboring even and odd levels seek each other. To analyze the role of d in the overcritical region the equations for the critical parameters of the source (49) (i.e., the set [V0 , a , d], for which the discrete levels approach the lower continuum) should be considered. For odd levels the ‘‘critical charges’’ are determined from the equation

√

sin[a (V0 − 1)2 − 1] = 0 ,

(59)

which coincides with the similar equation for a single potential well (28) up the replacement a → 2a. Since (59) does not depend on d, any change of d for fixed (V0 , a) does not yield the diving of odd levels into the continuum. At the same time, for even levels the equation for their ‘‘critical charges’’ takes the form

√

√

√

(V0 − 1)2 − 1 cos[a (V0 − 1)2 − 1] + 2d V0 sin[a (V0 − 1)2 − 1] = 0 .

(60)

So for even levels the ‘‘critical charges’’ depend on the distance between the sources. The parameter d can be easily found from (60), and so the dependence d(V0 , a) together with condition d > 0 defines the critical values of d for even levels in the potential (49). The regions (V0 , a), where the even levels diving into the lower continuum takes place by certain d > 0, are shown as shaded ones in Fig. 12. The non-shaded regions in Fig. 12 correspond to those sets of (V0 , a), when the Eq. (60) does not possess any solutions with d > 0. The solid and dashed curves in Fig. 12 determine the values (V0 , a), when dcr = 0, and so correspond to the critical charges of a separate potential well (compare with Fig. 1). This way there appear two essentially different regimes for behavior of ρvRac ,2 (x) and E Rv ac ,2 as functions of d. The first one corresponds to the situation, when on the whole interval 0 < d ⩽ ∞ no discrete level attains the lower continuum, nor does any one return back from the continuum

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

149

Fig. 12. The shaded regions correspond to those sets of (V0 , a), when by varying d it is possible to provide the diving of the lowest even level into the continuum. The solid lines correspond to even, while the dashed ones — to odd ‘‘critical charges’’ of a single source (1).

Fig. 13. (a): The renormalized charge densities ρvRac ,2 (x), which correspond to different distances between the centers of

the sources s = 1, 3.6, 6.4 for the case a = 1, V0 = 2 and (b): the dependence of the renormalized vacuum energy E Rv ac ,2 on s for the same a = 1, V0 = 2 (color online).

(the parameters (V0 , a) correspond to d < 0 in the Eq. (60)). In this case the integral vacuum charge Qv ac ,2 keeps its value, while the vacuum energy E Rv ac ,2 , the jumps in which are entirely due to creation–annihilation of discrete levels at the lower threshold, is a continuous function of d. This regime for ρvRac ,2 (x) and E Rv ac ,2 (s) is shown in Fig. 13, calculated for V0 = 2, a = 1, which correspond to the lowest unpainted region in Fig. 12. The induced charge densities, presented in Fig. 13a, correspond to d = 0 , 1.3 , 2.7, or in terms of the distance between potential wells centers s = 1 , 3.6 , 6.4. Numerical integration confirms that in this case the integral induced charge Qv ac does not depend on s and vanishes, since the parameters (V0 , a) lie in the subcritical region and so varying s does not lead to appearance of new levels at the lower threshold, while the dependence of the renormalized vacuum energy E Rv ac ,2 (s), as it follows from Fig. 13b, is given by a continuous curve. The second regime for ρvRac ,2 (x) and E Rv ac ,2 in dependence on the distance between sources is realized for (V0 , a) , which lie in the shaded regions in Fig. 12. For such (V0 , a) with growing d one (or several) discrete levels emerge at the threshold of the lower continuum by passing through

150

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

Fig. 14. (a): The renormalized charge densities ρvRac ,2 (x), which correspond to different distances between the centers of the sources s = 1, 4, 7.6 for the case a = 1, V0 = 4.08 and (b): the dependence of the renormalized vacuum energy E Rv ac ,2 on s for the same a = 1, V0 = 4.08 (color online).

the corresponding dcr . During this process the vacuum charge Qv ac each time grows by +|e|, while E Rv ac ,2 acquires a specific jump upwards by +1. The direct calculation of functions ρvRac ,2 (x) and E Rv ac ,2 (s) (see Fig. 14) for the parameters V0 = 4.08 and a = 1, which reside in the first from below dashed region in Fig. 12, confirms these effects completely. In Fig. 14a the charge densities ρvRac ,2 (x), corresponding to different distances between the centers of the sources s = 1 , 4 , 7.6, are shown. Numerical check shows that for s = 1 and s = 4 the integral vacuum charge Qv ac = −|e|, while for s = 7.6 it grows up to Qv ac = 0 due to emergence for scr ≃ 4.0709 (dcr ≃ 1.5354) of one even level at the lower threshold. Simultaneously for s = scr there takes place a jump in the vacuum energy E Rv ac ,2 by +1 (see Fig. 14b). 5. Casimir forces between two short-range Coulomb sources The Casimir interaction energy E int v ac for the system of two short-range Coulomb sources (49) is determined through the relation (3) with subsequent renormalization. Indeed here the efficiency of the method (35), based on the integration of the logarithmic derivative of the Wronskian along the WK contours (Fig. 2) compared to evaluation of E Rv ac via the sum of the phase integral and discrete levels (30), (32), shows up most clearly, since it provides for E int v ac more analytic details, at least for d ≫ 1. Upon subtraction of 2E Rv ac (V0 , a/2) from (53) the general structure of E int v ac (d) takes the form 2 E int v ac (d) = Iint (d) − Sint (d) + Λint (d) V0 ,

(61)

with Iint (d) being the contribution from the integral term, Sint (d) – from the negative discrete levels, while Λint (d) V02 – from the renormalization term, respectively. It is appropriate to start with the renormalization term, the asymptotics of which for large d can be explored most simply and in the general form. By means of (43) and (57) the renormalization coefficient Λint (d) = Λ(a, d) − 2λ(a/2) can be represented as

Λint (d) = a/π − 2λ2 (a + d) − 2λ2 (d) − 4λ2 (a/2) + 4λ2 (d + a/2) − 2(a/(2π ) − 2λ2 (a/2)) = 4λ2 (d + a/2) − 2λ2 (a + d) − 2λ2 (d) . In the next step, by means of the asymptotics

λ2 (a → ∞) = a/π − 1/16

(62)

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

151

there follows that Λint (d) decreases exponentially at large distances between sources. Further, inserting the definition of λ2 (42) into (62), one finds

Λint (d) = 4λ2 (d + a/2) − 2λ2 (a + d) − 2λ2 (d) = 4I2 (d + a/2) − 2I2 (a + d) − 2I2 (d) √ √ ∫ ∞ 2 2 (1 − e−2a 1+y )2 e−4d 1+y 1 dy ⩽ 0 . =− 2π 0 (1 + y2 )2

(63)

So the contribution to the interaction energy E int v ac (d) from the renormalization term turns out to be strictly negative and exponentially decreasing for d ≫ 1. The exact form of the asymptotics of Λint (d ≫ 1) can be found from the expression (62) via triple integration of the MacDonald function asymptotics in the way, quite similar to the evaluation of the asymptotics of λ(a → ∞) in Section 3, and takes the form e−2s

Λint (d ≫ 1) = √

sinh2 a

(

− √

π

s

sinh a(13 sinh a − 8a cosh a)

+

16s3/2

( +O

1

))

s5/2

.

(64)

Now let us consider the behavior of the integral term in (61) for d ≫ 1, at first without subtracting the contribution from separated wells. Upon integration by parts it can be written as follows I(d) = −

1

π

∞

∫

dy Re [ln (Jred (d, iy))] ,

(65)

0

where the ‘‘reduced’’ Wronskian Jred (d, ϵ ) = J(d, ϵ )/J0 (ϵ )

(66)

contains in the nominator the Wronskian J(d, ϵ ) for the double-well potential (49)

√ 2 e−2a

J(d, ϵ ) =

1−ϵ 2

√

1 − ϵ2

√

[

1−ϵ 2 2 f2 (

f12 (ϵ ) − e−4d

] ϵ) ,

(67)

in which f1 (ϵ ) =

√

√

1 − ϵ 2 cos(a (V0 + ϵ )2 − 1)

√ √ (68) −(ϵ 2 − 1 + ϵ V0 ) sin(a (V0 + ϵ )2 − 1)/ (V0 + ϵ )2 − 1 , √ √ 2 2 f2 (ϵ ) = V0 sin(a (V0 + ϵ ) − 1)/ (V0 + ϵ ) − 1 , √ while in the denominator the Wronskian J0 (ϵ ) = 2 1 − ϵ 2 , corresponding to the free case V0 = 0. The behavior of the integral (65) for large d is found via the expansion of the integrand of the following form

( ln [Jred (d, iy)] = ln

√ e f12 (iy)

−2a

1+y2

)

1 + y2

( )2 √ √ ( ) f2 (iy) 2 2 − e−4d 1+y + O e−8d 1+y . f1 (iy)

(69)

After substituting the expansion (69) into the integral (65) one obtains two first leading terms in the asymptotics of I(d) for d ≫ 1 I(d) ≃ −

+

1

π

1

[ (

∞

∫

π 0 ∫ ∞

dy Re ln

[

dy Re e 0

√ e f12 (iy)

√ −4d

1+y2

(

−2a

1+y2

)]

1 + y2

f2 (iy) f1 (iy)

)2 ]

(70)

.

152

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

Since the first term in (70) does not depend on d, the main term in the asymptotics of the integral term in E int v ac (d) for d ≫ 1 takes the form Iint (d) = I(d) − I(d → ∞) = −

≃

1

π

2

dy Re ln (1 + y )

π

√

[

1+y2

−4d

dy Re e

(

f2 (iy)

f12 (iy)

)2 ]

f1 (iy)

0

1+y2

e2a

0

∞

∫

√

[ (

∞

∫

1

)] Jred (d, iy) (71)

.

For large d the integrand in (71) decreases rapidly with growing y, hence, the main contribution to the integral is provided by small y. Therefore it turns out to be efficient to rewrite the expression (71) in the form Iint (d) ≃

e−4d

dy Re e

π

√

[

∞

∫

−4d(

1+y2 −1−y2 /2)

(

f2 (iy)

)2 ]

f1 (iy)

0

2

e−2dy ,

(72)

and thereafter to expand the square brackets in the integrand in the power series in y. All the integrals, emerging this way, can be calculated analytically. The final expansion of Iint (d) for d ≫ 1 reads

(

e−4d

Iint (d) = V02 √

2π d

A2 2

+

(

1

3A2

8d

8

)

(

+B +O

1

))

d2

,

(73)

where z0 =

( B = A3

√

V02 − 1 ,

A=

( ctg(az0 ) −3V02 1 − + z0

+

2

,

a sin2 (az0 )

)2 A−2−

(1 + z04 ) z03

ctg(az0 )

(74)

)

a z02

1 1 + z0 ctg(az0 )

(

sin (az0 )

1 − 2V02 (1 − az0 ctg(az0 ))

)

,

It should be specially noted that the formulae (74) work equally well both for V0 > 1 and V0 < 1. For V0 = 1 upon taking in (74) the limit z0 → 0 the expressions for A and B are replaced by A= B=−

a

2

45(1 + a)4

a 1+a

,

(45 + 135a + 255a2 + 210a3 + 68a4 + 8a5 ).

√

−4d So the asymptotics of the integral term in E int / d quite v ac (d) for d ≫ 1 turns out to be ∼ e similar to the behavior of the renormalization term (64). It should be mentioned that the expansion (73) can be used also for finite d in the case, when the each next term in the expansion (69) is much less than the previous one. At the same time, there might occur an alternative situation, as for the case a = 1, V0 = 8, considered below, when the coefficients A and B turn out to be quite large. The reason is that the zero denominator in A is nothing else, but the condition for existence of the level with ϵ0 = 0 in the single well. For a = 1, V0 = 8 the lowest level is ϵ0 ≃ 0.02085, and so by sufficiently small variation of the well parameters this level can be made strictly zero. It follows whence that in the general case the expansion given above does not hold for the case, when there exists in the well the level close to ϵ0 = 0, since in this case the expansion coefficients A and B become large. In the latter case it should be taken into account by expanding the square bracket in (72) in the power series in y that the expansion of the function f1 (iy) starts now from the linear in y term, since the first term of the series cos(az0 ) + sin(az0 )/z0 vanishes. As a result, for the case ϵ0 = 0 one

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

153

obtains Iint (d) = −

V02 − 1 (1 +

e−2d + O e−4d

(

a)V02

)

,

(75)

whence it follows that for this special case the rate of decrease of the integral term in (61) for d ≫ 1 becomes sufficiently less. It should be mentioned in addition that the multiplier before the leading exponent in (75) is strictly negative, since the zero level might appear in the single well only for V0 > 1. Now let us consider the (possible) contribution to (61) from negative discrete levels for d ≫ 1. The latter coincide with the corresponding zeros of the Wronskian J(d, ϵ ) and satisfy the equation

√

f12 (ϵ ) − e−4d

1−ϵ 2 2 f2 (

ϵ) = 0 .

(76)

For d → ∞ the Eq. (76) transforms into f1 (ϵ ) = 0, which obviously is the equation for degenerate by parity levels in the system with two infinitely separated wells, or, equivalently, for the levels of a single well. Let us consider one of the levels ϵ0 in the single well, for which f1 (ϵ0 ) = 0. In the limit d → ∞ the value ϵ0 serves as the zero approximation for corresponding even and odd levels in the double-well potential (49). To find the splitting of ϵ0 into the even and odd components for finite d ≫ 1, let us seek the solution of (76) in the form ϵ = ϵ0 + δϵ , where δϵ is a small correction to ϵ0 . Inserting this expansion into (76) and decomposing the l.h.s. in δϵ including the third order with account of f1 (ϵ0 ) = 0, one obtains a cubic equation

− A1 e

√ −4d 1−ϵ02

√ −4d 1−ϵ02

+ B1 e

δϵ + C1 δϵ 2 + D1 δϵ 3 = 0 ,

(77)

where A1 = f22 (ϵ0 ) , 2f2 (ϵ0 ) B1 = − √ 1 − ϵ02

[

C1 = f1′ (ϵ0 )

[

2dϵ0 f2 (ϵ0 ) +

]2

,

√

1 − ϵ02 f2′ (ϵ0 )

]

,

(78)

D1 = f1′ (ϵ0 )f1′′ (ϵ0 ) .

Applying further in (77) the method of successive iterations, one finds the splitting of the unperturbed level ϵ0 in the form √ −2d 1−ϵ02

δϵ = ±|K1 (a)|e

√ −4d 1−ϵ02

+ K2 (a, d)e

( +O e

) √ −6d 1−ϵ02

,

(79)

where q0 = V0 + ϵ0 , K1 (a) =

(80)

(1 − ϵ02 ) (1 − q20 ) V0 (ϵ0 + q0 + aq0

√

1 − ϵ02 )

,

K2 (a, d) =

(81) (1 − ϵ02 )3/2 (q20 − 1)2 2 V02 (ϵ0 + q0 + aq0

√

1 − ϵ02 )2

√

[ 2a2 q20 (1 − ϵ02 ) (ϵ0 q0 − 1) + (2 − ϵ02 − q20 ) (ϵ0 q0 + 1) + a 1 − ϵ02 (2ϵ0 q0 (q20 − 1) + (ϵ02 − 1) (2q20 − 1)) ] √ √ × 4dϵ0 + , 1 − ϵ02 (q20 − 1) (ϵ0 + q0 + aq0 1 − ϵ02 )

(82) whereby the upper sign in (79) corresponds to the odd level, while the lower — to the even one. Here is worth to note that for discrete levels in the single well like (1) there always holds the relation q0 > 1 (for details see [29]). So both K1,2 are always well-defined, since their denominators are strictly positive.

154

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

In the case of ϵ0 < 0 for sufficiently large d both levels ϵodd and ϵev en become also negative and so their total contribution to E int v ac (d) equals to √ −4d 1−ϵ02

ϵodd + ϵeven = 2 ϵ0 + 2 K2 (a, d) e

( ) √ −6d 1−ϵ02 +O e .

(83)

So in this case the contribution to E int v ac (d) for large d, caused by negative discrete level ϵ0 < 0 in the single well, takes the form Sint (d) = ϵodd + ϵev en − 2ϵ0

(

√ −4d 1−ϵ02

= 2 K2 (a, d) e

+O e

) √ −6d 1−ϵ02

(84)

.

At the same time, the zero level ϵ0 = 0 splits for finite d into a pair, where only the even one is negative, which gives the following term in E int v ac Sint (d) = ϵev en = −

V02 − 1 (1 +

a) V02

e−2d + O e−4d

(

)

.

(85)

It should be mentioned that the analysis performed above for the discrete √ levels contribution to the

interaction energy has the correct status only subject to condition d 1 − ϵ02 ≫ 1. It means that if the single well parameters are so that the level ϵ0 lies arbitrary close to the lower threshold, the expressions (84)–(85) will be valid only for distances, which provide the fulfillment of this condition. So the resulting behavior of E int v ac (d) for d ≫ 1 to a high degree turns out to be subject of the single well configuration. If there are only positive levels in the single well, the asymptotics of ( )

√

−4d E int / d due to the integral and renormalization terms. The strictly zero level v ac (d) should be O e

ϵ0 = 0 yields the contributions to Iint (d) and Sint (d) with twice less exponent rates (75) and (85), but in E int exactly cancel each other and so there remains the same exponential law v ac (d) (this terms ) −4d of decrease O e . In presence of negative levels in the spectrum of the single well the leading term in the asymptotics of E int v ac (d) becomes different, namely E int v ac (d) = −2 K2 (a, d) e

√ −4d 1−ϵ02

( ) √ −6d 1−ϵ02 +O e ,

(86)

where ϵ0 corresponds here to the closest to the ϵF = −1 negative √ energy level. It should be mentioned that in (86) ϵ0 can be arbitrarily close to ϵF = −1, hence 1 − ϵ02 – arbitrarily small int (but nonzero). In this case √ the exponential fall-down of E vac (d) takes place only at extremely large d

subject to condition d 1 − ϵ02 ≫ 1 and the Casimir interaction between such wells acquires some features of a long-range force. It should be mentioned, however, that the essence here is ( that under )−1/2 these conditions the exponential fall-down starts at extremely large distances d ≫ 1 − ϵ02 between sources, rather than by replacement of the exponential asymptotics by a power-like, what could happen only for a massless mediator similar to considered in [21,22]. The concrete type of interaction between the wells can be quite different subject to the single well parameters V0 and a, both in the asymptotics and for finite distances between the wells. In Fig. 15 the dependence E int v ac (d) is shown for a = 1 and V0 = 4.08, 7.4, 8, 10. As it follows from Fig. 15, for d ≫ 1 and V0 = 4.08, 8, 10 the interaction energy is positive (reflecting wells), whereas for V0 = 7.4 the energy at large distances becomes negative (attracting wells). Such behavior can be easily understood by means of the analysis presented above. Actually, for V0 = 4.08, 10 in the corresponding single well the lowest discrete level is negative (−0.9648 and −0.90811, respectively). As a result, for d ≫ 1 the behavior of E int v ac (d) is defined previously by the contribution from the discrete spectrum, which in this case has the form E int v ac (d) ≃ −2K2 (a, d)e

√ −4d 1−ϵ02

→ −4dϵ0

(1 − ϵ02 )3/2 (q20 − 1)2 V02 (ϵ0 + q0 + aq0

√

1 − ϵ02 )2

√ −4d 1−ϵ02

e

>0,

(87)

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

155

Fig. 15. Different types of dependence of interaction energy between two wells on the distance d between them for a = 1 and (a,b): V0 = 4.08; (c,d): V0 = 7.4; (e,f): V0 = 8; (g,h): V0 = 10.

156

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157

√

since in the coefficient K2 (a, d) under the condition d 1 − ϵ02 ≫ 1 the main term in the square bracket in (82) will be 4dϵ0 . So in presence of a negative level ϵ0 < 0 in the ‘‘initial’’ single well the interaction energy becomes positive for sufficiently large distances between wells. For V0 = 7.4, 8 the negative levels in the single well are absent, and so the behavior of E int v ac (d) for d ≫ 1 is defined by the following expression 2 E int v ac (d ≫ 1) = Iint (d) + V0 Λint (d) ≃

V02

e−4d

√

2π d

[ ( 1

1

2

1 + z0 ctg(az0 )

]

)2 −e

−2a

2

sinh (a)

.

(88)

The sign of E int v ac (d ≫ 1) depends on the sign of the square bracket in (88). For V0 = 8 it is positive, since for these values V0 and a the expression 1 + z0 ctg(az0 ) is close to zero, as it was already mentioned above. For the same reason the positive asymptotics is reached for sufficiently larger d (see Fig. 15f). For V0 = 7.4 the square bracket in (88) is negative, and hence, for d ≫ 1 the wells attract each other (Fig. 15d). 6. Conclusion To summarize, we introduced the ln [Wronskian] contour integration techniques for calculating the Casimir effect and obtained the fluctuation-induced force between two short-range Coulomb sources mediated by one dimensional massive Dirac fermions. The main result is that essentially non-perturbative vacuum QED-effects, including the effects of super-criticality, are able to add some new properties to Casimir forces between such sources, which turn out to be more diverse compared to the case of scattering potentials with scalar coupling to fermions, considered in [20–22]. In particular, we have shown that the interaction energy between two sources can exceed sufficiently large negative values and simultaneously reveal some features similar to a long-range force, which could significantly influence the properties of such quasi-one-dimensional QED-systems. The most intriguing circumstance here is that their mutual interaction is governed first of all by the structure of the discrete spectrum of a single source, in dependence on which it can be tuned to give an attractive, a repulsive, or an (almost) compensated Casimir force with various rates of the exponential fall-down, quite different from the standard exp(−2 ms) law. Let us mention once more that the essence of the long-range interaction between sources, which appears whenever the single well contains a level ϵ0 close to the lower threshold, (is that under )−1/2 these conditions the exponential between sources, rather than by fall-down starts at extremely large distances d ≫ 1 − ϵ02 replacement of the exponential asymptotics by a power-like behavior, what could happen only for a massless mediator. No less interesting is the pattern of interaction observed in the case of sources of opposite charge, as well as in the δ -limit with sources of negligible width, which can also be explored in detail within the presented ln [Wronskian] contour integration approach. However, the detailed study of the Casimir interaction in such systems lies beyond the scope of this paper and will be presented separately. Acknowledgments The authors are very indebted to Prof. P.K. Silaev and Dr. O.V. Pavlovsky from MSU Department of Physics for interest and helpful discussions. This work has been supported in part by the RF Ministry of Sc. & Ed. Scientific Research Program, Projects No. 01-2014-63889, A16-116021760047-5, and by RFBR Grant No. 14-02-01261. References [1] [2] [3] [4]

J. Reinhardt, Greiner W. Rep. Prog. Phys. 40 (1977) 219. W. Greiner, B. Muller, J. Rafelski, Quantum Electrodynamics of Strong Fields, Springer-Verlag, Berlin-Heidelberg, 1985. G. Plunien, B. Muller, W. Greiner, Phys. Rep. 134 (1986) 87. R. Ruffini, G. Vereshchagin, S.-S. Xue, Phys. Rep. 487 (2010) 1.

Y. Voronina, I. Komissarov and K. Sveshnikov / Annals of Physics 404 (2019) 132–157 [5] [6] [7] [8] [9] [10] [11] [12]

[13] [14] [15]

[16] [17] [18] [19]

[20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33]

157

W. Greiner, J. Reinhardt, Quantum Electrodynamics, Springer Science, Business Media, 2012. V.M. Kuleshov, et al., Phys.-Usp. 58 (2015) 785. J. Rafelski, J. Kirsch, B. Muller, J. Reinhardt, W. Greiner, arXiv:1604.08690 [nucl-th] 2016. S.I. Godunov, B. Machet, M.I. Vysotsky, Eur. Phys. J. C 77 (2017) 782. A. Davydov, K. Sveshnikov, Yu. Voronina, Internat. J. Modern Phys. A 32 (2017) 1750054. Yu. Voronina, A. Davydov, K. Sveshnikov, Theoret. Math. Phys. 193 (2017) 1647. Yu. Voronina, A. Davydov, K. Sveshnikov, Phys. Part. Nucl. Lett. 14 (2017) 698. M.I. Katsnelson, Phys. Rev. B 74 (2006) 201401; A.V. Shytov, M.I. Katsnelson, L.S. Levitov, Phys. Rev. Lett. 96 (2007) 236801; K. Nomura, A.H. Macdonald, Phys. Rev. Lett. 98 (2007) 076602; V.N. Kotov, V.M. Pereira, B. Uchoa, Phys. Rev. B 78 (2008) 075433; V.M. Pereira, V.N. Kotov, A.H. Castro Neto, Phys. Rev. B 78 (2008) 085101; I.F. Herbut, Phys. Rev. Lett. 104 (2010) 066404; Y. Wang, et al., Science 340 (2013) 734; Y. Nishida, Phys. Rev. B 90 (2014) 165411. A. Davydov, K. Sveshnikov, Yu. Voronina, Int. J Mod. Phys. A 33 (2018) 1850004; Int. J Mod. Phys. A 33 (2018) 1850005. Yu. Voronina, K. Sveshnikov, A. Davydov, P. Grashin, Physica E 106 (2019) 298–311; Physica E 109 (2019) 209-224. R. Barbieri, Nucl. Phys. A. 161 (1971) 1–11; V.P. Krainov, JETP 37 (1973) 406; A.E. Shabad, V.V. Usov, Phys. Rev. Lett. 96 (2006) 180401; Phys. Rev. D 73 (2006) 125021; Phys. Rev. D 77 (2008) 025001. V.N. Oraevskii, A.I. Rez, V.B. Semikoz, JETP 45 (1977) 428; B.M. Karnakov, V.S. Popov, JETP 97 (2003) 890; M.I. Vysotsky, S.I. Godunov, Phys. Usp. 184 (2014) 206. E.B. Kolomeisky, J.P. Straley, M. Timmins, Phys. Rev. A 78 (2008) 022104. F. Romeo, Eur. Phys. J. B 89 (2016) 242. A. Flachi, M. Nitta, S. Takada, R. Yoshii, Phys. Rev. Lett. 119 (2017) 031601. E. Elizalde, M. Bordag, K. Kirsten, J. Phys. A: Math. Gen. 31 (1998) 1743; A. Bulgac, A. Wirzba, Phys. Rev. Lett. 87 (2001) 120404; A. Erdas, Phys. Rev. D 83 (2011) 025005; E. Elizalde, S.D. Odintsov, A.A. Saharian, Phys. Rev. D (2011) 105023; S. Bellucci, A.A. Saharian, Phys. Rev. D 79 (2009) 085019; L.P. Teo, Phys. Rev. D 91 (2015) 125030, 085034; A. Flachi, L.P. Teo, Phys. Rev. D 92 (2015) 025032; A. Recati, J.N. Fuchs, C.S. Peca, W. Zwerger, Phys. Rev. A 72 (2005) 023616; J.N. Fuchs, A. Recati, W. Zwerger, Phys. Rev. A 75 (2007) 043615; M. Napiorkowski, J. Piasecki, J. Stat. Phys. 156 (2014) 1136; P. Jakubczyk, M. Napiorkowski, T. Sek, Eur. Phys. Lett. 113 (2016) 30006. P. Sundberg, R.L. Jaffe, Ann. Phys. 309 (2004) 442. D. Zhabinskaya, J.M. Kinder, E.J. Mele, Phys. Rev. A 78 (2008) 060103(R). D. Zhabinskaya, E.J. Mele, Phys. Rev. B 80 (2009) 155405. E.H. Wichmann, N.M. Kroll, Phys. Rev. 101 (1956) 843. M. Gyulassy, Nuclear Phys. A 244 (1975) 497. L.S. Brown, R.N. Cahn, L.D. McLerran, Phys. Rev. D 12 (1975) 581. A. Neghabian, Phys. Rev. A 27 (1983) 2311. P.J. Mohr, G. Plunien, G. Soff, Phys. Rep. 293 (1998) 227. U. Fano, Phys. Rev. 124 (1962) 1866. W. Greiner, Relativistic Quantum Mechanics. Wave equations, third ed., 2000, pp. 197–205. R. Rajaraman, Solitons and Instantons, North-Holland, Amsterdam, 1982. K. Sveshnikov, Phys. Lett. B 255 (1991) 255. C. Itzykson, J.-B. Zuber, Quantum Field Theory, Courier Corporation, 2012. Yu. Voronina, I. Komissarov, K. Sveshnikov, Theor. Math. Phys. (2018) (in press).