Jonas Lampart

Mathieu Lewin

Semi-classical Dirac vacuum polarisation in a scalar field

We study vacuum polarisation effects of a Dirac field coupled to an external scalar field and derive a semi-classical expansion of the regularised vacuum energy. The leading order of this expansion is given by a classical formula due to Chin, Lee-Wick and Walecka, for which our result provides the first rigorous proof. We then discuss applications to the non-relativistic large-coupling limit of an interacting system, and to the stability of homogeneous systems.

Introduction

The three-dimensional, massive Dirac operator in a scalar field φ is given by

D φ = -i 3 j=1 α j ∂ ∂x j + β(1 + φ) ,
where α j , β ∈ C 4×4 are a representation of the anti-commuting Dirac-matrices with α 2 j = 1, β 2 = 1 and we have chosen units such that = c = m = 1. Under appropriate assumptions on φ it is a self-adjoint operator with domain H 1 (R 3 , C 4 ) ⊂ L 2 (R 3 , C 4 ). It arises, for instance, in nuclear physics, where the interaction between nucleons is solely modelled by various meson fields, including the scalar σ-meson φ. According to the picture of the Dirac-sea, a particle, such as a nucleon, should be described by a function ψ in the positive spectral subspace for D φ in combination with the filled 'sea' of negative energy states. The state of this particle is then described by its one-particle density matrix γ := χ (-∞,0] (D φ ) + |ψ ψ| , 1 where χ I denotes the characteristic function of I. In this picture, the spectral projection χ (-∞,0] (D φ ) describes a state with no particles, the vacuum. The total energy is then (formally) given by

Tr(D φ γ) = - 1 2 Tr |D φ | + ψ, D φ , ψ .
The first term is the formal expression for the energy of the vacuum since the spectrum of D φ is symmetric with respect to zero. Of course, -(1/2) Tr |D φ | is not finite and must be regularised.

Similar objects arise if one considers instead of the scalar field βφ an electrostatic field φ. Using various methods of regularisation, these have been studied rigorously by many authors, especially from a variational point of view [START_REF] Chaix | From quantum electrodynamics to mean field theory: I. The Bogoliubov-Dirac-Fock formalism[END_REF][START_REF] Chaix | From quantum electrodynamics to mean field theory: II. Variational stability of the vacuum of quantum electrodynamics in the mean-field approximation[END_REF][START_REF] Bach | On the stability of the relativistic electron-positron field[END_REF][START_REF] Hainzl | Non-perturbative mass and charge renormalization in relativistic no-photon quantum electrodynamics[END_REF][START_REF] Hainzl | Existence of a stable polarized vacuum in the Bogoliubov-Dirac-Fock approximation[END_REF][START_REF] Hainzl | Self-consistent solution for the polarized vacuum in a no-photon QED model[END_REF][START_REF] Hainzl | The mean-field approximation in quantum electrodynamics: the no-photon case[END_REF][START_REF] Hainzl | Existence of atoms and molecules in the mean-field approximation of no-photon quantum electrodynamics[END_REF]. The charge density of the vacuum was derived in a perturbative regime by Hainzl [START_REF] Hainzl | On the vacuum polarization density caused by an external field[END_REF].

Of the several possible regularisation procedures we will follow [START_REF] Chin | A relativistic many-body theory of high density matter[END_REF][START_REF] Serot | The relativistic nuclear manybody problem[END_REF][START_REF] Walecka | Theoretical Nuclear and Subnuclear Physics[END_REF] and define the regularised vacuum energy by

E vac (φ) = Tr -1 2 |D φ | -R(φ) , (1)
where R(φ) is the fourth order Taylor expansion of -1 2 |D φ | at φ = 0. It will be proved in Proposition 2.1 below that - 1 2 |D φ | -R(φ) is trace-class under the sole assumption that φ ∈ H 1 (R 3 , R). This is a reasonable assumption as it it is equivalent to the finiteness of the energy of the σ-field

1 2 R 3 |∇φ(x)| 2 + M 2 |φ(x)| 2 dx .
If we consider the total energy, the sum of the energies of the σ-field and the Dirac particles, the stationary points solve a system of coupled equations. We discuss applications of our results, obtained for fixed φ, to this problem in Sections 4.1 and 4.2. The dynamics of such coupled systems were studied by Sabin [START_REF] Sabin | Global well-posedness for a nonlinear wave equation coupled to the Dirac sea[END_REF] with a scalar field solving a non-linear wave equation, and in [START_REF] Hainzl | Existence of global-in-time solutions to a generalized Dirac-Fock type evolution equation[END_REF][START_REF] Sabin | Global well-posedness for a nonlinear wave equation coupled to the Dirac sea[END_REF] for an electric field. The formula (1) is very complicated to evaluate, even numerically, for a given φ. In applications it is therefore important to obtain simple approximations of E vac in different regimes. For the study of nuclear matter, Walecka [START_REF] Serot | The relativistic nuclear manybody problem[END_REF][START_REF] Walecka | Theoretical Nuclear and Subnuclear Physics[END_REF], Lee-Wick [21, p. 2306] and Chin [START_REF] Chin | A relativistic many-body theory of high density matter[END_REF]Eq. (3.53)] calculated the vacuum energy per unit volume for a constant φ ∈ (-1, +∞) and obtained the simple formula

V(φ) = - 1 4π 2 (1 + φ) 4 log(1 + φ) -P (φ) , (2)
where

P (φ) = - 1 4π 2 φ + 7 2 φ 2 + 13 3 φ 3 + 25 12 φ 4
2 is the fourth-order Taylor polynomial of the first term. In this paper, we give the first rigorous derivation of the effective vacuum energy (2), in a semi-classical regime where the field φ varies slowly. Our main result is the following Theorem 1.1 (Derivation of the effective vacuum energy). Let φ ∈ C ∞ (R 3 , R) be such that

• |φ| < 1; • (1 + |x| 2 ) (1+|α|)/2 ∂ α x φ ∈ L ∞ (R 3 ) for every multiindex α ∈ N 3 . Then we have lim ε→0 ε 3 E vac φ(ε•) = R 3 V(φ(x)) dx , ( 3
)
where V is given by (2).

The convergence (3) holds under weaker assumptions on φ, but in Section 3 we actually prove a more precise result. Namely, we show the existence of a function V ε which approximates the true vacuum energy to any order, and which coincides with V to first order. That is,

ε 3 E vac φ(ε•) = R 3 V ε (φ)(x) dx + O(ε n ), ∀n ≥ 1
where

V ε = V + O(ε 2
). The function V ε is a series in ε, each term depending on a finite number of derivatives of φ. We thereby justify as well the higher-order expansions found in later works [START_REF] Blunden | Evaluation of Dirac sea effects in a finite system[END_REF][START_REF] Chan | Derivative expansion for the one-loop effective actions with internal symmetry[END_REF][START_REF] Cheyette | Derivative expansion of the effective action[END_REF][START_REF] Perry | Effects of the Dirac sea on finite nuclei[END_REF]. Explicit formulas of the type (2) can be obtained for all odd dimensions. We will focus on the most interesting case of three-dimensional space, the corresponding proofs for other dimensions can easily be obtained by adapting our arguments line-by-line.

As we said, the regularisation used in (1) is not the only possible choice. Following [START_REF] Gravejat | Construction of the Pauli-Villars-regulated Dirac vacuum in electromagnetic fields[END_REF], one can also consider the Pauli-Villars energy

E PV (φ) := - 1 2 Tr   4 j=0 c j |D m j ,φ |  
which corresponds to introducing new masses m j and coefficients c j such that c 0 = m 0 = 1 and

4 j=0 c j (m j ) k = 0, ∀k = 0, ..., 3.
Our method can be generalized to this case, leading to the convergence

lim ε→0 ε 3 E PV φ(ε•) = - 1 4π 2 R 3   4 j=0 c j (m j + φ) 4 log(1 + φ/m j )   .
This model still depends on the regularisation masses. In the limit where m j → ∞ for all j = 1, .., 3, one exactly recovers R 3 V(φ).

The advantage of the Pauli-Villars scheme is that it preserves gauge invariance. Since the σ-model studied in this paper has no such invariance, the simpler regularisation procedure (1) works just as well. Using the Pauli-Villars method will be important for the study of the Dirac vacuum in a slowly varying magnetic field in the forthcoming work [START_REF] Gravejat | Derivation of the Heisenberg-Euler Lagrangian[END_REF]. In this case, the effective energy obtained in the limit is the famous Heisenberg-Euler Lagrangian [START_REF] Heisenberg | Folgerungen aus der Diracschen Theorie des Positrons[END_REF][START_REF] Weisskopf | Über die Elektrodynamik des Vakuums auf Grund der Quantentheorie des Elektrons[END_REF][START_REF] Schwinger | On gauge invariance and vacuum polarization[END_REF].

