Eduardo De Rafael

Hadronic Vacuum Polarization in QCD and its Evaluation in the Euclidean

We discuss a new technique to evaluate integrals of QCD Green's functions in the Euclidean based on their Mellin-Barnes representation. We present as a first application the evaluation of the lowest order Hadronic Vacuum Polarization (HVP) contribution to the anomalous magnetic moment of the muon 1 2 (gµ -2)HVP ≡ a HVP µ . It is shown that with a precise determination of the slope and curvature of the HVP function at the origin from lattice QCD (LQCD), one can already obtain a result for a HVP µ which may serve as a test of the determinations based on experimental measurements of the e + e -annihilation cross-section into hadrons.

I Introduction

In this and forthcoming papers we shall be concerned with QCD two-point functions of color singlet local operators with possible insertions of soft operators. These Green's functions, weighted by appropriate known functions, when integrated over the full range of their Euclidean momenta dependence, govern the hadronic contribution to many electromagnetic and weak interaction processes which appear as low energy observables. Two well known examples, with no soft insertions, are:

1. The Hadronic Vacuum Polarization two-point function:

Π µν (q) = i d 4 x e iq•x 0|T (J µ (x)J ν (0)) |0 = (q µ q ν -q 2 g µν )Π(Q 2 ) , Q 2 = -q 2 ≥ 0 , (1.1) where J µ (x) denotes the hadronic electromagnetic current in the Standard Model. For the light quarks u, d, s,

J µ (x) = (-ie) 2 3 ū(x)γ µ u(x) - 1 3 d(x)γ µ d(x) - 1 3 s(x)γ µ s(x) . (1.2)
2. The Left-Right two-point function in the chiral limit:

Π µν LR (q) = 2i d 4 x e iq•x 0|T L µ (x)R ν (0) † |0 = (q µ q ν -g µν q 2 )Π LR (Q 2 ) , Q 2 = -q 2 ≥ 0 , (1.3)
where L µ (x) = d(x)γ µ 1 2 (1 -γ 5 )u(x) and R µ (x) = d(x)γ µ 1 2 (1 + γ 5 )u(x) .

(1.4)

To lowest order in the electromagnetic coupling, the Hadronic Vacuum Polarization (HVP) contribution to the anomalous magnetic moment of the muon 1 2 (g µ -2) HVP ≡ a HVP µ is governed by a weighted integral [1,2] of the hadronic photon self-energy function Π(Q 2 ) in the Euclidean, renormalized onshell at Q 2 = 0:

a HVP µ = α π 1 0 dx(1 -x) -Π Q 2 ≡ x 2 1 -x m 2 µ . (1.5)
This integral requires knowing the function Π(Q 2 ) all the way from Q 2 = 0 (x = 0) to Q 2 = ∞ (x = 1). In the second example, the π + -π 0 mass difference in the Standard Model, to lowest order in the electroweak coupling and in the chiral limit, is governed by a weighted integral of the function Π LR (Q 2 ) in the Euclidean over the full range 0 ≤ Q 2 ≤ ∞ (see ref. [3] and references therein):

(m 2 π + -m 2 π 0 )| SM = α π 3 4f 2 π ∞ 0 dQ 2 1 - Q 2 Q 2 + M 2 Z -Q 2 Π LR (Q 2 ) , (1.6)
The first term in the r.h.s. produces a mass difference of electromagnetic origin [4] while the second one corresponds to the small contribution induced by the electroweak Z gauge boson propagator [3].

Other examples of two-point functions with soft insertions are described e.g. in ref. [5] where references to the relevant literature can also be found. The purpose of this paper is to describe a new approach to evaluate Euclidean QCD integrals of this type. The problem of computing analytically Green's functions like Π(Q 2 ) and Π LR (Q 2 ) in QCD is due to the fact that QCD perturbation theory (pQCD) is only applicable to short-distances i.e. large Q 2 values and, in fact, in the case of Green's functions which are order parameters of spontaneous chiral symmetry breaking, like Π LR (Q 2 ), pQCD gives no contribution at all in the chiral limit. The Green's functions in question, however, obey dispersion relations which relate their values in the Euclidean to integrals over spectral functions defined in the Minkowski domain. The hadronic photon self-energy Π(Q 2 ) above is not an order parameter of spontaneous chiral symmetry breaking and it obeys the subtracted dispersion relation

Π(Q 2 ) = ∞ 4m 2 π ± dt t -Q 2 t + Q 2 1 π ImΠ(t) , (1.7)
1 while Π LR (Q 2 ) obeys an unsubtracted dispersion relation

Π LR (Q 2 ) = ∞ m 2 π ± dt 1 t + Q 2 1 π ImΠ LR (t) . (1.8)
Furthermore, the two spectral functions 1 π ImΠ(t) and 1 π ImΠ LR (t) are accessible to experiment. The first one is directly accessible to experiment via the one photon e + e -annihilation cross section into hadrons (m e → 0):

σ(t) = 4π 2 α t 1 π ImΠ(t) , (1.9)
and this is in fact the way which a HVP µ has been evaluated to a high degree of precision [START_REF] Davier | [END_REF]9,[START_REF] Davier | Update of the Hadronic Vacuum Polarization Contribution to the muon g -2[END_REF]. Similarly, the hadronic spectral function 1 π ImΠ LR (t), which is the difference:

1 π ImΠ LR (t) = 1 π ImΠ VV (t) - 1 π ImΠ AA (t) (1.10)
of Vector-Vector and Axial-Axial correlation functions, is accessible via e + e -annihilation and hadronic τ -decay data. QCD perturbation theory contributes both to the VV and to the AA correlation functions, but the contribution is the same in the chiral limit and vanishes in the difference.

The two examples above are, however, rather exceptional. In most cases, the associated spectral functions to the Green's functions one is interested in are not accessible to experimental determination. In that sense, both Π(Q 2 ) and Π LR (Q 2 ) provide excellent theoretical laboratories to test non perturbative approaches to the determination of the more general type of Green's functions we are concerned with.

This first paper is dedicated to the study of the hadronic photon self-energy Π(Q 2 ) and its contribution to a HVP µ . There is a persistent discrepancy between the latest experimental determination of the anomalous magnetic moment of the muon [START_REF] Bennett | The g-2 Collab.)[END_REF]:

a exp µ = 116 592 089(63) × 10 -11 , (1.11)
and its theoretical evaluation in the Standard Model (see e.g. ref [START_REF] Th | The Muon (g -2) Theory Value: Present and Future[END_REF] for a recent description of the various contributions as well as for earlier references.) The lowest order HVP contribution to a th µ , evaluated from a combination of experimental results on e + e -data: a HVP µ = (6.923 ± 0.042) × 10 -8 [START_REF] Davier | [END_REF] and a HVP µ = (6.949 ± 0.043) × 10 -8 [9] ,

(1.12)

gives at present the contribution with the largest error. [START_REF] Burger | [END_REF][START_REF] Chakraborty | HPQCD collaboration)[END_REF][START_REF] Sz | [END_REF] and references therein), may eventually serve as a test of the above results. This, and the planned experiments at Fermilab and JPARC to measure a µ in the near future, which aim at reducing the present uncertainty in a exp µ by a factor of four, has prompted a renewed activity on this topic.

This paper is organized as follows. In the next section we discuss the Mellin-Barnes representation of the HVP-function Π(Q2 ). Section III recalls various equivalent representations which can be used to evaluate a HVP µ . Section IV is dedicated to an application of Ramanujan's Master Theorem to the HVP-function and Section V to the Marichev Interpolation of the Mellin transform of the HVP-Spectral Function, which is in fact the approach that we propose. Our procedure to evaluate a HVP µ is discussed in Section VI. First in the case where only the first Mellin moment M(0) is known, i.e. when only the slope of Π(Q 2 ) at the origin is known; and also when both the first two Mellin Moments M(0) and M(-1) are known, i.e. when the slope and curvature of Π(Q 2 ) at the origin are known. We illustrate the Marichev Interpolation Approach with the example of Vacuum Polarization in QED, and we test it with a phenomenological Toy Model [START_REF] Lellouch | Private Comunication[END_REF] which reproduces the basic features of the hadronic spectral function. We then apply the same interpolation approach to the BHLS model of ref. [21] as well as to a recent LQCD evaluation [START_REF] Sz | [END_REF] of the first two moments M(0) and M(-1), and finally we conclude.