In physical quantities the parameter ε may correspond to /mc, where m is the mass of the Dirac particles, in which case the limit ε → 0 would be relevant to both the semi-classical limit → 0 and the non-relativistic limit c → ∞. For an interacting system, for which bound states exhibiting a certain scaling behaviour may exist, other identifications are possible. As an application, we discuss the simultaneous non-relativistic and large-coupling limit of the Dirac-sigma model in Section 4.1, where ε is /mc. Finally, in Section 4.2 we present some numerical results about the stability of a homogeneous system of Dirac particles interacting with a scalar field, a question related to nuclear matter.

The regularised vacuum energy

From now on we always work with the Dirac operator D φ including a semi-classical parameter ε

D φ := -iε 3 j=1 α j ∂ ∂x j + β(1 + φ) ,
which is of course equivalent to taking a slowly varying field φ ε (x) = φ(εx), by a unitary dilation. In this section we introduce the regularised vacuum energy and derive its basic properties. We suppose φ ∈ H 1 (R 3 , R), which ensures that βφ is infinitesimally D 0 -bounded and thus that this operator is self-adjoint with domain Dom(D φ ) = H 1 (R 3 , C 4 ). We define the regularised vacuum energy by the formula:

E vac (φ) := 1 2 Tr -|D φ | + 4 j=0 1 j! d j ds j s=0 |D sφ | . ( 4
)
Note that for small s, zero is not in the spectrum of D sφ so we can expect that |D sφ | is differentiable. Our proof that ( 4) is a correct definition relies on the integral formula

|D φ | = 1 π R D 2 φ D 2 φ + ω 2 dω = 1 2π R 2 - iω D φ + iω + iω D φ -iω dω , ( 5
)
where the integral (in the sense of Bochner) exists in the norm topology of L (H 2 , L 2 ), even if zero is in the spectrum of D φ , because in that norm the integrand decays like ω -2 for ω → ∞ and stays bounded at ω = 0.

Proposition 2.1. The regularised vacuum energy defined by equation ( 4)

is a con- tinuous function of φ ∈ H 1 (R 3 , R). It equals E vac (φ) = - 1 2π Tr R D 2 φ D 2 φ + ω 2 - 4 j=0 1 j! d j ds j s=0 D 2 sφ D 2 sφ + ω 2 dω
and is smooth on the open set

U := {φ ∈ H 1 (R 3 , R) : 0 / ∈ σ(D φ )}.
Before proving Proposition 2.1 let us prove two preparatory lemmas on the resolvent of D φ . Throughout, we will denote by S p (H ), 1 ≤ p ≤ ∞, the p-th Schatten-ideal in the bounded operators L (H ) on the Hilbert space H . The most important cases are the trace class (p = 1) and the compact operators (p = ∞).

Lemma 2.2. For every

ω ∈ R and 3 < q ≤ ∞ the map φ → βφ(D 0 + iω) -1 is continuous from H 1 (R 3 , R) to S q (L 2 (R 3 , C 4 )).
Proof. Using the Kato-Seiler-Simon (KSS)-inequality [28, Theorem 4.1]

f (x)g(-i∇ x ) Sq(L 2 (R 3 )) ≤ (2π) -3/q f L q (R 3 ) g L q (R 3 ) , (6)
and the Sobolev embedding H 1 (R 3 ) → L 6 (R 3 ) we have for 3 < q ≤ 6:

βφ(D 0 + iω) -1 Sq ≤ C(1 + ω 2 ) (3-q)/2q ε -3/q φ H 1 (p 2 + 1) -1/2 L q (7)
with some constant C > 0. Since S r → S q continuously for r > q this proves the claim.

Lemma 2.3. The set

U := {φ ∈ H 1 (R 3 , R) : 0 / ∈ σ(D φ )} is an open subset of H 1 (R 3 , R) and for every ω ∈ R the map φ → (D φ +iω) -1 from U to L (L 2 (R 3 , C 4 )
has infinitely many continuous derivatives.

Proof. If φ ∈ U and ξ ∈ H 1 (R 3 ) has sufficiently small norm, then βξD -1 φ has operator norm less than one. Hence D φ + βξ is invertible by a Neumann series and U is open. Then, the resolvent formula

(D φ + βξ ± iω) -1 -(D φ ± iω) -1 = -(D φ + βξ ± iω) -1 βξ(D φ ± iω) -1 (8)
holds and implies continuity. Letting ξ = sζ this allows us to calculate the derivative

Diff((D φ ± iω) -1 ), ζ = lim s→0 1 s (D φ + βsζ ± iω) -1 -(D φ ± iω) -1 = -(D φ ± iω) -1 βζ(D φ ± iω) -1 ,
which is obviously also continuous. Using [START_REF] Chin | A relativistic many-body theory of high density matter[END_REF] again shows that convergence to the derivative holds also in the operator norm, and the fact that the derivative is a product of differentiable functions gives smoothness.

Proof of Proposition 2.1. Consider |D sφ | given by equation [START_REF] Chan | Derivative expansion for the one-loop effective actions with internal symmetry[END_REF]. By Lemma 2.3 the integrand is a smooth function of s in a neighbourhood of s = 0. One easily checks that its derivatives at s = 0 are elements of L 1 (R, dω), L (H 2 , L 2 ) , so we can exchange integration and differentiation w.r.t. s, which establishes equality of the formula for E vac given in the Proposition with (4), in the sense that the expressions inside the trace are equal. In order to show that this expression is an operator of trace-class we first calculate

d j ds j s=0 (D sφ ± iω) -1 = (-1) j j! (D 0 ± iω) -1 βφ (D 0 ± iω) -1 j .
Then, using (5) and the resolvent formula (8), we can write the sum of the integrands as 1 2 (I(ω, φ) + I(-ω, φ)) with

I(ω, φ) : = -iω(D φ + iω) -1 + iω 4 j=0 1 j! d j ds j s=0 (D sφ + iω) -1 = iω(D φ + iω) -1 βφ (D 0 + iω) -1 5 . ( 9
)
Using the estimate [START_REF] Chin | Relativistic many-body studies of high density matter[END_REF] with q = 5 we find

I(ω, φ) S 1 ≤ iω(D φ + iω) -1 S∞ |(βφ (D 0 + iω) -1 | 5 S 1 ≤ C(1 + ω 2 ) -1 ε -3 φ 5 H 1 . ( 10
)
This shows that I(ω, φ) ∈ L 1 (R, S 1 ) for all φ ∈ H 1 , so the expression (4) is indeed well-defined. Continuity in φ is an immediate consequence of Lemma 2.2 and the Dominated Convergence Theorem. Since Diff(βφ(D 0 +iω) -1 )), ξ = βξ(D 0 +iω) -1 ) the derivatives of I(ω, φ) also satisfy an estimate of the form [START_REF] Dimassi | Spectral asymptotics in the semi-classical limit[END_REF], by the Hölder inequality for traces. Hence smoothness of E vac also follows from dominated convergence.

The semi-classical expansion of the vacuum energy

In this section we derive the formula (2) as the leading term in an expansion of E vac in powers of ε. The semi-classical expansion is naturally also an expansion in orders of derivatives acting on φ. As such, it reproduces the 'derivative expansions' discussed in the physics literature [START_REF] Blunden | Evaluation of Dirac sea effects in a finite system[END_REF][START_REF] Chan | Derivative expansion for the one-loop effective actions with internal symmetry[END_REF][START_REF] Cheyette | Derivative expansion of the effective action[END_REF][START_REF] Perry | Effects of the Dirac sea on finite nuclei[END_REF]. For this reason we make very strong smoothness assumptions on the field φ for its derivation. These can be somewhat relaxed if one is only interested in a specific order of the expansion.

We will also need to assume that |φ| < 1. The lower bound φ > -1 is certainly necessary due to the logarithm in formula (2), and the upper bound is needed to obtain the following series expansion:

Tr I(ω, φ) dω = Tr ∞ j=5 (-1) j+1 j βφ(D 0 + iω) -1 j dω . ( 11
)
To prove this formula, we use that βφ(D 0 + iω) -1 S∞ < 1 and write

Tr [I(ω, φ)] = iω Tr 1 + βφ(D 0 + iω) -1 -1 (D 0 + iω) -1 βφ(D 0 + iω) -1 5 = iω Tr (D 0 + iω) -1 βφ(D 0 + iω) -1 5 1 + βφ(D 0 + iω) -1 -1 = iω ∞ j=0 (-1) j Tr (D 0 + iω) -1 βφ(D 0 + iω) -1 j+5 .
Then we note that

ω d dω Tr βφ(D 0 + iω) -1 j = -iωj Tr (D 0 + iω) -1 βφ(D 0 + iω) -1 j ,
and integrate by parts, which is justified by the KSS-inequality (6) and gives [START_REF] Gravejat | Construction of the Pauli-Villars-regulated Dirac vacuum in electromagnetic fields[END_REF].

Now ∞ j=5 (-1) j+1 j x j = log(1 + x) -r(x)
, with a polynomial of order four r(x). As I(-ω) = I(ω) * we then find

E vac (φ) = - 1 4π Tr log 1 + βφ(D 0 + iω) -1 + r βφ(D 0 + iω) -1 (12)
+ Tr log 1 + (D 0 -iω) -1 βφ + r (D 0 -iω) -1 βφ dω .

Under sufficiently strong assumptions on φ we can represent this quantity as an asymptotic series in ε, using techniques of semi-classical pseudo-differential calculus.

We will mainly refer to [START_REF] Dimassi | Spectral asymptotics in the semi-classical limit[END_REF] for results on this topic. This book treats only scalarvalued symbols, as opposed to the C 4×4 valued symbols we use, so we will indicate in which sense the results there translate to our setting, with only minor modifications of the proofs.

Theorem 3.1. Let φ ∈ C ∞ be such that for every multiindex α ∈ N 3 there exists a constant c α , with c 0 < 1, satisfying

∂ |α| ∂x α φ(x) ≤ c α (1 + |x| 2 ) -(1+|α|)/2 .
Then E vac (φ) admits an asymptotic expansion at ε = 0 of the form

ε 3 E vac (φ) = R 3 n k=0 ε k V k (φ)(x) dx + O(ε n+1 ) ,
where V n (φ)(x) depends on the derivatives of order at most n of φ at x and

V 0 (φ)(x) = - 1 4π 2 1 + φ(x) 4 log(1 + φ(x)) -P (φ(x)) V 1 (φ)(x) = 0 .
Proof. In order to keep track of the dependence on ω we will use a pseudo-differential calculus depending on ω as an auxiliary parameter. That is, we define the symbol class S m to be those

a ∈ C ∞ (R 3 × R 3 × R, C 4×4 ) satisfying ∂ α ∂x α ∂ β ∂p β a(x, p, ω) C 4×4 ≤ C α,β 1 + |x| 2 m-|α| 2 1 + |p| 2 + ω 2 m-|β| 2
, for all multiindices α, β ∈ N 3 and appropriate constants C α,β . An ε-admissible

symbol of order m ∈ Z is then a function ε → a(ε) in C ∞ ([0, ε 0 ), S m ) with a N := d N
dε N ε=0 a ∈ S m-N , for which the remainder term of the Taylor expansion of order N at ε = 0 is of order

ε N +1 in the Fréchét topology of S m-N -1 . Let A = βφ(x)(D 0 + iω) -1 . Because (D 0 + iω) -1 is the quantisation of (α • p + β + iω) -1 there is an ε-admissible symbol a(ε) of order -1 such that A equals the ε-Weyl-quantisation Op W ε (a), i.e. (Af )(x) = 1 (2πε) 3 R 6 e iε -1 x-y,p a (x + y)/2, p, ω, ε f (y) dydp ,
for every Schwartz-function f (as follows from inspection of [START_REF] Dimassi | Spectral asymptotics in the semi-classical limit[END_REF]Proposition 7.7]). The leading term is

a 0 = βφ(x)(α • p + β + iω) -1 ∈ S -1 . Because φ ∞ < c 0 < 1, A has norm less than c 0 . Since F (z) = log(1 + z) -r(z)
(a 0 (x, p) -z) -1 C 4×4 ≤ 2 1 -c 0 , so (a 0 (x, p)-z) -1 ∈ S 0 , with bounds independent of z. This implies that (A-z) -1 = Op W ε (b(z))
, for an admissible symbol b(z) of order zero, with b 0 (z) = (a 0 (x, p)-z) -1 (cf. [10, (8.11)], note that the correct generalisation of a lower bound on a scalar symbol is an upper bound on the inverse). Thus

F (A) = 1 2πi Op W ε γ F (z)b(z) dz ,
so F (A) itself is the quantisation of an admissible symbol f of order zero. Its expansion is given by

f n = 1 2πi γ F (z)b n (z) dz ,
in particular f 0 = F (a 0 ) (the formulas for the higher order terms will differ from those for scalar-valued symbols). We now use this to evaluate the trace and the integral over ω, order-by-order in n. For n = 0 we have

f 0 = F (a 0 ) = a 5 0 ∞ j=0 (-1) j j + 5 a j 0 ∈ S -5 ,
since the sum is absolutely convergent. This implies that f 0 and its derivatives are elements of L 1 (R 6 , C 4×4 ), so (cf. [10, Theorem 9.4]) Op W ε (f 0 ) ∈ S 1 and

Tr[Op W ε (f 0 )] = 1 (2πε) 3 R 6 Tr C 4 [f 0 (x, p)] dpdx .
The leading order in the expansion of V is thus given by

V 0 (φ) = - 1 4π(2π) 3 R 3 R 3 Tr C 4 F (a 0 (x, p))] + F (a 0 (x, p) * ) dpdω = - 1 4π(2π) 3 Tr C 4 log(1 + a 0 )(1 + a * 0 ) -r(a 0 ) -r(a * 0 ) dpdω . ( 13
)
Now,

(1 + a 0 )(1 + a * 0 ) = 1 + a 0 + a * 0 + a 0 a * 0 = 1 + 2φ(x) + φ(x) 2 (p 2 + 1 + ω 2 ) -1 ,
so the C 4 -trace is trivial and expanding into power series and setting ξ = (p, ω) ∈ R 4 we obtain:

V 0 (φ) = - 1 8π 4 R 4 ∞ j=3 (-1) j+1 j (ξ 2 + 1) -j (2φ(x) + φ 2 (x)) j -r (φ(x), ξ) dξ .
The new remainder term r is a polynomial of fourth order in φ(x), and it is integrable in ξ, since the first term and the total integrand both are. The ξ-integral of the first term evaluates to

R 4 (ξ 2 + 1) -j dξ = π 2 (j -1)(j -2) ,
and the sum then yields

∞ j=3 (-1) j+1 j(j -1)(j -2) (2φ(x) + φ 2 (x)) j = (1 + 2φ(x) + φ 2 (x)) log(1 + 2φ(x) + φ 2 (x)) -r (φ(x)) = 2(1 + φ(x)) 4 log(1 + φ(x)) -r (φ(x)) ,
with another fourth-order polynomial r . With this we obtain

V 0 (φ) = - 1 4π 2 (1 + φ(x)) 4 log(1 + φ(x)) + r (φ(x)) + R 4 r (φ(x), ξ) dξ =:-P (φ(x))
.

By construction, P is a fourth order polynomial and since the first four derivatives of the integrand in equation ( 13) w.r.t. φ vanish, P (t) is the Taylor polynomial of

-1 4π 2 (1 + t) 4 log(1 + t).
We have thus found the claimed leading order term for E vac (φ) and it remains to estimate the error. First, note that f -4 n=0 ε n f n ∈ S -5 and that this implies a bound in L 1 (R, S 1 ) on the quantisation, because

Op W ε (1 + x 2 ) -5/2 (1 + p 2 + ω 2 ) -5/2 ∈ L 1 (R, S 1 ) ,
and for every σ ∈ S -5

Op W ε (1 + x 2 ) -5/2 (1 + p 2 + ω 2 ) -5/2 Op W ε (σ) ∞ ≤ C ,
by the Calderon-Vaillancourt Theorem [START_REF] Dimassi | Spectral asymptotics in the semi-classical limit[END_REF]Theorem 7.11]. This proves that the expansion of f gives an expansion of E vac (φ) of the claimed form, with

V n (φ)(x) := R R 3 Tr C 4 [f n (x, p, ω)] dpdω ,
for the terms of order at least five. It thus remains to estimate f n for 1 ≤ n ≤ 4. For brevity, we will only perform the proof for f 1 , as the arguments for higher order terms are essentially the same. The explicit form of b 1 , and thereby f 1 , is obtained from the formula for the Weyl-Moyal product (see [START_REF] Dimassi | Spectral asymptotics in the semi-classical limit[END_REF]Theorem 7.3]), which represents the product of operators on the space of admissible symbols. Denoting this product by • and by {•, •} the Poisson-bracket we have