II Mellin-Barnes Representation of the HVP-Function.

We shall extensively use the fact that the electromagnetic hadronic self-energy function Π(Q 2 ) in Eq. (1.7) has a useful representation in terms of the Mellin transform of its spectral function 1π ImΠ(t) defined as follows [START_REF] Rafael | [END_REF]:

M 1 π ImΠ(t) (s) ≡ M(s) = ∞ t0 dt t t t 0 s-1 1 π ImΠ(t) , t 0 = 4m 2 π ± , -∞ ≤ Re(s) < 1 , (2.1)
where we have normalized the spectral function t-variable to the two-pion threshold value 1 . With this normalization M(s) is dimensionless and monotonously decreasing in the range: -∞ ≤ Re(s) < 1.

In QCD, the Mellin transform M(s) is singular at s = 1 with a residue which is fixed by pQCD. The contribution from the three light u, d, s quarks gives

M(s) ∼ s→ 1 α π 4 9 + 1 9 + 1 9 N c 1 3 [1 + O(α s )] 1 1 -s . (2.
2)

The representation of Π(Q 2 ) in terms of M(s) follows from inserting the Mellin-Barnes identity

1 1 + Q 2 t = 1 2πi cs+i∞ cs-i∞ ds Q 2 t -s Γ(s)Γ(1 -s) (2.3)
in the integrand of the r.h.s. of the dispersion relation in Eq. (1.7). Then, one has

Π(Q 2 ) = -Q 2 ∞ t0 dt t 2 1 2πi cs+i∞ cs-i∞ ds Q 2 t -s Γ(s)Γ(1 -s) 1 π ImΠ(t) = - Q 2 t 0 1 2πi cs+i∞ cs-i∞ ds Q 2 t 0 -s Γ(s)Γ(1 -s) M(s) , c s ≡ Re(s) ∈]0, 1[ . (2.4)
The interest of this integral representation is encoded in the so called converse mapping theorem of ref. [22] (see also refs. [23,24,25] for applications in QED). This theorem relates the singularities in the complex s-plane of the integrand, i.e. the singularities of Γ(s)Γ(1 -s) M(s) in our case, to the asymptotic expansions of Π(Q 2 ) for Q 2 large and for Q 2 small. These relations are as follows:

• Expansion for Q 2 → ∞
In the r.h.s. of the fundamental strip defined by c s ≡ Re(s) ∈]0, 1[ in Eq. (2.4) i.e. for Re(s) ≥ 1, the most general singular expansion 2 of the function Γ(s)Γ(1 -s) M(s) is of the following type:

Γ(s)Γ(1 -s) M(s) ≍ p=1 k=0
a p,k (s -p) k+1 where here (p, k) ∈ N .

(2.5)

The corresponding asymptotic behaviour of Π(Q 2 ) for Q 2 large ordered in increasing powers of t 0 /Q 2 is then:

Π(Q 2 ) ∼ Q 2 →∞ - Q 2 t 0 p=1 k=0 (-1) k+1 k! a p,k t 0 Q 2 p log k Q 2 t 0 . (2.6)
From Eq. (2.5) there follows that the possible presence of powers of log Q 2 terms in this expansion is correlated to the possible singular behaviour of the Mellin transform M(s) at s = 1, 2, 3, • • • .

In particular, the leading behaviour of Π(Q 2 ) for Q 2 → ∞, which is controlled by pQCD, corresponds to p = 1 and k = 1, and

Π(Q 2 ) ∼ Q 2 →∞ -a 1,1 log Q 2 t 0 + • • • , (2.7)
where, from the light u, d, s quarks contribution,

a 1,1 = α π 2 3 1 3 N c [1 + O(α s )] . (2.8)
In other words, the singular behaviour of the Mellin transform in Eq. ( 2.1) at s = 1 is correlated with the leading log Q 2 behaviour of Π(Q 2 ) when Q 2 → ∞ and with the fact that the hadronic spectral function goes asymptotically to a constant:

1 π ImΠ(t) ∼ t→∞ α π 2 3 1 3 N c [1 + O(α s )] .
(2.9)

• Expansion for Q 2 → 0
In the l.h.s. of the fundamental strip i.e. for s ≤ 0, the most general form of the singular expansion of the function Γ(s)Γ(1 -s) M(s) is of the following type:

Γ(s)Γ(1 -s) M(s) ≍ p=0 k=0 b p,k (s + p) k+1 , (2.10)
and the corresponding asymptotic behaviour of Π(Q 2 ) ordered in increasing powers of Q 2 /t 0 is then:

Π(Q 2 ) ∼ Q 2 →0 - Q 2 t 0 p=0 k=0 (-1) k b p,k Q 2 t 0 p log k Q 2 t 0 . (2.11)
In QCD the Mellin transform in Eq. (2.1) for Re(s) < 1 is not singular and, therefore, there are no powers of log Q 2 in the expansion of Π(Q 2 ) for Q 2 → 0. The expansion in this region is a power series:

Π(Q 2 ) ∼ Q 2 →0 - Q 2 t 0 n=0 b n Q 2 t 0 n , (2.12)
and the coefficients b n are fixed by the moments M(s) at s = -n, with n = 0, 1, 2, 3, • • • . From the dispersion relation in Eq. (1.8) there follows that these moments correspond to successive derivatives of the HVP self-energy Π(Q 2 ) at the origin:

M(-n) = ∞ 0 dt t t 0 t 1+n 1 π ImΠ(t) = (-1) n+1 (n + 1)! (t 0 ) n+1 ∂ n+1 (∂Q 2 ) n+1 Π(Q 2 ) Q 2 =0
.

(2.13)

More precisely b n = (-1) n M(-n) . ( 2

.14)

We conclude that the determination of a few terms of the Taylor expansion of Π(Q 2 ) in LQCD, i.e. of a few derivatives of Π(Q 2 ) at Q 2 = 0, is equivalent to a determination of the Mellin transform of the physical spectral function at a few discrete values s = -n which we call the Mellin Moments.

III

The HVP Contribution to g µ -2.

There are several equivalent representations of a HVP µ which we next recall.

• The Standard Representation in terms of the Hadronic Spectral Function [26] a

HVP µ = α π ∞ t0 dt t 1 0 dx x 2 (1 -x) x 2 + t m 2 µ (1 -x) 1 π ImΠ(t) . (3.1)
This is the traditional representation for a determination of a HVP µ when using experimental data and/or phenomenological models. From this representation one can easily see that the integral over the Feynman parameter x is a function of m 2 µ t which decreases monotonously as t runs from the hadronic threshold t 0 = 4m 2 π to t = ∞. The contribution to a HVP µ is, therefore, dominated by the low-t behaviour of the spectral function. In fact, from the inequality

x 2 (1 -x) x 2 + t m 2 µ (1 -x) ≤ x 2 m 2 µ t , (3.2)
there follows a rigorous upper bound for a HVP µ [27]:

a HVP µ ≤ α π 1 3 ∞ t0 dt t m 2 µ t 1 π ImΠ(t) = α π 1 3 m 2 µ - ∂Π(Q 2 ) ∂Q 2 Q 2 =0 , (3.3)
where the equality in the r.h.s. results from the relation in Eq. (2.13) when n = 0. In other words, a HVP µ is bounded by the Slope of the HVP at the Origin i.e. by the First Mellin Moment M(0). As emphasized in ref. [START_REF] Rafael | [END_REF]: "the comparison between the determinations of M(0) from LQCD and experimental results should provide an important first test".

Using the dispersion relation in Eq. (1.7) one can rewrite the parametric representation in Eq. (3.1) as follows:

a HVP µ = α π ∞ t0 dt t 1 0 dx x 2 (1 -x) x 2 + t m 2 µ (1 -x) 1 π ImΠ(t) = α π 1 0 dx(1 -x) ∞ t0 dt t x 2 1-x m 2 µ t + x 2 1-x m 2 µ 1 π ImΠ(t) = α π 1 0 dx(1 -x) -Π x 2 1 -x m 2 µ , (3.4)
resulting in the representation for a HVP µ quoted in Eq. (1.5) i.e.

• The Representation in terms of the Euclidean Photon Self-Energy [1, 2]

a HVP µ = α π 1 0 dx(1 -x) -Π x 2 1 -x m 2 µ , Q 2 ≡ x 2 1 -x m 2 µ .
Trading the Feynman parametric x-integration by an integration over the Euclidean Q 2 -variable results in a more complicated expression

a HVPL µ = α π ∞ 0 dQ 2 Q 2 Q 2 4m 2 µ + Q 2   4m 2 µ + Q 2 -Q 2 4m 2 µ + Q 2 + Q 2   2 [-Π(Q 2 )] , (3.5)
which is the one proposed in ref. [28] for LQCD determinations of a HVPL µ . This requires, however, an interpolation procedure to evaluate Π(Q 2 ) in the Q 2 -regions where there is not a direct LQCD evaluation. It is precisely this interpolation procedure which is the main concern of this paper.

• The Representation in terms of the Adler Function

a HVP µ = α π 1 2 1 0 dx x(2 -x)A Q 2 ≡ x 2 1 -x m 2 µ (3.6)
where

A(Q 2 ) = -m 2 µ ∂Π(Q 2 ) ∂Q 2 , (3.7)
follows from the one in Eq. (1.5) integrating by parts and using the fact that Π(0) = 0.

In terms of the Euclidean Q 2 -variable:

a HVP µ = α π ∞ 0 dω G(ω) A(ωm 2 µ ) , (3.8)
where

G(ω) = 1 4 (2 + ω)(2 + ω -ω) √ 4 + ω) -2 and ω = Q 2 m 2 µ . (3.9)
In this representation the upper bound [27] in Eq. (3.3) follows from the positivity of the spectral function 1 π ImΠ(t) which gives rise to the inequality