b 1 (z) = -b 0 (z)a 1 b 0 (z) -(b 0 (z) • (a 0 -z)) 1 b 0 (z) = -(a 0 -z) -1 a 1 (a 0 -z) -1 -i 2 {(a 0 -z) -1 , a 0 }(a 0 -z) -1 . Now for every σ ∈ S m we have [σ, (a 0 -z) -1 ] = -(a 0 -z) -1 [σ, a 0 ](a 0 -z) -1 ∈ S m-1 ,
so the first term of b 1 , for example, can be written as

b 0 (z)a 1 b 0 (z) = a 1 b 0 (z) -2 + [a 0 , a 1 ]b 0 (z) -3 + [a 0 , [a 0 , a 1 ]]b 0 (z) -4 + S -5 .
On the other hand, we have

1 2πi γ F (z)σ(a 0 -z) -dz = σ 1 2πi γ F (z) 1 ( -1)! d -1 dz -1 (a 0 -z) -1 dz = σ (-1) -1 ( -1)! d -1 dz -1 F (a 0 ) ,
and since d dz F (a 0 ) ∈ S -5+l for ≤ 5 this implies that

1 2πi γ F (z)(a 0 -z) -1 a 1 (a 0 -z) -1 dz ∈ S -5 .
The same argument, applied to the second term of b 1 , gives f 1 ∈ S -5 , which proves the desired bound.

To complete the proof, we still need to check that V 1 = 0. We have

Tr C 4 1 2πi γ F (z)(a 0 -z) -1 a 1 (a 0 -z) -1 dz = Tr C 4 [a 1 F (a 0 )]
and

a 1 (x, p) = βφ(x) • (αp + β + iω) -1 1 = -i 2 {βφ, (αp + β + iω) -1 } = -i 2 ∇ x ∇ p a 0 .
Thus, for every x with φ(x) = 0,

Tr C 4 [a 1 F (a 0 )](x, p) = -i 2 ∇ p Tr C 4 [(∇ log φ(x))F (a 0 (x, p))] ,
the dp-integral of which vanishes. For the second term of b 1 , the C 4 -trace equals zero, because ∇ x a 0 commutes with a 0 and thus

Tr C 4 (∇ x (a 0 -z) -1 )(∇ p a 0 )(a 0 -z) -1 = = -Tr C 4 (a 0 -z) -1 (∇ x a 0 )(a 0 -z) -1 (∇ p a 0 )(a 0 -z) -1 = Tr C 4 (a 0 -z) -1 (∇ x a 0 )(∇ p (a 0 -z) -1 ) = Tr C 4 (∇ p (a 0 -z) -1 )(∇ x a 0 )(a 0 -z) -1 ) .
This concludes the proof of Theorem 3.1

Applications

An interacting Dirac Klein-Gordon model

In this section we discuss the implications of the results of Section 3 for a simple model of a particle, described by a Dirac equation, interacting with a classical field, described by a Klein-Gordon equation. Such models arise in the mean-field approximation of models in nuclear physics. They are often treated with the additional 'no sea' approximation, which amounts to neglecting the contribution of negative energy states. In such models, there are usually additional vector-fields to consider. However, due to gauge-invariance of the Dirac equation, the leading order of the vacuum energy is solely due to the scalar field. With the 'no sea' approximation, the static model amounts to finding critical points of the energy

E 0 (ψ, σ) := R 3 ψ(x) -icα • ∇ + β(c 2 + gσ(x)) ψ(x) + 1 2 c 2 |∇σ(x)| 2 + 1 2 M 2 c 4 |σ(x)| 2 dx , under the constraint that ψ ∈ L 2 (R 3 , C 4 ) be a normalised function in the positive spectral subspace of -icα • ∇ + β(c 2 + gφ(x))
, where we have introduced the physical constants, c (speed of light), M (mass of the field), and g (particle-field coupling constant). If we rescale lengths by ε = √ c, choose large couplings of the form g = g 0 ε -7/2 and set φ(x) = ε -3/2 g 0 σ(εx) this becomes

E 0 (ψ, φ) = ε -4 R 3 ψ(x) -iεα • ∇ + β(1 + ε 2 φ(x)) ψ(x) + 1 2 g -2 0 ε 2 ε 2 |∇φ(x)| 2 + M 2 |φ(x)| 2 dx .
Writing ψ = (ϕ, χ) with ϕ, χ ∈ H 1 (R 3 , C 2 ), the Euler-Lagrange equations with a Lagrange multiplier 1 -ε 2 µ are given by

       -iεσ • ∇χ = -ε 2 µϕ -ε 2 φϕ -iεσ • ∇ϕ = (2 -ε 2 µ + ε 2 φϕ)χ -ε 2 ∆φ + M 2 φ = -g 2 0 (ϕ 2 -χ 2 ) .
For ε small enough, these equations are known [START_REF] Lewin | Uniqueness and non-degeneracy for a nuclear nonlinear Schrödinger equation[END_REF]Theorem 3] to have a solution

(ϕ ε , χ ε , φ ε ), for which ϕ ε → vψ NLS in H 2 (R 3 , C 2 )
, where v ∈ S 3 ⊂ C 2 is an arbitrary vector and ψ NLS is the unique positive solution to the non-linear Schrödinger equation

-1 2 ∆ψ -g 0 M 2 |ψ| 2 ψ + µψ = 0 . ( 14
)
By a scaling argument, there exists a unique µ for which ψ NLS L 2 = 1, and the corresponding solution (ϕ ε , χ ε , φ ε ) is a critical point of E 0 satisfying the constraints.

If we apply the same procedure to the energy E including the contribution of the 'sea' of negative energy states, defined as in Section 2, we obtain

E(ψ, φ) = E 0 (ψ, φ) + ε -4 E vac (ε 2 φ) .
Remember that we have shown in Proposition 2.1 that E vac (φ) is smooth in an H 1neighbourhood of φ = 0, and that its first four derivatives at this point vanish.

Hence the contribution of E vac (ε 2 φ) will be small, of order ε 3 , for ε → 0, despite the pre-factor of ε -4 . It is thus easy to show, by the same implicit-function argument used in [START_REF] Lewin | Uniqueness and non-degeneracy for a nuclear nonlinear Schrödinger equation[END_REF], that this energy also has a critical point (ψ ε , φ ε ), exhibiting the same limit, as ε → 0, as without the vacuum contribution. For the term E vac (ε 2 φ) both an expansion around φ = 0 and the semi-classical expansion of Section 3 are appropriate. The second is justified by the fact that the solution to the NLS equation ( 14) is smooth and exponentially decreasing (see [START_REF] Strauss | Existence of solitary waves in higher dimensions[END_REF]), so the hypothesis of Theorem 3.1 are fulfilled in this limit.

The stability of homogeneous systems

In this section we discuss the effects of vacuum polarisation on the stability of homogeneous systems. For uniform nuclear matter, stability was already studied in the early works [START_REF] Chin | An equation of state for nuclear and higherdensity matter based on relativistic mean-field theory[END_REF][START_REF] Chin | Relativistic many-body studies of high density matter[END_REF][START_REF] Chin | A relativistic many-body theory of high density matter[END_REF] and many later contributions, but with an emphasis on the dominant effect of the vector mesons and neglecting vacuum polarisation.

In this section we present formal derivations and numerics, we will indicate how some of these can be made rigorous. We take the state of the system to be of uniform density and chemical potential µ > 0,

γ µ = χ (-∞,µ] (D φ ) .
The energy of the positive energy states is then formally given by

E 0 (µ, φ) = 1 2g 2 R 3 M 2 |φ| 2 + |∇φ| 2 dx + Tr(χ (0,µ] (D φ )) ,
where M is the mass of the field and g the coupling constant. The total energy is the sum E 0 (µ, φ) + E vac (φ). For a constant field, the energy per unit volume becomes

E(µ, φ) = M 2 2g 2 |φ| 2 + 4 (2π) 3 p 2 ≤µ-(1+φ) 2 p 2 + (1 + φ) 2 dp + V(φ) , ( 15
)
with V(φ) given by equation ( 2). We define φ vac ∈ (-1, 0] to be the local minimum of [START_REF] Hainzl | Self-consistent solution for the polarized vacuum in a no-photon QED model[END_REF] in this interval and φ 0 ∈ (-1, 0] to be that of E -V (see Figure 1). Note that φ vac is only a local minimum, as V is unbounded from below for positive values of φ. We will now study the response of the field φ vac to a local perturbation η of this system by calculating the response function T , which is formally given by

d 2 dδ 2 δ=0 E(µ, φ vac + δη) = 1 (2π) d R 3 |η(p)| 2 T (p, µ) dp ,
and compare this to the corresponding results for φ 0 . This expression is only formal since φ vac + δη / ∈ H 1 and E vac (φ vac + δη) is not defined. In order to give a proper definition, one should consider a system in a finite box and take the thermodynamic limit of the difference E(µ, φ vac + δη) -E(µ, φ vac ), with η ∈ H 1 . We will not perform this construction here, the function T will give a finite expression on the right hand side for sufficiently regular η.