- d dω Π( ωm 2 µ ≤ - d dω Π ωm 2 µ | ω=0 , (3.10)
and the fact that the function

G(ω) ≡ 1 4 (2 + ω) 2 + ω - √ ω √ 4 + ω -2 , (3.11)
is positive and monotonously decreasing. Therefore

a HVP µ ≤ α π - d dω Π ωm 2 µ | ω=0 ∞ 0 dω 1 4 (2 + ω) 2 + ω - √ ω √ 4 + ω -2 = α π 1 3 - d dω Π ωm 2 µ | ω=0 = α π 1 3 ∞ 4m 2 π dt t m 2 µ t 1 π ImΠ(t) . (3.12)
The function G(ω) has the following asymptotic behaviours:

G(ω) ∼ ω→0 1 2 - √ ω + ω - 5 8 ω 3/2 + 1 4 ω 2 + O[ω 5/2 ] , (3.13)
and

G(ω) ∼ ω→∞ 1 2ω 2 - 2 ω 3 + 7 ω 4 - 24 ω 5 + O[ 1 ω 11/2 ] . (3.14)
• The Mellin-Barnes integral representation [START_REF] Rafael | [END_REF] Inserting the Mellin-Barnes expression for Π(Q 2 )that we obtained in Eq. (2.4) in the Euclidean representation of Eq. (3.5), and explicitly integrating over the x-parameter, results in a Mellin-Barnes representation for a HVP

µ a HVP µ = α π m 2 µ t 0 1 2πi cs+i∞ cs-i∞ ds m 2 µ t 0 -s F (s) M(s) , (3.15)
where F (s) denotes the product of Gamma-functions

F (s) = -Γ(3 -2s)Γ(-3 + s)Γ(1 + s) . ( 3
M(s) = - d ds M(s) = ∞ t0 dt t t t 0 s-1 log t 0 t 1 π ImΠ(t) (3.17) at s = -1, -2, • • • which give negative contributions: a HVP µ = α π m 2 µ t 0 1 3 M(0) + m 2 µ t 0 25 12 -log t 0 m 2 µ M(-1) + M(-1) + m 2 µ t 0 2 97 10 -log t 0 m 2 µ M(-2) + 6 M(-2) + m 2 µ t 0 3 208 5 -log t 0 m 2 µ M(-3) + 28 M(-3) + O m 2 µ t 0 4    . (3.18)
The bulk of the overall contribution to a HVP µ comes in fact from just the first few terms. The first term is the upper-bound of ref. [27] with successive fast improvements from the following terms. This expansion, when tested with phenomenological models [START_REF] Lellouch | Private Comunication[END_REF]21], reproduces the full answer to a good accuracy with just the terms in the first two lines.

The Mellin Moments M(n) and their derivatives M(n) at n = 0, -1, -2, • • • provide tests for a comparison between experimental results, phenomenological models and combined LQCD-Padé determinations. They can also be used as an alternative way to evaluate integrals like the ones we shall encounter later in Eqs. (5.25) and (5.30).

IV Ramanujan's Master Theorem

Let us reconsider the Taylor series expansion in Eq. (2.12) which, as we have shown, is governed by the Mellin Moments M(-n), n = 0, 1, 2, • • • , i.e.

- t 0 Q 2 Π(Q 2 ) ∼ Q 2 →0 M(0) - Q 2 t 0 M(-1) + Q 2 t 0 2 M(-2) - Q 2 t 0 3 M(-3) + • • • . (4.1)
This expansion provides the basis for an application of Ramanujan's Master Theorem [START_REF] Berndt | Ramanujan's Notebooks, Part I[END_REF] to our case. It states that:

∞ 0 d Q 2 t 0 Q 2 t 0 s-1 M(0) - Q 2 t 0 M(-1) + Q 2 t 0 2 M(-2) + Q 2 t 0 3 M(-3) + • • • = Γ(s)Γ(1 -s)M(s) . (4.2)
This integral identity, which follows from the inverse transform in Eq. (2.4), is at the basis of the approach that we are going to use.

For pedagogical purposes let us first apply it to the simple case of Vacuum Polarization in QED.

IV.1 Application to Vacuum Polarization in QED

Consider the Euclidean behaviour of vacuum polarization in QED for a fermion of mass m which, to lowest order in α, is given by the simple Feynman parametric integral

Π QED (Q 2 ) = - α 2π 1 0 dy(1 -y 2 ) log 1 + Q 2 4m 2 (1 -y 2 ) . (4.3)
For Q 2 -small it has the Taylor series expansion

- 4m 2 Q 2 Π QED (Q 2 ) ∼ Q 2 →0 α 2π n=0 Q 2 4m 2 n (-1) n n + 1 1 0 dy(1 -y 2 ) 2+n , (4.4)
and since

1 0 dy(1 -y 2 ) 2+n = √ π 2 Γ(3 + n) Γ( 7 2 + n) , (4.5)
it can be expressed as follows

- 4m 2 Q 2 Π QED (Q 2 ) ∼ Q 2 →0 n=0 (-1) n Q 2 4m 2 n α 2π 1 n + 1 √ π 2 Γ(3 + n) Γ( 7 2 + n) . (4.6)
Recall that, as discussed earlier, the successive derivatives of Π QED (Q 2 ) at Q 2 → 0 are given by the Mellin Moments M QED (-n), n = 0, 1, 2, • • • , of the QED spectral function.

The application of Ramanujan's Master Theorem to this Taylor series is straightforward. First, it states that the Mellin transform of

-4m 2 Q 2 Π QED (Q 2
) is related to the full Mellin transform of the spectral function M QED (s) as follows:

∞ 0 d Q 2 4m 2 Q 2 4m 2 s-1 - 4m 2 Q 2 Π QED (Q 2 ) = Γ(s)Γ(1 -s)M QED (s) , (4.7)
and, furthermore, that the full Mellin transform function M QED (s) can be simply obtained by the replacement n → -s in the n-dependent coefficient of the previous Taylor series. By simple inspection of the Taylor series we conclude, without having to do any integral, that

M QED (s) ≡ ∞ 4m 2 dt t t 4m 2 s-1 1 π ImΠ QED (t) = α π 1 3 1 1 -s 3 √ π 4 Γ(3 -s) Γ( 7 2 -s) , (4.8)
where 1 π ImΠ QED (t) is the lowest order QED spectral function

1 π ImΠ QED (t) = α π 1 3 1 + 2m 2 t 1 - 4m 2 t θ(t -4m 2 ) . (4.9)
The function M QED (s) thus obtained is the analytic continuation to the complex s-plane of the function defined in the region s < 1 (i.e. the fundamental strip) where the Mellin Moments M QED (-n), n = 0, 1, 2, • • • are well defined by direct integration of the QED spectral function.

We shall come back to this simple example later on.

IV.2 Ramanujan's Theorem and the HVP-Function

Ramanujan's theorem applied to the HVP self-energy function in Eq. (4.2) guarantees that the incorporation of more and more moments M(-n) at integer n-values n = 0, 1, 2, • • • converges to the full Mellin function M(s). The fact that these moments are numerically accessible to LQCD (at least for the low n-values [START_REF] Sz | [END_REF]) via the determination of the derivatives of the Euclidean hadronic self-energy function Π(Q 2 ) at Q 2 = 0, provides an interesting starting point towards an alternative evaluation of the HVP contribution to a HVP µ from first principles. Ramanujan's Theorem, however, does not tell us which is the best Interpolating Function we should use to approximate the exact M(s) function when one only knows numerically a few M(-n) moments. Padé Approximants to Π(Q 2 ) [29], or the method of conformal polynomials [30], cannot be the answer because they fail to reproduce the pQCD behaviour at s = 1 of M(s). Padé Approximants to Π(Q 2 ) in the low-Q 2 -region, up to a "reasonable" Q 2 0 -value from which onwards the pQCD prediction for Π(Q 2 ) takes over (see e.g. ref. [31] and references therein), is a possible way to proceed but to our knowledge it has not been proved to be the best interpolation procedure.

These considerations have prompted us to investigate alternative approaches based on functional interpolations of Mellin Moments which respect known properties of QCD, in particular the fact that M(s) is singular at s = 1. In the next section we present a new technique in this direction inspired by Marichev's class of Mellin transforms [START_REF] Marichev | Handbook of Integral Transforms of Higher Transcendental Functions: Theory and Algorithmic Tables[END_REF] (see also ref. [START_REF]Transforms and Applications Handbook, Third Edition[END_REF], Chapter 12), shown to be applicable to a large class of functions. Although we cannot prove that this approach is the best interpolation procedure for the QCD Green's functions that we are concerned with, it turns out to be surprisingly successful when tested with the previous QED example, as well as with phenomenological models of the hadronic spectral function which we later discuss.

V Marichev's Interpolating Approach.

The most general form of a Mellin transform of Marichev's class is a fraction involving products of Gamma-functions:

M(s) = C i,j,k,l Γ(a i -s)Γ(c j + s) Γ(b k -s)Γ(d l + s) , (5.1)
where C and a i , b k , c j and d l are real constants and the variable s appears only with a ± coefficient.

In our case, these constants will be adjusted so as to reproduce as well as possible properties that we know of the QCD Mellin transform of the physical spectral function in Eq. (2.1). In particular, the choice of the ± signs in the Γ-functions of the interpolating expressions that we shall consider must result in a function monotonously decreasing in the range 1 > Re(s) ≥ -∞. This excludes interpolations which produce poles and/or zeros in this region.

The inverse Mellin transform of a function of the Marichev class is a generalized hypergeometric function which can be reconstructed using the Slater procedure [START_REF] Slater | Generalized Hypergeometric Functions[END_REF]. This way, one can obtain the underlying spectral function, as well as the underlying self-energy function in the Euclidean, corresponding to a given Marichev interpolation 3 . The fact that practically all the known functions in Mathematical Physics can be expressed as generalized hypergeometric functions gives a strong support to the interpolation approach that we are advocating.

Let us first illustrate how the Marichev interpolating approach works in the simple case of vacuum polarization in QED that we discussed before and where the same requirement of monotonously decreasing for the Mellin transform also applies:

V.1 Application to Vacuum Polarization in QED

The QED Mellin transform in Eq. (4.8) is indeed of the Marichev class. We shall show below that, in this case, the interpolation method we propose converges very fast to the exact result.

• Assume that we only know the asymptotic behaviour of the QED spectral function i.e. As already discussed, this implies that M(s) has a pole at s = 1 with α π 1 3 as the residue at the pole and, therefore, fixes what we shall call in this case the First Marichev Interpolation to a simple ratio of Gamma-functions:

M QED (s) ⇒ M (1) (s) = α π 1 3 1 1 -s = α π 1 3 Γ(1 -s) Γ(2 -s) , (5.3)
the upper-script (1) in M (1) (s) meaning that we have only used as information the value of M QED (s) at s = 1.

• The next step will use the information that, besides the singular behaviour at s = 1 we also know the slope of Π(Q 2 ) at Q 2 = 0 which, as previously discussed, is equivalent to say that we know the first Mellin Moment at s = 0. In QED:

M QED (0) = α π 1 3 4 5 , (5.4)
and with this information we can now improve our ansatz to a Second Marichev Interpolation:

M QED (s) ⇒ M (1)
(0) (s) = α π 1 3 1 1 -s Γ(b -1) Γ(b -s) (5.5)
which, at s = 1, does not change the singular behaviour of M QED (s) and it satisfies the requirement of being a function monotonously decreasing in the range 1 > Re(s) ≥ -∞. The new parameter b can then be fixed from the identity

M QED (0) (1)
(0) = M QED (0) ⇒ b = 9 4 , (5.6)
and therefore

M QED (s) (1),(2) = α π 1 3 1 1 -s Γ(5/4) Γ(9/4 -s) . (5.7)
• A further improvement results when we add the information that we also know the term of

O M 2 Q 2 in the asymptotic expansion of Π QED (Q 2 ) when Q 2 → ∞.
According to our discussion in Section II this is equivalent to say that

M QED (2) = α π 1 3 - 3 2 , (5.8)
and it allows us to consider an improved Third Marichev Interpolation function satisfying the requirement of being monotonously decreasing in the range 1 > Re(s) ≥ -∞:

M QED (s) ⇒ M QED (s) (1),(2) (0)
= α π

1 3 1 1 -s Γ(c -1) Γ(c -s) Γ(d -s) Γ(d -1) . (5.9)
The new parameters c and d are fixed by matching this ansatz to the physical values of M QED (s) at s = 0 and s = 2, which implies the equations:

d -1 c -1 = 4 5 , and c -2 d -2 = 3 2 , (5.10)
with the results c = 7 2 , and d = 3 , (

and therefore

M (1), (2) (0)
(s) = α π

1 3 1 1 -s Γ(5/2) Γ(3 -s) Γ(7/2 -s) . (5.12)
Quite remarkably, we find that this Third Marichev Interpolation already coincides with the Exact Result for M QED (s) in Eq. (4.8)! (recall that Γ(5/2) = 3 √ π 4 ). If we try to improve the Third Marichev Interpolation with further information, e.g. the knowledge of M QED (-1), we find that the new input value coincides exactly with the one predicted by M (1),(2) (0) (s) at s = -1 and, therefore, there is no room for further improvement.

We have found that the Marichev interpolation approach that we are advocating, when applied to vacuum polarization in QED, reproduces the exact expression of the Mellin transform of the spectral function with just the information provided by the values of three Mellin Moments. Encouraged by this remarkable success we propose to apply the same approach to vacuum polarization in QCD which will be discussed in the next subsection; but, before we leave this QED example, we still want to comment on another issue: the calculation of the QED vacuum polarization contribution to the anomalous magnetic moment of an external fermion. We explain this in the following sub-subsection. With s = 1 2 -iτ as a choice, and after some simplifications, the previous integral becomes:

a QED (VP) = α π 2 1 4 1 2π +∞ -∞ dτ          1 1 + τ 2 - 10 4 + τ 2 + 40 25 + τ 2 π 2 [cosh(πτ )] 2 ℜ(τ ) -iτ 1 1 + τ 2 - 5 4 + τ 2 + 16 25 + 4τ 2 π 2 [cosh(πτ )] 2 = α π 2 × 0.01568742185910 , (5.14)
which, to this remarkable accuracy, agrees with the exact analytic result (see e.g. ref. [24]):

a QED (VP) = α π 2 119 36 - π 2 3 .
(5.15)

Notice that the imaginary part of the integrand in the second line of Eq. (5.14) gives zero contribution to the integral. The shape of the function R(τ ) in the first line, which is symmetric under τ → -τ , defines the real part of the integrand and it is shown in Fig ( 1).

V.2 Application to Vacuum Polarization in QCD

We can now go back to QCD. As already stated in Eq. (2.2) the QCD Mellin transform M(s) is singular at s = 1 with a residue which is fixed by pQCD. For three light u, d, s quarks and neglecting O(α s ) corrections this fixes the pQCD-Marichev Interpolation to the simple function

M(s) ⇒ M (1) (s) = C 1 1 -s , (5.16)
where C will denote the overall constant

C ≡ α π 1 3 N c 2 3 .
(5.17)

Contrary to the QED case discussed before, the O(1/Q 2 ) term in the asymptotic expansion of Π(Q 2 ) for large Q 2 in QCD vanishes in the chiral limit which requires that

M(2) = 0 , (5.18)
and the corresponding interpolation is then

M (1,2) (s) = C 1 (1 -s)Γ(2 -s) . (5.19)
This is as much as we shall use from the short-distance behaviour of QCD.

V.2.1 First Marichev Interpolation when M(0) is known

Let us now construct the QCD equivalent of the Marichev Interpolation which in the QED example above already reproduced the exact result. This corresponds to the case where, besides M(1) and M(2), we also know the slope at the origin of the HVP-function, i.e. we know M(0). The corresponding interpolation, which in the QCD case we shall call the First Marichev Interpolation has the same functional form as the QED one in Eq. (5.9), i.e.

M(s) ⇒ M

(1,2)

(0) (s) = C 1 1 -s Γ(c -1) Γ(c -s) Γ(d -s) Γ(d -1) , ( 5

VI Tests with a Phenomenological Model

In order to test the interpolation approach proposed above, we shall apply it to a phenomenological model of the hadronic spectral function, which we call the Toy Model [START_REF] Lellouch | Private Comunication[END_REF]. The model has been constructed to simulate the basic features of the phenomenological spectral function, but with fixed parameters (i.e. no errors) so as to be handled mathematically as an exact function. The Toy Model spectral function in the low energy range below 1 GeV 2 is shown in Fig. (2). Although this Toy Model should not be confused with the experimental determination of the spectral function, beautifully shown e.g. in ref. [START_REF] Davier | Update of the Hadronic Vacuum Polarization Contribution to the muon g -2[END_REF], it reproduces nevertheless rather well its phenomenological features and, for our purposes, it can be considered as a good theoretical test laboratory. The Toy Model Spectral Function below 1 GeV 2

We first observe that, by contrast to the complex structure of the spectral function of the Toy Model shown in Fig. (2), its Mellin transform, which is shown in Fig. (3), has an extraordinarily smooth shape. In particular, the value at s = 0, which corresponds to the slope of the Π HVP (Q 2 )function at the origin and which we shall use later is M ToyM (0) = 0.7057904 × 10 -3 .

(6.1)

The value at s = -1, which corresponds to the curvature of the same Π HVP (Q 2 )-function at the origin (the second derivative) and which we shall also use later is

M ToyM (-1) = 0.1151594 × 10 -3 . (6.2)
For 0 ≤ Re(s) < 1 the Mellin transform continues to rise monotonously to become singular at s = 1 as predicted by pQCD.

Using the standard representation of a HVP µ in Eq. (3.1) we find that the Toy Model predicted value of the HVP contribution to the muon anomaly is a HVP µ (ToyM) = 6.812175 × 10 -8 , (6.3) somewhat smaller than the determinations using e + e -data in Eqs. (1.12) which, in fact, have some extra contributions included; but, as already said, the purpose of the model is not to reproduce experimental results but rather to be used as a testing theoretical laboratory. The Toy Model above has four active flavours: the three light quarks u, d, s and the heavy charm quark with mass M c . On the other hand, the Marichev interpolations that we have discussed in the previous section are for three light flavours u, d, s. Therefore, in order to compare it with the Toy At this level it is interesting to compare the Mellin transform of the Toy Model with those corresponding to the two Marichev interpolations in Eqs. (5.24) and (5.26). This is shown in Fig. (4) below. One can see a net improvement from the first interpolation (the green curve) to the second one (the blue curve) which is already quite close to the Toy Model one (the red curve). VI.2 An Application to the LQCD Prediction of ref. [START_REF] Sz | [END_REF].

Let us now apply the Marichev interpolation technique to recent LQCD results. The lattice QCD BMWc collaboration has recently published results on the first two moments M(0) and M(-1). Their numbers [START_REF] Sz | [END_REF]:

Π 1 [GeV -2 ] = 0.0999(10)(9)(23)(13) and Π 2 [GeV -4 ] = -0.181(6)(4)(10)(2) , (6.22) when expressed in the normalization (t 0 = 4m 2 π ± ) of our Mellin Moments, with the charm contribution subtracted, and with their errors added quadratically, correspond to the values: M(0) [START_REF] Sz | [END_REF] = (0.704 ± 0.021) × 10 -3 and M(-1) Mellin Transforms corresponding to the Marichev Interpolations using as an input the central values of the LQCD results of ref. [START_REF] Sz | [END_REF].

Blue is the first interpolation, Red the second one. Figure [START_REF] Th | The Muon (g -2) Theory Value: Present and Future[END_REF] shows the shapes of the real parts of the functions ℜ 1,2 0 (τ ) and ℜ 1,2 0,-1 (τ ) in Eqs. (6.11) and (6.12) using the central values of the LQCD results in ref. [START_REF] Sz | [END_REF]. It is interesting that, although these shapes differ in detail from the ones corresponding to the Toy Model in Fig. (5) and from the QED one in Fig. (1) they are qualitatively rather similar. This is probably due to the fact that they all have in common the Gamma-function structure of the Marichev general ansatz in Eq. (5.1).

The corresponding predictions for a HVP µ using the first and second interpolations, with the values of the moments in Eq. (6.23) in which the charm contributions have been subtracted, are:

. 16 )HVP µ in powers of m 2 µ

 162 Applying the converse mapping theorem to this representation, results in a series expansion of a t0 , with coefficients which are governed by the values of M(s) at s = 0, -1, -2, • • • and give positive contributions; and by the values of the first derivative of M(s):

V. 1 . 1 FFigure 1 :

 111 Figure 1: Shape of the Function R(τ ) in the integrand of Eq. (5.14) (first line).

. 20 )

 20 but with the parameters c and d restricted now to satisfy the two QCD constraints at s = 2 and s = 0: c = 2 and d = 1 + 1 C M(0) . (5.22)

  Figure 2:

Figure 4 :

 4 Figure 4: Mellin Transforms of the Hadronic Spectral Function Red: Toy Model without Charm Green: Marichev's Interpolation with only M(0) as input Blue: Marichev's Interpolation with M(0) and M(-1) as input.

Figure 5 : 1

 51 Figure 5:Shape of the Functions ℜ 1,2 0 (τ ) (green curve)and ℜ 1,2 0,-1 (τ ) (blue curve)

  Figure 6:

Figure ( 6 )

 6 Figure (6) shows the Mellin transforms obtained with the central values of the numbers above; here the blue curve corresponds to the First-Interpolation, the red curve to the improved Second-Interpolation. These curves are rather similar to the ones in Fig. (4) which test the Toy Model.

(

  Second ) = (6.81 ± 0.30) × 10 -8 .(6.26)

Notice that in ref.[START_REF] Rafael | [END_REF] the chosen normalization scale is the muon mass.

The singular expansion (or singular series) of a meromorphic function is a formal series collecting the singular elements at all poles of the function (a singular element being the truncated Laurent's series at O(1) of the function at a given pole) and it is conventionally denoted by the symbol ≍[22].

The Slater procedure applied to several examples, as well as the convergence of the Marichev interpolation, will be discussed in a forthcoming paper[START_REF] Greynat | The Mellin-Barnes Approach to Hadronic Vacuum Polarization and g µ -2[END_REF].

I am very grateful to Santi Peris for reminding me of this fact.

Acknowledgments.

I am very much indebted to David Greynat for many helpful comments and suggestions on the topics discussed in this paper; in particular for bringing to my attention Ramanujan's Master Theorem. I have also benefited from discussions with Laurent Lellouch, who has kindly provided the code of his Toy Model, and with Marc Knecht and Santi Peris. I wish to thank David Greynat and Laurent Lellouch for a careful reading of the successive versions of the manuscript.

With A denoting the quantity:

the interpolation in question is then:

(5. 24) One can now proceed to the determination of the corresponding prediction for a HVP µ inserting this

(s) interpolating Mellin transform in Eq. (5.13) and evaluating numerically the integral with e.g. the choice c s = 1 2 :

where F (s) is the function defined in Eq. (3.16).

V.2.2 Second Marichev Interpolation when M(0) and M(-1) are known

The previous interpolation can be improved once we know both the slope and the curvature of Π HVP (Q 2 ) at the origin i.e. when M(0) and M(-1) are known, in which case we use as an ansatz the following Second Marichev Interpolation:

with the e and f parameters fixed by the matching equations:

(0,-1) (-1) = M(-1) .

(5.27)

In terms of the quantities A (defined in Eq. (5.23)) and the ratio:

we find

, and e = Rf .

(5.29)

The corresponding prediction for a HVP µ is then given by the numerical evaluation of the integral (which is in fact a Fourier-like transform):

(5.30)

We wish to emphasize that the first and second Marichev functions in Eqs. (5.24) and (5.26) are unique in the sense that, with the information provided, they are the most general Mellin transforms satisfying the criteria stated after Eq. (5.1).

We have now all the ingredients to test these First and Second Marichev Interpolations with phenomenological models, and then to apply them to the evaluation of a HVP µ using as an input the recent LQCD determinations [START_REF] Sz | [END_REF] of M(0) and M(-1). Model, we have to subtract from it the charm quark c contribution. This we do by observing that the effective charm-quark contribution is well described by a constituent charm quark model with a spectral function

with (we shall use for M c its central value)

This gives a contribution to the Mellin transform:

and the Mellin transform of the Toy Model spectral function to compare with should, therefore, be: Furthermore, the contribution to the muon anomaly from the charm spectral function in Eq. (6.4) with M c = 1.275 GeV, using e.g. the standard representation in Eq. (3.1), is

which also has to be subtracted as well from the one in Eq. (6.3). More precisely, the predictions for the muon anomaly using the Marichev interpolation will be compared to the value

where for convenience we have factorized the overall constant C in Eq. (5.17). The shapes of the real parts of these two functions of τ are shown in Fig. (5). The green curve corresponds to ℜ 1,2 0 (τ ), the blue curve to ℜ 1,2 0,-1 (τ ) The result we get for a HVP µ using the First Marichev Interpolation given by Eq. (5.24), i.e. the one corresponding to the curves in green in Figs. ( 4) and ( 5), is a HVP µ (first) = 6.25021 × 10 -8 . (6.13)

It reproduces the Toy Model value at the 6.6% level. Not competitive enough for a comparison with the experimental results, but a net improvement with respect to the upper bound [27] value: and reproduces the Toy Model value at the 0.6% level. We find this very encouraging! VI.1 Test with the BHLS Model of ref. [21].

One may perhaps suspect that the reason for the success of the previous results is due to the particular choice of the Toy Model as a reference. We have, therefore, also considered another phenomenological model as an alternative reference: the so called BHLS-Model of ref. [21], and applied the same method to it. For that we choose the entries corresponding to what the authors of ref. [21] call Data Direct.

The values quoted for the first two moments are:

M(0) BHLS = (10.1307 ± 0.0745) × 10 -5 , (6.16)

M(-1) BHLS = (0.23507 ± 0.00185) × 10 -5 , (6.17 Using these moments as an input, the result for a HVP µ from the first Marichev approximation, the one which only requires M(0) as an input, is a HVP µ (first) = 626.12 × 10 -10 .

(6.20)

Using the second Marichev approximation, which requires both M(0) and M(-1) as input, we find

in agreement with the BHLS value in Eq. (6.18) at the 1% level. With the charm contribution in Eq. (6.9) subtracted to the BHLS value in Eq. (6.18) the agreement is at the 0.6% level, much the same as in the case of the Toy Model. Shape of the Functions ℜ 1,2 0 (τ ) (blue curve)and ℜ 1,2 -1,0 (τ ) (red curve) using as an input the central values of the LQCD results of ref. [START_REF] Sz | [END_REF].

The error in a HVP µ (First ) is an average of the two limits of error induced by the error in M(0). The error in a HVP µ (Second ) is sensitive to both the errors in M(0) and M(-1); it has been estimated by evaluating a HVP µ with the limit of errors in M(0) in Eq. (6.23) keeping the central value of M(-1), then evaluating a HVP µ with the limit of errors in M(-1) in Eq. (6.23) keeping the central value of M(0) and finally averaging the partial errors quadratically. We find these results very encouraging to pursue with more accurate LQCD determinations of M(0) and M(-1) and, if possible, with the determination of higher moments.

VII Conclusions.

We conclude from this work that with a precise determination of the first Mellin Moment M(0) i.e., with a precise determination of just the slope of the HVP function at the origin accessible to LQCD, one can already obtain a result for a HVP µ which provides a first rough test of the determinations using experimental data. Notice that in the First Marichev Interpolation, besides the eventual determination of M(0), one only uses as further information two well known properties of QCD: asymptotic freedom and the fact that in the chiral limit there is no 1/Q 2 term in the OPE of Π(Q 2 ). With such limited input there is no prediction from Padé approximants one can compare with.

The Second Marichev Interpolation of the Mellin Transform of the hadronic spectral function which we have developed above results in a much more accurate determination. It includes as an input the determinations of the first two moments M(0) and M(-1), i.e. the determination of the first two derivatives of the HVP function Π(Q 2 ) at Q 2 = 0 accessible to LQCD. The test with the Toy Model above results in a determination of a HVP µ with an accuracy of 0.6% which is very encouraging. The application to the determination of the M(0) and M(-1) moments from LQCD [START_REF] Sz | [END_REF] points towards a very promising future in this direction.

It would be very helpful to be able to test the Marichev Interpolation Approach with real experimental data. In that respect we encourage our colleagues of refs. [START_REF] Davier | [END_REF]9] to publish the values of a few moments: M(0), M(-1), M(-2) , • • • of the same physical spectral function which they use for their determination of a HVP µ .


Goffredo Chirco

Thibaut Josset

Statistical mechanics of covariant systems with multi-fingered time

Keywords: numbers: 05, 20, -y, 04, 20, Cv

Recently, in [Class. Quantum Grav. 33 (2016) 045005], the authors proposed a new approach extending the framework of statistical mechanics to reparametrization-invariant systems with no additional gauges. In this work, the approach is generalized to systems defined by more than one Hamiltonian constraints (multi-fingered time). We show how well known features as the Ehrenfest-Tolman effect and the Jüttner distribution for the relativistic gas can be consistently recovered from a covariant approach in the multi-fingered framework. Eventually, the crucial role played by the interaction in the definition of a global notion of equilibrium is discussed.

I. INTRODUCTION

Thermodynamics and general relativity are well understood essential pillars of modern physics. Individually, these theories provides an extremely accurate description of nature. However, when considered together, their interplay raises unexpected questions: Why are the laws of black holes mechanics equivalent to the laws of thermodynamics [START_REF] Bardeen | The four laws of black hole mechanics[END_REF]? Why is Hawking radiation thermal [START_REF] Hawking | Black hole explosions?[END_REF]? What are the degrees of freedom contributing to black hole entropy [START_REF] Bekenstein | Black holes and entropy[END_REF]? Is the Einstein equation an equation of state [START_REF] Jacobson | Thermodynamics of spacetime: The Einstein equation of state[END_REF]? Is gravity an entropic force [START_REF] Verlinde | On the origin of gravity and the laws of Newton[END_REF]?

Despite much effort, the nature of the interplay between thermodynamics and general relativity remains mostly obscure. In particular, no fully (general) relativistic statistical mechanics -taking spacetime's degrees of freedom into account-has yet been found. One of the main obstacles in this sense is due to the central role played by time in statistical mechanics. General relativity has forced us to reconsider the basic notions used in mechanics such as space, time, observable, reference frame. In particular, in a general covariant framework, the notion of time is considerably weakened 1 . This necessarily leads to a serious modification of the notion of thermodynamic equilibrium.

In [START_REF] Chirco | Statistical mechanics of reparametrization-invariant systems. It takes three to tango[END_REF] the authors discussed a relational reformulation of the notions of statistical state and equilibrium, for systems characterised by reparametrization invariance with no extra gauge symmetry. Such systems, defined by a vanishing Hamiltonian, have one-dimensional orbits, with no preferred parametrization. In other words, the time parameter appearing in the equations of motion is unphysical. In this context, the notion of statistical state was shown to rely on a suitable split of the system into two parts, one on which measurements are performed, and one used as a clock. From a relational perspective, the ideal Newtonian time of classical mechanics turned out to be just a peculiar example of "clock system".

In this paper, we move a step forward along the same line, by considering the case of reparametrizationinvariant systems with additional gauge symmetries. Such systems are defined by more than one vanishing constraint, so the orbits are multi-dimensional, their parametrization being still unphysical. For instance, a field theory like general relativity has one Hamiltonian constraint per point in space, and thus as many (unphysical) time parameters. Such systems are typically referred to as multi-fingered time systems.

We investigate the statistical mechanics of systems with multi-fingered time, making use of the formal approach defined in [START_REF] Chirco | Statistical mechanics of reparametrization-invariant systems. It takes three to tango[END_REF] on single constraint systems. More precisely, we will argue for the crucial role of interaction in the construction of the statistical mechanics of these systems.

The paper is organized as follows. Section II reviews the Hamiltonian formulation of relativistic mechanics and the statistical mechanics deduced from it, in the case of a single constraint. Section III probes a simple system with a two-fingered time, and derive the Ehrenfest-Tolman effect. Section IV recovers the statistical mechanics of a gas of special relativistic particles, starting from a completely covariant mechanical description. Finally, section V summarizes our results and briefly discusses the possible implications for gravity.

II. RELATIVISTIC SYSTEMS

Classical mechanics provides a very good description of the evolution of observables in Newtonian time, but it becomes inappropriate for relativistic systems. Therefore, generalized Hamiltonian mechanics has been developed by Dirac [START_REF] Dirac | Generalized hamiltonian mechanics[END_REF]. Following [START_REF] Rovelli | Quantum gravity[END_REF], we give here a brief summary of the formalism of generalized Hamiltonian mechanics, arXiv:1606.04444v1 [gr-qc] 14 Jun 2016 which will then be used throughout the paper.

A. Mechanics

Let us consider a relativistic system described in a n-dimensional configuration space C, coordinatized by partial observables (q α ), whose dynamics is given by a reparametrization-invariant action L(q α , qα )dτ . Any such singular Lagrangian leads to a vanishing Hamiltonian, relating the different conjugate momenta p α = ∂L ∂ qα . If there are extra gauge symmetries, one gets a set of constraints {C i = 0 , 1 ≤ i ≤ k}. We call extended phase space the 2n-dimensional symplectic space X = T * C = {(q α , p α )}, equipped with ω X = dp α ∧ dq α , and presymplectic surface the 2nk subspace Σ defined by the constraints.

The dynamics can be recast in the action p α qα -N i C i (q α , p α ) dτ , where {N i } are Lagrange multipliers. The orbits are then given by Hamilton equations

dq α dτ = N i (τ ) ∂C i ∂p α , dp α dτ = -N i (τ ) ∂C i ∂q α . (1)
Different choices for the functions {N i } corresponds to gauge-equivalent representations of the same kdimensional motion. Finally, the physical phase space Γ is simply the set of motions, and it has dimension 2n-2k (see figure 1).

X ⌃ ⇡

B. Statistical mechanics of reparametrization-invariant systems

The possibility of doing statistical mechanics of reparametrization-invariant systems, without additional gauge (k = 1), has been recently discussed in [START_REF] Chirco | Statistical mechanics of reparametrization-invariant systems. It takes three to tango[END_REF]. This new approach is based on a formulation of the ergodic hypothesis, without any reference to Newtonian time. In particular, it shows that, for a system made of two noninteracting subsystems

X = X a × X b , ω X = ω X a + ω X b ,
and

C = C a + C b = 0 , (2)
the microcanonical ensemble for one subsystem (say S b ), using the other one (here S a ) as a clock, is well defined and simply reads2

δ(C b -I) dµ X b δ(C b -I) dµ X b , (3)
where dµ X b is the Liouville measure on X b , ω X b and I indicates the value of the function C b (preserved by the dynamics). It is important to notice that (3) does not depend on any specific choice of clock in S a , nor of its dynamics (given by C a ). Given a well defined notion of microcanonical ensemble, it is then possible to apply the standard tools of statistical mechanics. More precisely, if the subsystem of interest in turn splits into weakly interacting components

X = X a × X b × X c , ω X = ω X a + ω X b + ω X c , and C = C a + C b + C c + V bc = 0, (4)
one can generalize the notions of entropy and temperature, respectively

S α (I α ) = def log δ(C α -I α )dµ X α , 1 T α I = def dS α dI α .
The equilibrium between S b and S c , with respect to S a , is then simply given by the (generalized) zeroth law

T b I = T c I . (5)
In brief, the existence of absolute Newtonian time and energy can be reformulated in a relational way, by means of a suitable split of the whole system into subsystems. Doing so naturally extends statistical mechanics and thermodynamics to reparametrization-invariant systems (cf. [START_REF] Chirco | Statistical mechanics of reparametrization-invariant systems. It takes three to tango[END_REF] for more details). Let us explore how this approach could be extended to systems with more than constraints.

III. EHRENFEST-TOLMAN EFFECT

In 1930, Tolman jointly applied the mass-energy equivalence (special relativity) and the equivalence principle (general relativity) to heat [START_REF] Tolman | On the weight of heat and thermal equilibrium in general relativity[END_REF]. He found that, for a fluid at equilibrium in a non-uniform gravitational field, the (proper) temperature is not constant3 . In a second paper [START_REF] Tolman | Temperature equilibrium in a static gravitational field[END_REF], Tolman and Ehrenfest gave the following characterization of thermodynamical equilibrium, in a static spacetime:

T ( x) -g 00 ( x) = cste, ( 6
)
where T is the temperature measured by a local observer at rest and √ -g 00 is the norm of the timelike Killing vector field in stationary coordinates. Equation ( 6) is now called Ehrenfest-Tolman effect, and has been re-derived in various manners [START_REF] Balazs | On thermodynamic equilibrium in a gravitational field[END_REF][START_REF] Ebert | Carnot cycles in general relativity[END_REF][START_REF] Rovelli | Thermal time and the Tolman-Ehrenfest effect: temperature as the speed of time[END_REF][START_REF] Haggard | Death and resurrection of the zeroth principle of thermodynamics[END_REF].

Despite the apparent conflict between ( 5) and ( 6), we will see that our formalism, when applied to multifingered time systems, naturally leads to Ehrenfest-Tolman effect.

A. Two non-relativistic systems

Let us consider two non-relativistic systems S a and S b in some small regions4 of a stationary spacetime. Let us denote τ a and τ b the (local) proper times. Assuming for the time being that S a and S b do not interact, the global system is described in the extended phase space

X = X a × X b = {τ a , p τ a , q a , p a , τ b , p τ b , q b , p b }, ω = dp τ a ∧ dτ a + dp a ∧ dq a + dp τ b ∧ dτ b + dp b ∧ dq b ,
by the two constraints

C a = p τ a + H a (q a , p a ) = 0 C b = p τ b + H b (q b , p b ) = 0 . ( 7
)
Denoting n a +1 (resp. n b +1) the dimension of the configuration space of S a (resp. S b ), the orbits of the global system are two-dimensional surfaces, and the physical phase space Γ has dimension 2(n a + n b ), as expected from Newtonian mechanics.

B. Gauge choices

It may seem unusual to describe the dynamics of a nonrelativistic system by the mean of two time observables. Let us pick a single variable t: without loss of generality, τ a = t and τ b = f (t). The symplectic form simply reads

ω = d(p τ a + f (t)p τ b ) ∧ dt + dp a ∧ dq a + dp b ∧ dq b .
Denoting p t = p τ a +f (t)p τ b the conjugate momentum to t, the reduced system is then described in the extended phase space X = {t, p t , q a , p a , q b , p b }, by the single constraint

C = p t + H a (q a , p a ) + f (t)H b (q b , p b ) = 0. ( 8
)
From the point of view of mechanics, any such gauge choice is leading to a well defined constrained system. In particular, the 2(n a + n b )-dimensional physical phase space remains unchanged, no information has been lost. On the other hand, the affine gauge choice

f (t) = αt + β (9)
is the only one leading to a time-independent Hamiltonian in [START_REF] Chirco | Statistical mechanics of reparametrization-invariant systems. It takes three to tango[END_REF]. In relational terms, it is the only choice leading to a split of the system between a (Newtonian) clock component, and the other partial observables. Thus, it is preferred from a statistical mechanics perspective, to satisfy requirement [START_REF] Hawking | Black hole explosions?[END_REF]. However, the rate α is still arbitrary. In other words, time and energy scales can still be chosen independently for the two systems S a and S b , thus it would be meaningless to compare their temperatures.

In the next paragraph, we illustrate how the interaction between S a and S b selects a preferred gauge choice.

C. Interaction and statistical state

Let us consider the modified constraint

C a = p τ a + H a (q a , p a ) + V (q a , q b , τ a , τ b ) = 0 C b = p τ b + H b (q b , p b ) = 0 .
In Newtonian space-time, it would be natural to assume that the (instantaneous) interaction V depends on the time variables only through τ aτ b . In a stationary spacetime, V depends only on τ a √ -g a 00

-τ b √ -g b 00
, where g a 00 (resp. g b 00 ) is the gravitational field seen by the subsystem S a (resp. S b ), in stationary coordinates. After gauge fixing, the single constraint (8) splits if and only if α = -g b 00 / -g a 00 , and in that case

C = p t +H a (q a , p a )+ -g b 00 -g a 00 H b (q b , p b )+V (q a , q b ) = 0.
Thus, by requiring condition (4) to be satisfied, the interaction selects a preferred gauge (see figure 2). One can now write the zeroth law (5) in terms of generalized temperatures, and finally deduce Ehrenfest-Tolman relation

T a -g a 00 = T b -g b 00
for the proper temperatures.

To summarize, it is clear that there is no way to define a temperature for the global system without considering the interaction. Ehrenfest-Tolman effect is due to the change in the weight of energy that is moved in space, from one subsystem to the other one.

S a q a τ a √ -g a 00 S b q b τ b -g b 00 Interaction V gauge fixing τ b = f (τ a )
Figure 2. Two non-relativistic systems S a and S b , evolving in a stationary, locally uniform, gravitational field, can be described either with their own (local) time observables (gauge redundancy) or using a single time variable. The interaction V selects a preferred gauge, leading to an equilibrium state which satisfies Ehrenfest-Tolman effect.

IV. RELATIVISTIC GAS

After a long controversy about the correct statistical mechanics and thermodynamics in special relativity, the generalization of Maxwell distribution initially proposed by Jüttner in 1911 [START_REF] Jüttner | Das maxwellsche gesetz der geschwindigkeitsverteilung in der relativtheorie[END_REF],

γ( v) 5 e -mγ( v)/k B T , where γ( v) = 1 √ 1-v 2 ,
seems to be confirmed by recent numerical simulations [START_REF] Cubero | Thermal equilibrium and statistical thermometers in special relativity[END_REF]. Using the momentum p = mγ( v) v, the probability distribution reads

1 Z e -β µ pµ d 3 pd 3 x, (10)
where β µ is a time-like 4-vector, reducing to 1 k B T , 0, 0, 0 in the center of mass frame.

A. Free special relativistic particles

The generalized Hamiltonian formalism presented in section II is very convenient to keep manifest the symmetries of special relativistic systems. In particular, a free particle of mass m is characterized by p µ p µ + m 2 = 0. A gas of non-interacting N such particles is simply described in the 8N -dimensional extended phase space

X = × α {x (α)µ , p (α) µ },
by the set of constraints

C (α) = p (α) µ p (α)µ + m 2 = 0, 1 ≤ α ≤ N. (11)
Motions are N -dimensional surfaces, and the physical phase space has dimension 6N .

B. Gauge choices

Let us choose a single time-like variable t, related to the coordinates of each particle by

t = f α x (α) , 1 ≤ α ≤ N.
As in section III, a generic gauge choice would lead to a time-dependent Hamiltonian (i.e. a constraint that does not split), and requirement (2) for statistical mechanics would not be satisfied.

As long as there is no interaction, it is possible to choose different Lorentz frames for different particles: t = x (α)0 , where x (α) = Λ(α)x (α) . In that case, the single constraint is

C = p t + α p (α) 2 + m 2 = 0, (12)
where p

(α) i = [Λ(α)] i µ p (α) µ
is the conjugate momentum to x (α)i . The constraint [START_REF] Einstein | Zur theorie des statischen gravitationsfeldes[END_REF] splits into a clock part and a time-independent Hamiltonian. However, as the Lorentz transformations {Λ(α)} can be chosen arbitrarily, the temperature of the gas is ill-defined.

C. Gas of weakly interacting particles

Though, the picture changes when a coupling between the particles is turned on.

Let us consider a (scattering) interaction term5 of the form

V x (α)µ -x (β)µ x (α) µ -x (β) µ
. It is clear that the interaction part of the Hamiltonian constraint is timeindependent only if one chooses the same Lorentz frame for all the particles. So, without loss of generality, let us take Λ(α) = I. Then, [START_REF] Einstein | Zur theorie des statischen gravitationsfeldes[END_REF] reduces to the classical formulation of special relativistic particles:

C = p t + α p (α) 2 + m 2 + α =β V [ x(α) -x(β) ] = 0,
where V insures thermalization of the gas, but remains small compared to the kinetic energy. For any value of p t = -E (and total 3-momentum6 P ), one can write the microcanonical state (3), and deduce Jüttner distribution [START_REF] Rovelli | Quantum gravity[END_REF].

The equilibrium statistical state is not selected by the Lorentz invariance of each free constraint [START_REF] Tolman | On the weight of heat and thermal equilibrium in general relativity[END_REF], but rather by the Lorentz invariance of the interaction term V .

V. DISCUSSION

It is known since Boltzmann's work that a weak, mixing, interaction is necessary for reaching equilibrium. However, in general, the details of the interaction do not matter much, and its analysis is often relegated to the ergodic hypothesis. In the case of relativistic systems described by several constraints, the role of the interaction is dramatically enhanced, as it entails whether the constraint can be split or not. Therefore, the interaction determines the subsystems for which it is meaningful to talk about equilibrium, and the corresponding clock systems. This led us to recover well known features (Ehrenfest-Tolman effect and Jüttner distribution) from a covariant approach.

The interplay between interaction and gauge has been studied in [20], and its implications for statistical mechanics were suggested in the context of thermal time hypothesis [START_REF] Chirco | Coupling and thermal equilibrium in general-covariant systems[END_REF]. In this work, we have taken a step further by showing how equilibrium states are selected in the whole set of possible statistical states.

Trying to apply this formalism to general relativity, one faces a technical problem: how to find a nice gauge and identify a split of the system. This issue is very similar to the problem of time in quantum gravity. However, in the context of statistical mechanics one may relax requirement (4) into

C = C a + C b + C c + V bc + V abc = 0,
where V abc V bc . In other words, the only physically reasonable assumption is that the different components are "thermalizing" faster than exchanging "energy" with the clock system.

VI. ACKNOWLEDGMENTS

Figure 1 .

 1 Figure 1. Generalized Hamiltonian mechanics.

The "problem of time" plays a central role for several distinct issues[START_REF] Kuchař | Time and Interpretations of Quantum Gravity[END_REF][START_REF] Rovelli | Analysis of the distinct meanings of the notion of time, in different physical theories[END_REF], from the choice of an internal time variable in quantum gravity, to the emergence of a directional global time in cosmology.

If C b is not the only conserved quantity for S b , the statistical state should include additional "δ-functions".

This fact was actually glimpsed by Einstein[START_REF] Einstein | Zur theorie des statischen gravitationsfeldes[END_REF], even before General Relativity was completed!

So that the local gravitational field is approximately uniform.

V can be thought as an effective potential coming from a (Lorentz invariant) interaction with an (relativistic) ambient field.

In the absence of external potential, the total 4-momentum is conserved, insuring Lorentz covariance of the statistical state.

We are grateful to Carlo Rovelli for careful reading of the draft and useful comments.