The response of the vacuum can be obtained by derivation of the formula (cf. ( 12))

E vac (φ+δη) = - 1 2π R Tr log 1+δβη(D 0 +βφ+iω) -1 -r(φ+δη, ω) dω . ( 16
)
Calculations similar to those of [START_REF] Gravejat | Construction of the Pauli-Villars-regulated Dirac vacuum in electromagnetic fields[END_REF]Lemma 4.2], the details of which are given in the Appendix, give the result

T vac (p) =(2π) 3 V (φ vac ) + 2πp 2 1 0 u(1 -u) log 1 + u(1 -u)p 2 (1 + φ vac ) 2 + u(1 -u)p 2 du + 2πp 2 1 0 φ vac u(1 -u) (1 + u(1 -u)p 2 ) 2 2 -6u + 2up 2 -8u 2 p 2 + 6u 3 p 2 + φ vac (1 + 2up 2 -7u 2 p 2 + 10u 3 p 2 -5u 4 p 2 ) du .
The complicated form of this expression is largely due to the greater number of regularising terms as compared to [START_REF] Gravejat | Construction of the Pauli-Villars-regulated Dirac vacuum in electromagnetic fields[END_REF]. The response of the positive energy states is given by

-2πiT + (µ, p) = cµ R 3 Tr C 4 β(αq + β(1 + φ vac ) + z) -1 β(α(q -p) + β(1 + φ vac ) + z) -1 dqdz ,
where c µ is a path in the resolvent set of the matrices αp + β(1 + φ vac ) enclosing the interval (0, µ]. The total response function is then given by

T (µ, p) = M 2 2g 2 + p 2 2g 2 + T + (µ, p) + T vac (p) .
Numerical evaluation of these quantities, introducing the Fermi momentum k F := µ 2 -(1 + φ vac ) 2 and choosing the parameters g = 10 and M = 1 2 (this choice is suggested by [START_REF] Reinhard | The relativistic mean-field description of nuclei and nuclear dynamics[END_REF]Table 3]), gives the plots shown in Figure 2. We observe that the vacuum polarisation has a twofold stabilising effect: First, it stops φ from reaching the point -1 where the gap in the spectrum of D φ closes, and second it renders the minimizer φ more stable to local perturbations compared to φ 0 , which can be unstable at high Fermi momentum. It is an interesting open problem to make these observations rigorous. Tr r(φ + δη, ω)

= Tr -βη(D + iω) -1 2 + 2φβ(D + iω) -1 βη(D + iω) -1 2 -φ 2 β(D + iω) -1 βη(D + iω) -1 2 + 2 β(D + iω) -1 2 βη(D + iω) -1 2 ,
which is only formal as the terms are not individually trace-class, but can be made rigorous by inserting ±r(δη, ω) which is the remainder for a Dirac operator D 0 with mass 1 + φ. (

Now let f, g ∈ C (R 3 , C 4 × C 4 )
This gives us (writing f (q) = (αq + β + iω) -1 ) and g(q) = (αp + β(1 + φ) + iω) -1 )

T vac (p) = 1 2π
Tr C 4 βg(q)βg(q -p) -βf (q)βf (q -p) + 2φ(βf (q)) 2 βf (q -p) -φ 2 (βf (q)) 2 (βf (q -p)) 2 + 2(βf (q)) 3 βf (q -p) dqdω .

Using that g(q) = (αq + β(1 + φ) -iω)/(q 2 + (1 + φ) 2 + ω 2 ) the spinor trace of the first term can be evaluated, yielding Tr C 4 βg(q)βg(q -p) = 4(1 + φ) 2 -4ω 2 + 4qp -4q 2 (q 2 + (1 + φ) 2 + ω 2 )(q -p) 2

+ (1 + φ) 2 + ω 2 ) . ( 18
)
The first term involving f gives the same result with 'mass' one instead of 1 + φ, and for the remaining ones we have Tr C 4 βf (q) 1+j βf (q -p) 1+k = θ j,k (ω, q, p) (q 2 + 1 + ω 2 ) 1+j ((q -p) 2 + 1 + ω 2 ) 1+k , with the polynomials θ 1,0 =4 + 8pq -12q 2 -12ω 2 θ 2,0 =4 + 12pq -24q 2 -4(pq)q 2 + 4q 4 -24ω 2 + 4pqω 2 + 8q 2 ω 2 + 4ω 4 θ 1,1 =4 -4p 2 + 24pq -24q 2 + 4p 2 q 2 -8(pq)q 2 + 4q 4 -24ω 2 + 4p 2 ω 2 -8pqω 2 + 8q 2 ω 2 + 4ω q l i e -s(q 2 i -2q i (1-u)p i ) dq i ,

for ω and the three components of q. Evaluating these gives the terms j!k!π 2 s j+k-1 u j (1 -u) k Θ j,k (s, p, u)e -s(m+u(1-u)p 2 ) , with m = 1, or m = 1 + φ for the evaluation of ( 18), and

Θ 0,0 = -8 s + 4m 2 + 4u(1 -u)p 2 Θ 1.0 = -24 s + 4 -4p 2 + 16p 2 u -12p 2 u 2 Θ 2,0 = 1 s 24 s -48 + 8p 2 -24p 2 u + 24p 2 u 2 + 4 -4p 2 + 24p 2 u -24p 2 u 2 + 4p 4 u 2 -8p 4 u 3 + 4p 4 u 4 Θ 1,1 = 1 s 24 s + 12p 2 -36p 2 u + 24p 2 u 2 -48 -12p 2 + 36p 2 u -4p 4 u -24p 2 u 2 + 12p 4 u 2 -12p 4 u 3 + 4p 4 u 4 .
After gathering all of these terms (with the correct pre-factors) we split T vac (p) = T vac (0)+T sing (p)+T reg (p), where T sing contains all those summands whose integrand vanishes at p = 0 and diverges at s = 0. We then have T vac (0) = (2π) 3 V (φ) , since p = 0 corresponds to a perturbation by a constant field. Calculating the ds integrals for the other terms yields T sing (p) = 2π

  is a holomorphic function on the open unit disk, F (A) can be computed using Cauchy's formula. Let γ := {z : |z| = (1 + c 0 )/2}. Then we have for z ∈ γ

  Figure 1:φ vac (solid line) and φ 0 (dashed line) as func-tions of k F = µ 2 -(1 + φ 0 ) 2 and k F = µ 2 -(1 + φ vac ) 2 (respectively).

Figure 2 :

 2 Figure 2: Response functions T (µ, p) relative to φ 0 (dashed line, without vacuum corrections) and φ vac (solid line, including these) as functions of k F , at fixed p = 0.01 (left) and p = 1 (right).

We have d 2 dδ 2 δ=0

 2 Tr log 1 + δβη(D 0 + βφ + iω) -1 = -Tr βη(D 0 + βφ + iω) -1 )

  be sufficiently regular and decay at infinity. Then the operator ηg(D 0 + iω)ηf (D 0 + iω) has a continuous integral kernel andTr ηg(D 0 + iω)ηf (D 0 + iω) = 1 (2π) 3 |η(p)| 2 Tr C 4 g(αq + β + iω)f (α(q -p) + β + iω) dpdk .

4 .ξ 0 ∞ 0 s

 400 j ζ k e -aξ-bζ dξdζ =1 1+j+k u j (1 -u) k e -s(ua+(1-u)b) dsdu leads us to the Gaussian integrals R ω 2l e -sω 2 dω , R

= 2πp 2 1 0u( 1 - 1 0

 111 -u)p 2 e -su(1-u)p 2 e -s(1+φ) 2 -e -s s dsdu u) log 1 + u(1 -u)p 2 (1 + φ) 2 + u(1 -u)p 2 du ,andT reg (p) = 2πp 2 φu(1 -u) (1 + u(1 -u)p 2 ) 2 2 -6u + 2up 2 -8u 2 p 2 + 6u 3 p 2 + φ(1 + 2up 2 -7u 2 p 2 + 10u 3 p 2 -5u 4 p 2 ) du .

Acknowledgment. We acknowledge financial support from the European Research Council under the European Community's Seventh Framework Programme (FP7/2007-2013 Grant Agreement MNIQS 258023).

Appendix: Calculation of the vacuum response function

We start by taking the derivative of the integrand in [START_REF] Hainzl | Existence of atoms and molecules in the mean-field approximation of no-photon quantum electrodynamics[END_REF], writing φ = φ vac ,


Nassim Halli
email: nassim.halli@aselta.com

Henri-Pierre Charles

Jean-François Mehaut

Performance comparison between Java and JNI for optimal implementation of computational micro-kernels

Keywords: D.3.4 [Programming Languages]: Processors-Optimization, Run-time environments; D.2.8 [Software Engineering]: Metrics-Complexity measures, Performance measures Java, JNI, HPC, Performance, Vectorization

General purpose CPUs used in high performance computing (HPC) support a vector instruction set and an out-of-order engine dedicated to increase the instruction level parallelism. Hence, related optimizations are currently critical to improve the performance of applications requiring numerical computation. Moreover, the use of a Java run-time environment such as the HotSpot Java Virtual Machine (JVM) in high performance computing is a promising alternative. It benefits from its programming flexibility, productivity and the performance is ensured by the Just-In-Time (JIT) compiler. Though, the JIT compiler suffers from two main drawbacks. First, the JIT is a black box for developers. We have no control over the generated code nor any feedback from its optimization phases like vectorization. Secondly, the time constraint narrows down the degree of optimization compared to static compilers like GCC or LLVM. So, it is compelling to use statically compiled code since it benefits from additional optimization reducing performance bottlenecks. Java enables to call native code from dynamic libraries through the Java Native Interface (JNI). Nevertheless, JNI methods are not inlined and require an additional cost to be invoked compared to Java ones. Therefore, to benefit from better static optimization, this call overhead must be leveraged by the amount of computation performed at each JNI invocation. In this paper we tackle this problem and we propose to do this analysis for a set of micro-kernels. Our goal is to select the most efficient implementation considering the amount of computation defined by the calling context. We also investigate the impact on performance of several different optimization schemes which are vectorization, out-of-order optimization, data alignment, method inlining and the use of native memory for JNI methods.

INTRODUCTION AND MOTIVATION

The initial motivation for this study is to improve the performance of a Java application which intensively uses a small set of computational micro-kernels. Java has become an important general-purpose programming language and industry expresses great interest in using it for scientific and engineering computation [START_REF] Moreira | Java programming for high-performance numerical computing[END_REF][START_REF] Taboada | Java for high performance computing: Assessment of current research and practice[END_REF]. High performance scientific Java, however, suffers from the lack of optimized libraries and low-level control [START_REF] Moreira | Java programming for high-performance numerical computing[END_REF]. This is in part due to the machineindependent bytecode representation of Java which prevents from target specific optimization at the source code level. In a dynamically compiled environment, the optimizations are delegated to the JIT compiler which is in charge of translating bytecode into machine code at run-time, using dynamic information. JIT optimization like vectorization are essentials for performance and still remain a major concern [START_REF] Jantz | Performance potential of optimization phase selection during dynamic jit compilation[END_REF][START_REF] Kulkarni | Automatic construction of inlining heuristics using machine learning[END_REF][START_REF] Paleczny | The java hotspottm server compiler[END_REF]. They enable programmers to maintain a machineindependent code and save them from writing several versions for different target architectures. That is why in Java applications pure Java methods are often preferred to native ones. Though, the JIT compiler suffers from two main drawbacks. First, the JIT is a black box for developers. We have no control over the generated code nor any feedback from its optimization phases like vectorization. Secondly, the time constraint narrows down the degree of optimization [START_REF] Nuzman | Vapor simd: Auto-vectorize once, run everywhere[END_REF][START_REF] Jantz | Performance potential of optimization phase selection during dynamic jit compilation[END_REF] compared to static compiler like GCC or LLVM. As a result, it can be interesting to use statically compiled code to benefit from a deeper optimization for performance bottlenecks that usually represent a small part of the application. Java applications enable to call native code from dynamic libraries through the Java Native Interface (JNI). Nevertheless, JNI methods are not inlined and require an additional cost to be invoked compared to Java ones. Thus, to benefit from better static optimization, this call overhead must be leveraged by the amount of computation performed at each JNI invocation. This is what we call the flop-per-invocation Considering these aspects (summarized in Table 1) selecting the most efficient implementation between JNI and Java requires an advanced analysis that takes into account the flop-per-invocation. In this paper we tackle this problem and we propose to do this analysis for a set of micro-kernels in order to select the most efficient implementation for a given range of flop-per-invocation. We also investigate how the flop-per-invocation impacts the performance of several different optimization schemes which are vectorization, outof-order optimization, data alignment, method inlining and the use of native memory for JNI methods. The study is performed using the Java HotSpot Virtual Machine over a Sandy Bridge x86-64 architecture which supports the AVX (Advanced Vector Extensions) instruction set for vectorization. Furthermore, we use GCC for static compilation of the micro-kernels. Section 2 provides a background about performance analysis and Section 3 provides a background about the optimization considered for this study. Section 4 presents our benchmark methodology and experimental conditions. In Section 5 we expose and discuss the obtained results.

PERFORMANCE METRICS 2.1 Arithmetic intensity

We call flop-per-invocation the number of floating point operations (flop) performed by the kernel at each invocation. The performance of a kernel implementation is measured in flop/s (i.e. the number of floating point operations performed per second). We have the following equation:

F = AI × M (1)
Where F is the flop-per-invocation, M is the memory-perinvocation in byte corresponding to the amount of input and output memory processed by the kernel. AI is the arithmetic intensity in flop/byte which is the ratio of the flopper-invocation by the memory-per-invocation. The arithmetic intensity allows to locate micro-kernel bottlenecks [START_REF] Williams | Roofline: An insightful visual performance model for multicore architectures[END_REF]. If the bottleneck is the memory bandwidth then the kernel is called memory-bound. Otherwise the bottleneck is the CPU throughput and the kernel is called CPUbound. Indeed, we can write the following inequality:

Π > AI × β ( 2
)
Where Π is the CPU peak performance in flop/s and β is the peak memory bandwidth in byte/s. If Equation 2 is satisfied then the kernel is memory-bound (otherwise CPUbound). Since modern architectures use a multi-level cache hierarchy, the memory bandwidth depends on the data location over this hierarchy. This means that the memory bandwidth depends on the memory-per-invocation which impacts cache efficiency therefore increasing the memoryper-invocation leads to reducing the memory bandwidth.

Performance profile

We define the performance profile of a kernel implementation as the performance function of either the flop-perinvocation, or the memory-per-invocation (since there are linearly related). By plotting the performance profile of a kernel implementation for a wide range of memory-perinvocation, we can observe, for example, the memory bandwidth bottleneck.

Theorical impact of the invocation cost

A simple model to describe the invocation cost impact on performance when the flop-per-invocation decreases is given by the following equation:

P = PmaxF I + F (3)
Where P is the performance, F is the flop-per-invocation, Pmax is the peak performance reached by the kernel and I is a parameter which describes the invocation cost in flop.

We define the decrease-factor DF as followed:

DF = P Pmax = F I + F (4)
Figure 1 shows the decrease-factor as a function of the flopper-invocation for different values of I. We can see that for relatively small flop-per-invocation, the invocation cost has a large impact on performance. Accordingly, we must consider both the flop-at-invocation and the JNI invocation cost to select the most efficient implementation.

CODE OPTIMIZATION

We consider two kinds of optimization, The first ones are asymptotic optimizations that are significant when the flopper-invocation grows. In the second kind we have the optimizations that reduce the invocation cost of a method and which are significant for lower flop-per-invocation. They are described respectively in Section 3.1 and 3.2.

Asymptotic optimization

Vectorization and out-of-order optimization

Since microprocessor performance is limited by a maximum clock frequency the trend goes towards increasing processor throughput by handling multiple data in parallel. The principle of vectorization is to pack data into a vector register then perform one operation over all the data. In this study we consider the AVX instruction set available in the 64-bit Intel Sandy Bridge architecture. AVX enables to use 256 bits registers called YMM. It allows to perform four double precision operations at one time. Additionally to vector units, the out-of-order (o-o-o) execution can increase the level of parallelism by dispatching independent instructions over parallel execution units available in different ports. For example the Sandy Bridge has six execution ports. Ports 0, 1 and 5 are for arithmetic and logic operations and support 256 bit vector operation. Ports 2 and 3 are two identical memory-load ports while port 4 is for memory-store. Increasing the level of parallelism with o-o-o is done by breaking long dependency chains. Hence, the o-o-o engine can evenly distribute independent operations between the different execution ports. Mixing o-o-o and vectorization leads to an optimal instruction level parallelism. Table 3 shows an example of source code optimization using vectorization and o-o-o for a horizontal sum kernel.

The HotSpot Virtual Machine run-time uses two compilers, the client compiler and the server compiler. The server version [START_REF] Paleczny | The java hotspottm server compiler[END_REF], which is the focus in this study, is a highly optimizing bytecode compiler which targets long running applications and aims to improve asymptotic performance.

The server compiler supports auto-vectorization, it internally uses Super-word Level Parallelism (SLP) [START_REF] Larsen | Exploiting superword level parallelism with multimedia instruction sets[END_REF]. SLP detects groups of isomorphic instructions in basic blocks. This results in a lightweight vectorization approach suitable for Just-In-Time compilation. However, this approach cannot optimize across loop iterations and cannot vectorize loopcarried idioms like reduction. As a result, the JIT can not mix o-o-o and vectorization. Additionally the JIT compiler can be seen as a black box. Due to its portability specification Java does not provide machine-specific intrinsics like AVX intrinsics and of course, we can not inline assembly code directly in the Java source code. Additional metainformation, for example to ensure the compiler that pointer won't alias, cannot be passed to the JIT compiler. Finally, the JIT does not provide any feedback about its vectorization decisions which are invisible to programmers. Thus, to ensure that the code has been properly vectorized the generated code must be examined. Table 2 summarizes the Java vectorization drawbacks.

Data alignment

The proper alignment of packed data is required to observe vectorization benefit. In the considered case of a Sandy Bridge architecture, packets are 32-byte wide (which corresponds to 4 packed double). To ensure a benefit, packets must be aligned on 16 bytes which is the higher granularity for load and store instructions. With GCC we can use aligned AVX instructions if data packets are 32-byte aligned and pointers can be properly aligned using the posix memalign system call. In Java we cannot explicitly align primitive arrays and the JIT is in charge of aligning memory access. To do so, compilers perform loop peeling i.e. they begin to load or store packed data at the first aligned offset. However in some cases when one loop iterates over several arrays which are offset relative to each other, peeling cannot resolve alignment for all the arrays. This causes a serious performance penalty which remains unpredictable. As we want to control the alignment of a Java object and its memory layout, we can set statically a base offset to ensure aligned memory access. Figure 2 shows an example of a Java double array explicitly aligned on 32-byte by setting a base offset at 2 instead of 0.

Object header

Optimizations that reduce the invocation cost

Java method inlining

Inlining is a core optimization especially in dynamically compiled environment. It allows to achieve significant optimization by specializing a method considering its calling context. Inlining also eliminates the invocation cost. In this study, we only consider inlining as an optimization that eliminates the invocation cost. To do so we consider kernels which show a constant complexity through inlining. Thereby the comparison between JNI and Java is also consistent because the experiment is independent of the calling context.

JNI and native memory

JNI is a specification which allows applications to call statically compiled methods from dynamic libraries and to manipulate Java objects from there. JNI suffers from an additional invocation cost [START_REF] Grimmer | An efficient native function interface for java[END_REF][START_REF] Kurzyniec | Efficient cooperation between java and native codes jni performance benchmark[END_REF][START_REF] Stepanian | Inlining java native calls at runtime[END_REF]. Invoking a native target function requires two calls. The first is from the Java application to the JNI wrapper which sets up the native function parameters. The second is from the JNI wrapper to the native target function. But the most significant source of overhead occurs during the invocation of JNI callbacks. An example is to get an array pointer from the Java heap. A callback pays a performance penalty because two levels of indirection are used: one to obtain the appropriate function pointer inside a function table and one to invoke the function. Yet, these callbacks remains necessary to work with Java objects while avoiding memory copies that are even more expensive. In addition to the Java heap, where objects are allocated and automatically removed by the garbage collector, Java allows programmers to allocate chunks of memory out of the Java heap. This kind of memory is called native memory. We explore the use of native memory as an optimization technique in order to avoid JNI callbacks and so reduce its invocation cost.

BENCHMARK METHODOLOGY

For each kernel we provide several implementations that use the optimization techniques detailed in Section 3. For each implementation, we measure and plot its performance profile (defined in Section 2.2). We aim to select the best implementation for a given range of flop-per-invocation and also to analyze a given optimization.

Experimental Conditions

Performance results are measured on an Intel(R) Sandy Bridge Core(TM) i5-2500 3.30GHz CPU. The Linux system version is 2.6.3 and we use the OpenJDK Runtime Environment 1.8.0, 64-Bit Server VM (build 25.0-b70, mixed mode). For JNI implementation we use GCC 4.4.7. The CPU has three levels of cache, for each we give its size and approximate latency:

• The L1 data cache is 32kB, its approximate latency is around 4 cycles.

• The L2 cache is 256kB, its approximate latency is around 12 cycles

• The L3 cache is 6MB, its approximate latency is around 30 cycles

The theoretical performance peak is 41.6 Gflop/s and we measure the approximate memory bandwidth peak around 40 GB/s. According to Equation 2 in Section 2 a kernel is memory-bound (resp. CPU-bound) if its arithmetic intensity is lower (resp. greater) than 1 flop/byte.

Kernels sample and notations

We select basic micro-kernels, described in Table 5, for our experimentation. Tested kernels computes in double precision. Vectorization is always profitable since data are already packed into memory.

For each kernel we provide several implementation. An implementation is described by the following label:

T ype InvocationOpt AsymptoticOpts

Where T ype is the type of the implementation, InvocationOpt is the optimization that reduces the invocation cost and AsymptoticOpts are the asymptotic optimizations. Table 4 details the values taken by each attributes.

cpu-bound (N > 4)

Horner data-1st.

Computes a 64-degree polynomial on N data using the Horner method. Loops over the data array at first.

Measurements

Measurements are performed using timers inside a caller method. The caller executes the method we want to benchmark with a given number of iterations then returns the mean time spent executing the method. We also iterate over the caller to get the best mean time and then calculate the performance in flop/s performed by the method. Measurements are taken in a steady state considering the warm-up phase to compile the hotpots, but also the memory state. Since the caller invokes the method over the same data set, the memory bandwidth reaches its maximum for the given memory complexity. The caller arithmetic intensity is almost equal to the arithmetic intensity of the method multiplied by the method invocation count performed by the caller. Theoretically, by inlining at source code, we could increase the memory bandwidth and peak performance by swapping the loop iterating over the data with that iterating over the method computation. However, by default compilers do not perform this kind of optimization which can disturb numerical precision by switching floating point operations.

RESULTS AND ANALYSIS

For each kernel we provide the performance profile for several different implementations (Figures 3, 4, 5 and6). Performance profiles are plotted as a function of the memoryper-invocation. Thereby we can analyze the result by accounting for the cache hierarchy.

Most efficient implementation

The performance profile allows to locate the most efficient implementation for a given range of memory-per-invocation. This range is supposed to be known for a particular calling context.

The most efficient implementation of the Array addition kernel is the vectorized and inlined Java version (java inline vect) for the whole range of memory-per-invocation considered. This is an example where the JIT works as well as GCC, hence there is no benefits from using JNI in such a case. Concerning the Horizontal sum kernel, the best implementation depends on the considered range of memory-per-invocation. Since native memory breaks the Java language, we don't consider JNI with native memory as eligible. Thus, until 2kB the best implementation is the inlined Java version with an out-of-order optimization (java inline ooo). From 2kB the most efficient implementation is the vectorized and outof-order optimized JNI version (jni vect ooo). This kernel is an example where using static compilation leads to benefits but only when the amount of computation is sufficient to cover the invocation cost. The performance profile allows to quantify this threshold.

For the Horner kernel we provide two different algorithms. The coefficient-1st iterates over the polynomial coefficients at first then over the input values. The data-1st iterates over the data-values at first then over the coefficient values. The most efficient implementation considering both algorithm is provided by the JNI version which is vectorized and out-oforder optimized (jni vect ooo). This is true for the whole range of memory-per-invocation considered.

Vectorization

The efficiency of vectorization is related to the memory bandwidth. To observe a substantial speed-up, the bottleneck must be computational. For memory-bound kernels, we observe a performance decrease along with the memory-perinvocation that reduces the memory bandwidth.

Considering the Array addition kernel for a range of memoryper-invocation between L2 size and L3 size, all vectorized implementations are almost equal with a benefit of at least 30%. For memory-per-invocation greater than L3 size there is no longer a vectorization benefit since the memory bandwidth is the major bottleneck. Vectorization is really significant if the memory-per-invocation is lower than the L1 size and optimal for a memory-per-invocation of 16kB. Java implementations of the Horizontal sum kernel are not vectorized since the JIT doesn't support the vectorization of reduction idioms. Since the kernel is memory bound, vectorized JNI implementations are sensible to the memory bandwidth and from a memory-per-invocation greater than the L3 size we have no more vectorization benefits. This threshold is greater for this kernel than the Array addition one because its arithmetic intensity is also greater. Contrary to the array addition and the horizontal sum kernels, the Horner algorithm has an arithmetic intensity resolved at run-time since it depends on the polynomial degree and the number of points to evaluate. For the range of inputs considered, the kernel is CPU-bound (the polynomial degree is fixed to 64 and the number of points greater than 4). However, the coefficient-1st version is affected by the memory bandwidth like a memory-bound kernel. The coefficients fills a constant amount of memory that is relatively small and fits into the L1 data-cache (512 bytes for a 64-degree polynomial). However, for each coefficient this version needs to load all the input data which causes a lot of cache misses when the amount of data exceeds the L1 size. And this drastically reduces the memory bandwidth. The data-1st version iterates over the coefficient array through a dependency chain for each data. Data is loaded only one time and since the coefficient array fits onto the L1 datacache the memory bandwidth is optimal. However the data-1st version is not vectorized by the JIT because it exposes a reduction idiom. Vectorization of the coefficient first implementation leads to a significant speed-up factor of 2 but only for a memory-per-invocation lower than L3 size. Beyond that size, vectorization has no impact due to the memory bottleneck.

Out-of-order optimization

Out-of-order optimizations are performed by dividing a dependency chain into several independent sub dependencychains. Kernels which expose a dependency chain are the Horizontal sum and the Horner data-1st. The number of divisions to achieve an optimal out-of-order execution depends on the CPU execution units requested. For example, we need 4 dependency-chains for the Horizontal sum and 8 for the Horner data-1st. Mixing out-of-order and vectorization leads to an optimal instruction level parallelism and performance. Since the JIT can not vectorize such idioms, JNI proves to be an efficient alternative in such case. Considering the Horner data-1st kernel, mixing vectorization and out-of-order allows performance to reach the CPU peak.

Call-bound part

For all implementations, the performance penalty when the memory-per-invocation decreases is initially due to the constant cost of the loop in the measurements. In this decreas-ing part, implementations are call-bound since the invocation cost is the main bottleneck. Considering the Array addition kernel, the unvectorized Java implementation shows that this cost is rather limited since its performance remains relatively stable. The JNI and Java invocation cost is observed comparing the un-inlined version of the Java method with the JNI method which use native memory. Disabling inlining for Java methods shows that the invocation cost is limited since inlining yields a mean improvement of about 20%. Finally the comparison between native memory and heap memory for JNI methods shows that the call bound is mainly due to the invocation of JNI callbacks, first, to get a pointer from the Java array inside the heap and then to release it when the work is done. This causes a penalty of about 40%. The observation is the same for the Horizontal sum kernel where using native memory leads to a benefit of 60%.

Data Alignment

Misalignment performance penalty occurs when four-double packets are not aligned on 16-bytes. Peeling is used by compilers to begin vectorization at the first aligned packet. In Java we can't explicitly align data at allocation, moreover the garbage collector can move the data and as a result change their alignment. In Java double arrays are at least 8-byte aligned. Considering the Array addition kernel, if one array is aligned on 8-byte and the second one on 16byte then peeling is inefficient and memory access will be unaligned for one array. As showed in Figure 7, this leads to a performance decrease of around 45% which is totally unpredictable since alignment is resolved at allocation (invisible to Java programmer) and can be modified by the garbage collector.

RELATED WORKS

There is so far no performance comparison between pure Java and JNI that takes into account both the overhead of JNI calls and potential deeper optimization provided by the static compilation. Nuzman et al. presented a split autovectorization framework [START_REF] Nuzman | Vapor simd: Auto-vectorize once, run everywhere[END_REF] combining dynamic compilation with an off-line compilation stage aiming at being competitive with static compilation while conserving application portability. Parri et al. [START_REF] Parri | Returning control to the programmer: Simd intrinsics for virtual machines[END_REF] designed an API called jSIMD that uses JNI as a bridge to map Java code to SIMD instructions using vectorized data of various types. Regarding JNI performance issues, Grimmer et al. implemented the Graal Native Function Interface (GNFI) for the Graal Virtual Machine [START_REF] Grimmer | An efficient native function interface for java[END_REF] as an alternative to JNI. GNFI aims to mitigate all the disadvantages met using JNI both concerning programming flexibility than performance. Stepanian et al. [START_REF] Stepanian | Inlining java native calls at runtime[END_REF] proposed an approach for the IBM TR JIT compiler to widen the compilation span by inlining native code. Finally, Kurzinyec and Sunderam [START_REF] Kurzyniec | Efficient cooperation between java and native codes jni performance benchmark[END_REF] studied the performance of different JNI implementations for several different JVM.

CONCLUSION

In this paper we have presented a performance analysis for a set of micro-kernels considering the JIT vectorization limitation and the JNI invocation cost. By plotting the performance profile for several different implementations of a kernel, we have aimed to select the most efficient implementation for a specific amount of computation. We have showed that one major performance issue in Java concerns reduction kernels that are not vectorized. As a consequence Java implementations may suffer from severe performance penalties compared to JNI ones. By using native memory we showed that JNI suffers from a major performance penalty coming from callbacks use to access data inside the Java heap.

Figure 1 :

 1 Figure 1: Theoretic performance profile : decreasefactor function of the flop-per-invocation for different value of the invocation cost I

Figure 2 :

 2 Figure 2: Java double array 32-byte aligned. The first 32-byte aligned packet begins at the offset 2

Figure 4 :

 4 Figure 4: Horizontal sum performance profile

Figure 6 :

 6 Figure 6: Horner data-1st performance profile

  /s] j a v a j a v a _ v e c t _ u n a l i g n j a v a _ v e c t

Figure 7 :

 7 Figure 7: Array addition performance for a memoryper-invocation equal to 16kB

Table 1 :

 1 Performance drawbacks and benefits overview

		Drawbacks	Benefits
	JIT	Lower level of optimization	No call overhead
	JNI	Higher call over-head	Higher level of optimization

Table 2 :

 2 Java vectorization drawbacks summary

			HotSpot JIT	server	GCC
	Pointer Alias-ing	Possibly unresolved	Resolved with run-time checks or re-strict keyword
	Data ment	align-	Possibly unresolved	Alignment at alloca-tion
	Reduction id-iom	Unvectorized at the current state	Vectorized
	Combine Out-
	of-orders and	Not supported	Supported
	vectorization
	Auto-				Vectorization deci-
	vectorization	No feedback		sions and profitabil-
	Feedback			ity threshold
	Source vectorization code	Not supported	AVX Intrinsics or in-line assembly

Table 3 :

 3 Core vectorization and out-of-order optimizations for a horizontal sum reduction kernel in double precision

		No vectorization	Vectorization (AVX intrinsics)
	No o-o-o	f o r ( i =0; i <n;++ i ) { sum+=a [ i ] ; }	f o r ( i =0; i <n ; i +=4){ p=m m 2 5 6 l o a d u p d (&( a [ i ] ) ) ; sum=mm256 add pd ( sum , p ) ; }
			f o r ( i =0; i <n ; i +=16){
			p0=m m 2 5 6 l o a d u p d (&( a [ i ] ) ) ;
		f o r ( i =0; i <n ; i +=4){	p1=m m 2 5 6 l o a d u p d (&( a [ i + 4 ] ) ) ;
		sum0+=a [ i ] ;	p2=m m 2 5 6 l o a d u p d (&( a [ i + 8 ] ) ) ;
	o-o-o	sum1+=a [ i + 1 ] ; sum2+=a [ i + 2 ] ; sum3+=a [ i + 3 ] ; }	p3=m m 2 5 6 l o a d u p d (&( a [ i + 1 2 ] ) ) ; sum0=mm256 add pd ( sum0 , p0 ) ; sum1=mm256 add pd ( sum1 , p1 ) ; sum2=mm256 add pd ( sum2 , p2 ) ;
			sum3=mm256 add pd ( sum3 , p3 ) ;
			}

Table 4 :

 4 Implementation label description

	T ype	InvocationOpt	AsymptoticOpts
			• vect: vectorization
		• inline: inlining for	with aligned data
	• java: pure Java meth-ods + JIT	Java methods • native: native memory for JNI	• vect unalign: vec-torization with un-aligned data
	• jni: methods + JNI GCC	methods • empty if no opti-	• ooo: optimization out-of-order
		mization	• empty if no opti-
			mization

Table 5 :

 5 Micro-kernels Sample

	kernels description	arithmetic intensity	properties
	Array addition. Adds an array inside the other	1/16	memory-bound
	Horizontal sum. Sums the values of an array	1/8	memory-bound, unvectorized by the JIT
	Horner	coefficient-1st.
	Computes a 64-degree poly-
	nomial on N data using the	192N 8(64+2N +1)
	Horner method. Loops over
	the coefficients array at first.