Stephan Narison
email: snarison@yahoo.fr

Di-Gluonium Sum Rules, I = 0 Scalar Mesons and Conformal Anomaly

Keywords: QCD Spectral Sum Rules, Perturbative and Non-perturbative QCD, Exotic hadrons, Masses and Decay constants

.

Introduction

Gluonia / Glueballs bound states are expected to be a consequence of QCD [5]. However, despite considerable theoretical and experimental efforts, there are not (at present) any clear indication of their signature.

The difficulty is also due to the fact that the observed candidates can be a strong mixing of the gluonia with the qq light mesons or some other exotic mesons (four-quark, hybrid states). In particular, the (I = 0) Isoscalar Scalar Channel channel is overpopulated beyond the conventional scalar nonet such that one suspects that some of these states can be of exotic origins and may (for instance) contain some large gluon component in their wave functions.

In some previous works [1][2][3][4][6][7][8][9][10][11][12][13] based on QCD spectral sum rules (QSSR) à la SVZ [14][15][16] combined with some low-energy theorems (LET), we have tried to understand the (pseudo)scalar gluonia channels.

Since then, some progresses have been accomplished from experiments and from some other approaches which we shall briefly remind below.

After this short reminder, we update and improve our previous QSSR analysis. The impact of these novel results on our present understanding of the peculiar (I = 0) isoscalar scalar channel is discussed.

2. σ/f 0 (500) and f 0 (980) from scatterings data To quantify the observables of the low-lying scalar states, we have phenomenologically studied the ππ, KK, γγ [17,18], K e4 [19], J/ψ, φ radiative and D s semileptonic decays [20,21] data with the aim to understand the internal substructure of the σ/f 0 (500) and f 0 (980) mesons, 1 .

• Masses and hadronic widths from ππ → ππ, KK and Ke 4 data Using a K-matrix [26,27] analysis of the ππ elastic scattering data below 0.7 GeV, we found for the complex pole mass and the residue using one bare resonance ⊕ one channel [17] :

M pole σ [MeV] 422 -i290. (1)
We extend the previous analysis until 900 MeV and use two bare resonaces [σ(500), f 0 (980)] ⊕ two open channels (ππ, KK) parametrization of the data [18,19]. To the data used previously, we add the new precise measurement from NA48/2 [START_REF] Bloch-Devaux | The NA48/2 collaboration[END_REF] on Ke 4 (K → ππeν e ) data for the ππ phase shift below 390 MeV.

Then, averaging this result with the previous one in Eq.1, one obtains the final estimate :

M pole σ
[MeV] 452(12) -i260 (15) , |g σπ + π -| = 1.12(31) GeV,

|g f0K + K -| |g f0π + π -| = 2.58(1.34) , (2)
which agrees with the ones based on the analytic continuation and unitarity properties of the amplitude to the deep imaginary region [29,30] :

M pole σ
[MeV] 444 -i272 , and M pole σ

[MeV] 489 -i264 .

(3)

1 A more complete and comprehensive discussion on production processes and decays of scalar mesons can be found in the recent reviews of [22][23][START_REF] Amsler | PDG Review on scalar mesons below 2 GeV and Review on nonqq mesons in[END_REF][START_REF] Zyla | [END_REF].

2 However, in order to compare these results with the QCD spectral sum rules ones where the analysis is done in the real axis, one has to introduce the On-shell or Breit-Wigner (os) mass and width where the amplitude is purely imaginary at the phase 90 0 :

ReD((M os σ ) 2 ) = 0 =⇒ (M os σ , Γ os σ ) = (920, 700) MeV ,

where D is the propagator appearing in the unitary ππ amplitude. Similar result has been obtained in [31,32] where the wide resonance centered at 1 GeV has been interpreted in [31] as a resonance interfering destructively with the f 0 (980) and f 0 (1500).

For the f 0 (980), one obtains [18,19]: (34) -i18 (11) , |g f0π + π -| = 2.65 (10) GeV,

M pole f0 [MeV] 981
|g σK + K -| |g σπ + π -| = 0.37(6) , (5)
and

M os f0 M pole f0 , Γ os f0 Γ pole f0 , (6)
where Eq. 6 comes from the fact that f 0 (980) is narrow.

• σ/f 0 (500) and f 0 (980) γγ widths from γγ → ππ, KK scatterings

We proceed as above and deduce the average (in units of keV):

Γ dir σ pole = 0.16(4) , Γ resc σ pole = 1.89 (81) , Γ tot σ pole = 3.08(82) (7) for the direct, rescattering and total (direct ⊕ rescattering) γγ widths. For the on-shell mass, these lead to :

Γ dir σos = 1.2(3) keV, Γ resc σ pole = Γ resc σos .

For the f 0 (980), one obtains in units of keV :

Γ dir f0 | pole Γ dir f0 | os = 0.28 (1), Γ res f0 = 0.85 (5), Γ tot f0 = 0.16 (1) .

• Production of σ/f 0 (500) and f 0 (980) from some other channels Isoscalar Scalar production in some other channels have been reviewed in [22][23][START_REF] Amsler | PDG Review on scalar mesons below 2 GeV and Review on nonqq mesons in[END_REF][START_REF] Zyla | [END_REF], where the structures:

σ/f 0 (500) , f 0 (980) , f 0 (1370) , f 0 (1500) (10) have been also found in different processes (J/ψ-decays, central production in pp and e + e -, pp and pn annihilation at rest, πp and Kp).

We notice the production of the σ/f 0 (500) from J/ψ → ωπ + π - [34] with the extracted complex pole mass and width:

M pole σ = 541(39) -i252(42) , (11)
and of the f 0 (980) from J/ψ → γπ 0 π 0 by BESIII [33] and from J/ψ → φπ + π -and J/ψ → φK + K - [34] with a ratio of coupling :

g f0K + K -/g f0π + π -= 4.21 ± 0.32 . (12)
The f 0 (980) is also produced from φ → (π + π -, K + K -) + γ, a glue filter channel at KLOE [35] with the couplings :

g f0K + K -= (3.97 ∼ 4.74) GeV , g f0π + π -= -(2.22 ∼ 1.82) GeV , (13)
where in these last two experiments the presence of the σ/f 0 (500) improves the quality of the fit. The non-vanishing of g f0π + π -indicates that the f 0 (980) cannot be a pure ss state.

Production of some other (I = 0) isoscalar scalar mesons

To these data, we add the production of the new structures which can have gluon component [START_REF] Zyla | [END_REF] :

f 0 (1765) , f ) (2100) , f 0 (2020), f 0 (2200) . ( 14
)
observed from J/ψ → γπ + π -, K + k -, 4π radiative decays, πp, pp, pp → ηη , ηη, 4π [START_REF] Zyla | [END_REF] and recently from

η c → [f 0 (1.7) → K + K -], f 0 (2.1) → π + π -]η by BABAR [36].

Impacts of the previous phenomenological analysis and data

One can notice that the I = 0 scalar channel is overpopulated and goes beyond the usual nonet expectations. From the previous phenomenological analysis, one may expect that:

• σ/f 0 (500) meson
It cannot be a pure four-quark (ūu + dd) or / and ππ molecule state as it has non-vanishing coupling to KK.

However, from the large value of its γγ rescattering width (see e.g. [17]), one may be tempted to interpret it as ππ or/and KK molecule state but the KK molecule mass is expected to be around 1 GeV due to SU (3) breakings.

From the size of its direct on-shell width of about 1 keV to γγ [17], it caanot be a pure (ūu + dd) state which has a γγ-width expected to be about 4 keV from QCD spectral sum rules [11] and quark model [START_REF] Babcock | [END_REF] but larger than for a pure gluonium state of about (0.2-0.6) KeV [2,17]. It cannot also be a four-quark state as the model predicts much smaller (0.00..keV) direct γγ widths [10,38].

• f 0 (980) meson
A large ss component would lead to a direct γγ width of 0.4 keV which are comparable with the data but its non-zero coupling to ππ rules out this possibility. An estimate of a mass of the ss state from QSSR leads to a value around 1.4 GeV [11] which is (relatively) too high for the f -0(980).

• f 0 (1.37) meson
The large decay width of this meson to 2(ππ) in a S-wave may signal a large gluon component in its wave functions [1,2].

• f 0 (1.37) and f 0 (1.5) mesons

Its decay into η η and ηη and the presence of a 2(ππ) S indicate a large glun component in its wave function.

Theoretical expectations

• In the I = 0 isoscalar scalar channel, we have [1,2] explored the possibility that the σ/f 0 (500) meson can be the lowest scalar gluonium state σ B (hereafter, the subindex B refers to unmixed glyonium state) having a mass around 1 GeV. From our approach the σ is composed mainly with gluons and it couples strongly and universally to ππ and KK explaining its large width while its "direct" decay to γγ is expected to be about (0.2-0.6) keV .

• A maximal mixing of this gluonium with a (ūu + dd) state may explain the feature of the observed σ/f 0 (500) and f 0 (980) : narrow width in ππ but strong coupling to KK [4] 2 .

• We have also argued [1,2] that the contribution of the σ B to the subtracted (SSR) and unsubtracted (USR) sum rules is necessary for resolving the apparent inconsistency between these two sum rules where, in addition, a higher mass gluonium G(1.5) is needed which we have identified [1] with the G(1.6) found from the GAMS data [39]. The contribution of σ B to the SSR compensates the large contribution of the two-point subtraction constant to the SSR without appealing to more speculative direct instanton non-perturbative contributions.

• Another striking feature of the σ B is its analogy with its chiral partner the η which plays a crucial role for the U (1) A anomaly by its contribution to the topological charge (subtraction constant of the U (1) A twopoint correlator) [6,[42][43][44][45]. The σ B , being the dilaton particle of the conformal anomaly [1,46,47,67,68],

is expected to be associated to the trace of the energy-momentum tensor θ µ µ :

θ µ µ = 1 4 β(α s )G µν a G a µν + (1 + γ m (α s )) u,d,s m i ψi ψ i , (15)
with : γ m = 2a s + • • • is the quark mass anomalous dimension and a s ≡ α s /π. β(α s ) is the β-function normaized as :

β(α s ) = β 1 a s + β 2 a 2 s + β 3 a 3 s + • • • (16)
with:

β 1 = - 11 2 + n f 3 , β 2 = - 51 4 + 19 12 n f , β 3 = 1 64 -2857 + 5033 9 n f - 325 27 n 2 f ( 17
)
for n f number of quark flavours. For n f = 3 used in this paper, one has :

β 1 = -9/2, β 2 = -8, β 3 = -20.1198 . ( 18
)

Prospects

In this paper, we shall :

• Use selected low and high degree moments of QCD spectral sum rules (QSSR) within the standard SVZ expansion without instantons and parametrize the spectral function beyond the minimal duality ansatz "one resonance ⊕ QCD continuum" in order to extract the masses and decay constants of the scalar gluonia.

• Extract the value of the conformal charge ψ G (0) (value of the two-point correlator at zero momentum) in order to test the Low Energy Theorem (LET) [68] estimate.

• Discuss some phenomenological implications of our results for an attempt to understand the complex spectrum of the I = 0 observed scalar mesons.

The QCD anatomy of the two-point correlator

We shall work with the gluonium two-point correlator :

ψ G (q 2 ) = 16i d 4 x e iqx 0| θ µ µ G (x) θ µ µ † G (0)|0 (19)
built from the gluon component of the trace of the energy-momentum tensor in Eq. 15.

• The standard SVZ-expansion

Using the Operator Product Expansion (OPE) à la SVZ , its QCD expression can be written as :

ψ G (q 2 ) = β 2 (α s ) 2 π 2 0,1,2,... C 2n O 2n . ( 20
)
where C 2n is the Wilson coefficients calculable perturbatively while O 2n is a short-hand notation for the non-perturbative vacuum condensates 0|O 2n |0 of dimension 2n.

The unit perturbative operator (n = 0) Its contribution reads :

C 0 ≡ -Q 4 L µ C 00 + C 01 L µ + C 02 L 2 µ
with :

C 00 = 1 + 659 36 a s + 247.480a 2 s , C 01 = -a s 9 4 + 65.781a s , C 02 = 5.0625a 2 s , (21)
where the NLO (resp. N2LO) contributions have been obtained in [52] (resp. [53]) and has been adapted for a sum rule use in [54].

L µ ≡ Log(Q 2 /µ 2 )
where µ is the subtraction point. We shall use for 3 flavours : Λ = 340(28) MeV (22) deduced from α s (M Z ) = 0.1182 (19) from M χ c0,b0 -M ηc,η b mass-splittings [61,62], τ -decays [55,56] and the world average [START_REF] Zyla | [END_REF]57]. We shall use the running QCD coupling to order α 2 s :

a s (µ) = a (0) s 1 -a (0) s β 2 β 1 LL µ + a (0) s 2 β 2 β 1 2 LL 2 µ -LL µ -1 + β 3 β 1 , (23)
where :

a (0) s ≡ 1 -β 1 Log (µ/Λ) and LL µ ≡ Log 2 Log (µ/Λ) . (24)
The dimension-four gluon condensate (n = 2)

Its contribution reads :

C 4 O 4 ≡ (C 40 + L µ C 41 ) α s G 2 : C 40 = 2πa s 1 + 175 36 a s , C 41 = - 9 2 πa 2 s . (25)
We shall use the value :

α s G 2 = (6.35 ± 0.35) × 10 -2 GeV 4 (26)
determined from light and heavy quark systems [60][61][62].

The dimension-six (n = 3) gluon condensate

Its contribution reads (see [63] for the α s correction for n f = 0 and where a 1/π misprint is corrected):

C 6 O 6 = C 60 + L µ C 61 g 3 f abc G a G b G c /Q 2 : C 60 = a s , C 61 = - 29 4 a 2 s ( 27
)
with [64][65][66] :

g 3 f abc G a G b G c = (8.2 ± 1.0) GeV 2 α s G 2 , (28)
which notably differs from the instanton liquid model estimate [67][68][69] used in [2].

g 3 f abc G a G b G c ≈ (1.5±0.5) GeV 2 α s G 2 [14,
The dimension-eight (n = 4) gluon condensate

Its contribution reads :

C 8 O 8 = C 80 α 2 s G 4 /Q 4 : C 80 = 4πα s (29)
with :

α 2 s G 4 ≡ 14 α s f abc G a µρ G b νρ 2 -α s f abc G a µν G b ρλ 2 (0.55 ± 0.01) α s G 2 2 (30)
from a modified factorization with 1/N 2 c corrections (13/24) [70] and factorization (9/16) [67]. (0.55 ± 0.01) measures the deviation from factorization which has been found to be largely violated for the four-quark condensates [41,56,71].

• Beyond the standard SVZ-expansion

The tachyonic gluon mass (n = 1) To these standard contributions in the OPE, we can consider the one from a dimension-two tachyonic gluon mass contribution introduced by 3 which phenomenologically parametrizes the large order terms of the perturbative QCD series [74]. Its existence is supported by some AdS approaches [75][76][77]. This effect has been calculated explicitly in [78] :

C 2 O 2 = C 21 L µ λ 2 Q 2 : C 21 = 3 a s (31)
where λ 2 is the tachyonic gluon mass determined from e + e -→ hadrons data and the pion channel [78][79][80]:

a s λ 2 -(6.0 ± 0.5) × 10 -2 GeV 2 . ( 32
)
The direct instanton (n ≥ 5/2)

In an instanton liquid model [68,69], the direct instanton contribution is assumed to be dominated by the single instanton-anti-instanton contribution via a non-perturbative contribution to the perturbative Wilson coefficient [81,82]:

ψ(Q 2 )|Ī -I = 32β 2 Q 4 ρ 4 K 2 ρ Q 2 2 dn(ρ) (33)
where K 2 (x) is the modified Bessel function of the second kind. At this stage this classical field effect is beyond the SVZ expansion where the later assumes that one can separate without any ambiguity the perturbative Wilson coefficients from the non-perturbative condensate contributions.

Besides the fact that its contributes in the OPE as 1/Q 5 , i.e it acts like other high-dimension condensates not taken into account in the OPE, the above instanton effect depends crucially on the (model-dependent) overall density n = ∞ 0 dρ n(ρ) and on its average size ρ = (1/n) ∞ 0 dρ ρ n(ρ) which (unfortunately) are not quantitatively under a good control (ρ ranges from 5 [68] to 1.94, [65,66] and 1.65 GeV -1 [69], while n ≈ (0.5 ∼ 1.2) fm -4 ). They contribute with a high power in ρ to the spectral function Imψ(t), which can be found explicitly in [81,82], and behave as :

Imψ G (t)|Ī -I t→∞ ∼ n ρ √ t -5 t→0 ∼ n ρ √ t 4 . (34)
As such effects are quite inaccurate and model-dependent, we shall not consider them explicitly in the analysis. Instead, an eventual deviation of our results within the standard SVZ expansion from some experimental data or/and or some other alternative estimates (Low-Energy Theorems (LET), Lattice calculations,...) may signal the need of such (beyond the standard OPE) effects in the analysis. One should mention that the approach within the standard SVZ OPE and without a direct instanton effect used in the :

-U (1) A channel has predicted successfully the value of the topological charge, its slope and the η -mass and decay constant [6,13].

-Pseudoscalar pion and kaon channels have reproduced successfully the value of the light quark masses where we have also explicitly shown [83] that the direct instanton effect induces a relatively small correction contrary to some other strong claims [84] in the literature.

The Inverse Laplace transform sum rules

From its QCD asymptotic behaviour ∼ (-q 2 ) 2 Log(-q 2 /µ 2 ) one can write a twice subtracted dispersion relation:

ψ G (q 2 ) = ψ G (0) + q 2 ψ G (0) + q 4 π ∞ 0 dt t 2 Imψ G (t) (t -q 2 -i ) . (35)
• Following standard QSSR techniques [14,16], on can derive from it different form of the sum rules. In this paper, we shall work with the Exponential or Borel [14,85,86] or Inverse Laplace transform [87] finite energy sum rule(LSR) 4 :

L c n (τ ) = tc 0 dt t n e -tτ 1 π Imψ G (t) : n = -1, 0, 1, 2, 3 (36)
and the corresponding ratios of sum rules :

R c n+l n (τ ) ≡ L c n+l (τ ) L c n (τ ) , ( 37
)
where τ is the Laplace sum rule variable. In the duality ansatz :

1 π Imψ G (t) = 2 G f 2 G M 4 G δ(t -M 2 G ) + θ(t -t c )"QCD continuum" , (38)
where the f G are the lowest resonances couplings normalized as f π = 93 MeV and M G their masses while the "QCD continuum" comes from the discontinuity Imψ G (t)| QCD of the QCD diagrams from the continuum threshold t c . In the "one narrow resonance ⊕ QCD continuum" parametrization of the spectral function:

R c n+l n (τ ) M 2 G . (39)
• To get L -1 , we find convenient to take the Inverse Laplace transform of the first superconvergent 2nd derivative of the correlator. It reads explicitly:

L c -1 (τ ) = β 2 (α s ) 2 π 2 τ -2 n=0,2,••• D -1 n + ψ G (0) , (40)
with:

D -1 0 = C 00 + 2C 01 (1 -γ E -L τ ) + 3C 02 1 -π 2 /6 + (-1 + γ E + L τ ) 2 (1 -ρ 1 ) D -1 2 = C 21 λ 2 τ (1 -ρ 0 ), D -1 4 = -C 40 - C 41 36 55 + 6(γ E + L τ ) α s G 2 τ 2 , D -1 6 = -C 60 -C 61 (-1 + γ E + L τ ) g 3 f abc G a G b G c τ 3 , D -1 8 = C 80 2 (1.1 ± 0.5) α s G 2 2 τ 4 , (41)
with :

L τ ≡ -Log(τ µ 2 ) , ρ n ≡ e -tcτ 1 + t c τ + • • • (t c τ ) n n! (42)
is the QCD continuum contribution from the discontinuity of the QCD diagrams to the sum rule. We omit some eventual continuum contributions from the non-perturbative condensates which are numerically tiny.

ψ G (0) is the value of the two-point correlator at zero momentum and can be fixed by a LET as [67,68]:

ψ G (0) - 16 π β 1 α s G 2 = (1.46 ± 0.08) GeV 4 , (43)
which we shall test later on.

• To get L 0 , we take the Inverse Laplace transform of the 1st superconvergent 3rd derivative of the two-point correlator. In this way, we obtain :

L c 0 (τ ) = β 2 (α s ) 2 π 2 τ -2 n=0,2,••• D 0 n , (44)
with:

D 0 0 = 2C 00 -2C 01 (3 -2γ E -2L τ ) -6C 02 1 -3γ E + γ 2 -π 2 /6 + (-3 + 2γ E )L τ + L 2 τ (1 -ρ 2 ) D 0 2 = - C 21 2 λ 2 τ (1 -ρ 1 ), D 0 4 = -C 41 α s G 2 τ 2 , D 0 6 = C 60 -C 61 (γ E -L τ ) g 3 f abc G a G b G c τ 3 , D 0 8 = C 80 (0.55 ± 0.01) α s G 2 2 τ 4 , (45)
• The other higher degrees sum rules L c n (τ ) for n ≥ 1 can be deduced from the n th τ -derivative of L c 0 (τ ):

L c n (τ ) = (-1) n d n dτ n L c 0 (τ ) . (46)
• These superconvergent sum rules obey the homogeneous renormalization group equation (RGE):

- ∂ ∂t + β(α s )α s ∂ ∂α s L c n (e t τ, α s ) = 0 , (47)
with t ≡ (1/2)L τ with the renormalization group improved (RGI) solution :

L c n (e t τ, α s ) = L c n (t = 0, ᾱs (τ )) , ( 48
)
where ᾱs is the QCD running coupling.

• In the following analysis, we shall work with the family of sum rules having degrees less or equal to 4. Then, we shall select the sum rules which present stability (minimum or inflexion point) in the sum rule variable τ and in the continuum threshold t c such that we can extract optimal information from the analysis.

9. τ and t c behaviours of the QCD side of the LSR L c n Before doing the phenomenological analysis of the LSR, we study the τ and t c behaviour of their QCD expressions. For instance, we show the subtracted (SSR) L c -1 and unsubtracted (USR) L c 0 sum rules in Figs 1 and2, where we notice that the subtraction constant shifts the USR minimum (optimization point) at lower values of τ . This feature has leaded to the inconsistencies of the gluonium mass from these sum rules within a"one resonance ⊕ QCD continuum" parametrization of the spectral function as noted in [1,2] which can be cured by working instead with "two resonances ⊕ QCD continuum"

• The other USR moments L c 2 and L c 4 stabilize for τ of about 1.5 GeV -2 while L c 1 and L c 3 do not present τ -stability. To understand these τ -behaviours :

We write explictly the OPE expression of L 0 and L 1 to lowest order but includes the α s correction of the D = 4 contribution as this contribution vanishes to LO. Normalized to (2/π 2 )β 2 , they read:

L 0 = τ -3 1 + 9 2 πa 2 s α s G 2 τ 2 + a s g 3 f abc G a G b G c τ 3 + 4πα s α 2 s G 4 τ 4 , L 1 = 3τ -4 1 - 9 6 πa 2 s α s G 2 τ 2 -4πα s α 2 s G 4 τ 4 , (49)
where the changes of sign of the condensate contributions destabilize L 1 .

The inclusion of the huge PT radiative corrections shows that they also participate to the τ -behaviour of the QCD expression due to the non-trivial reorganization of these PT terms in each Laplace transformed moment.

• From this explicit analysis, one expects that, a priori:

One may extract reliably the decay constant of the gluonium from L c 1 . However, we shall see later on that the extraction of the decay constant from L c 3 can present τ -stability due to the exponential weight e M 2 G τ brought by the resonance contributions which compensates the decrease of the QCD expression

[L c n ∼ τ -(n+3)
] for increasing τ . In the specific example of L c 2 , we shall that the exponential factor shifts the τ -minimum to lower τ -values (Fig. 10) at which the PT corrections are much smaller than in Fig. 4, while for L c

3 , there appears a plateau inflexion point (Fig. 11) which render reliable the extraction of f G1 .

The mass of the gluonium can be reliably obtained from the ratios R c nl of the USR L c n,l which optimize at about the same value of τ and where the PT corrections tend to compensate. 10. The masses of the two lightest gluonia σ B and G 1

Ï ‡ 0.2 0.4 0.6 0.8 1.0 1.2 0.5 1.0 1.5 2.0 2.5 3.0 t @GeV -2 D L -1 @GeV 4 D Ú • 6 ‡ 5 4 Ï 3 L= 340 MeV, <aS G 2 >=6.35x10 -2 GeV 4 tc@GeV 2 D, <g 3 G 3 >=8.2 GeV 2 <aS G 2 >
t @GeV -2 D L 0 @GeV 6 D Ú • 6 ‡ 5 4 Ï 3 L= 340 MeV, <aS G 2 >=6.35x10 -2 GeV 4 tc@GeV 2 D, <g 3 G 3 >=8.2 GeV 2 <aS G 2 >
In order to get the lowest ground state gluonia masses, we shall work with different forms of USR. We shall not use for this case the observables involving L -1 in order to avoid the dependence of the result on the input value of ψ G (0) which we shall test later on. Among the possible ratios of sum rules, we have selected R c 20 and R c 42 which appear to give the most stable results versus the variation of different parameters.

• Mass of the lightest ground state (hereafter named σ B ) from R c 20

In so doing, we work with the ratio of sum rules R c 20 . We show its τ -behaviour for different values of t c in Fig. 8.

The case : one resonance ⊕ QCD continuum

We show in Fig. 7 the result with one resonance ⊕ QCD continuum parametrization of the spectral function. The τ -stability starts from t c = 2.5 GeV 2 while t c -stability is reached around 5.5 GeV 2 . We deduce in this range:

M σ B = 1088(2) τ (54) tc (57) Λ (1) G 2 (2) G 3 = 1088(79) MeV, ( 50
)
where the errors are mainly due to Λ and t c .

The case : two resonances ⊕ QCD continuum

We show in Fig. 8 the effect of a 2nd resonance with the parameters in [2] :

M G1 = 1.5 GeV and f G1 = 0.39 GeV. ( 51
)
It affects the σ B mass result by a negligible amount of about 22 MeV. We deduce the conservative result from t c = 3 to 5.5 GeV 2 and use an input error of 0.1 GeV -2 for the localization of τ :

M σ B = 1085(5) τ (59) tc (89) Λ (1) G 2 (2) G 3 (1) λ 2 (39) M G (53) f G = 1085(126) MeV, (52)
where the errors due to the non-perturbative condensates are regligible.

• Mass of the medium ground state gluonium G 1 from R c 42 Among the different combinations of sum rules, we find that R c 42 can give a reliable prediction of the medium ground state G 1 because L c 2 and L c 4 stabilize at about the same value of τ

(1.2 ∼ 1.5) GeV -2
(see Figs 4 and6). We show in Fig. 9

t @GeV -2 D M G 1 =HR 42 L 1ê4 @GeVD Ú 7 6 ‡ 5 4 Ï 3 <a S G 2 >=6.35x10 -2 GeV 4 , <g 3 G 3 >=8.2 GeV 2 <a S G 2 > t c @GeV 2 D, L= 0

The case : two resonances ⊕ QCD continuum

We obtain for t c = 3 to 6 GeV 2 :

M G1 = 1524(36) tc (7) τ (115) Λ (2) G 2 (2) G 3 (0) λ 2 (6) Mσ (11) fσ = 1524(121) MeV , (53)
where the error comes mainly from the one of Λ. We have used f σ B = 0.47 as anticipated from the next section.

The case : one resonance ⊕ QCD continuum A similar analysis leads to:

M G1 = 1515(123) MeV ( 54
)
indicating that the σ B effect to the value of the mass is quite small (about 9 MeV).

Decay constants f σB and f G1

To extract these decay constants, we work simultaneously with L c 2 , L c 3 and L c -1 and use an iteration procedure. The two former are expected to be more sensitive to f G1 while the third to f σ B .

We shall take M σ B from Eq. 52 and M G1 = 1.548 GeV anticipated from Eq. 53. We start the iteration by using f σ B =800 MeV from [1,2]. Then, we compare the value of f G1 from L c 2 , L c 3 . Then, we use the common solution of f G1 from the two sum rules into L c -1 to extract f σ B . We continue the iterations until we reach convergent solutions. The analysis corresponding to the final value of the decay constants are shown in Figs. 10 to 13.

-In Figs. 10 and reffig:fg12, we show the τ and t c behaviours of f G1 given the value of f σ B . -In Fig. 10, we show the results from L c 2 , L c 3 at the τ minimas or plateau from which we deduce the common solution for f G1 for fixed f σ B which allows to fix :

t c (3.3 ± 0.3) GeV 2 . ( 55
)
-In Fig. 13, we show the τ and t c behaviours of f σ B where we take the value τ = 0.5 GeV -2 at the inflexion point which is more pronounced for higher values of t c .

Then, we obtain :

f G1 = 394(45) tc (2) τ (22) Λ (0) G 2 (1) G 3 (23) Mσ B (9) fσ B (94) M G 1 = 394(109) MeV from L c 3 f σ B = 563(36) tc (64) τ (81) Λ (9) G 2 (11) G 3 (10) Mσ (3) M G 1 (107) f G 1 = 563(154) MeV from L c -1 (56)
One can note that value of f G1 agrees perfectly with the one in [2] while f σ B is in the range of the one obtained in [1,2]. 12. First radial excitation σ B of σ B

0.2 0.4 0.6 0.8 1.0 0.0 0.5 1.0 1.5 2.0 t @GeV -2 D f G 1 @GeVD Ú 4 4 ‡ 3.8 3.5 Ï 3.3 <aS G 2 >=6.35x10 -2 GeV 4 , <g 3 G 3 >=8.2 GeV 2 <aS G 2 > tc@GeV 2 D, L= 0
The σ B mass

We attempt to extract its mass from R c 20 which we have used to estimate M σ B . In so doing, we subtract the contributions of the σ B and G 1 at take the result for very large values of t c i.e by replacing the QCD continuum contribution by the σ B . In this way, we obtain the result in Fig. 14 We obtain for t c (6 ∼ 8) GeV 2 : at the τ -minimum :

M σ B = 1121(25) tc (10) τ (110) Λ (0) G 2 (3) G 3 (1) Mσ B (12) fσ B (7)M G1 (23) f G 1 = 1124(117) MeV, (57)
where the contributions of σ B and G 1 have been included after iteration from 2 (M σ B = 1124) to 3 (value in Eq. 57) resonances values of the σ B and G 1 masses and decay constants.

The σ B decay constant 13 GeV

Ï ‡ 0.2 0.4 0.6 0.8 1.0 0.5 1.0 1.5 2.0 t @GeV -2 D f s' B @GeVD ‡ • 8 Ï 6 L= 0.34 , <aS G 2 >=6.35x10 -2 GeV 4 tc@GeV 2 D, <g 3 G 3 >=8.2 GeV 2 <aS G 2 >, Ms' B =1.
We extract the decay constant from L c -1 using the same range of t c as for M σ B and taking τ 0.4

GeV -2 . We obtain:

f σ B = 648(75) tc (43) τ (124) Λ (20) G 2 (7) G 3 (24) Mσ B (137) fσ B (0)M G1 (56) f G 1 (33) M σ B = 648(216) MeV from L c -1 . (58)
Anticipating the estimate in Section 20 by extracting the decay constant from f 0 (1.37) → 2(ππ) S data through a LEV-SR, we shall work instead with the more precise value:

f σ B = 329(30) MeV. ( 59
)

Iterative improvments of the previous estimates

In the following, we study the effect of the σ B which is below M G1 on the previous estimates

• The mass M σ B
We include the contribution of the σ B which is below the G 1 mass. Using R c 20 as in previous section, the curves are similar to the ones in Fig. 8. The mass in Eq. 52 becomes :

M σ B = 1070(126)(1) f σ B (1) M σ B MeV , (60)
which is our final estimate.

• The mass M G1

Including the contribution of the σ B into R c 42 , the mass in Eq. 53 becomes :

M G1 = 1548(118)(15) f σ B ( 13
) M σ B = 1548(120) MeV , (61)
which is our final estimate.

• The decay constants f σ B and f G1

We repeat the procedure used previously by extracting simultaneously f σ B and f G1 from L c -1 and from L c 2 , L c 3 by including the σ B contribution. The value of t c at which the solutions from L c 2 , L c 3 meet is sligthly shifted from 3.3 to 3.2 GeV 2 . We obtain the final estimate:

f σ B = 456(154)(8) M σ B ( 12
) f σ B = 456(157) MeV, (62)
and: We deduce for t c (4 ∼ 8) GeV 2 at the τ -minimum:

f G1 = 365(109)(17) M σ B (1) f σ B = 365(110) MeV, (63)
M G 1 = 1563(23) tc (2) τ (138) Λ (3) G 2 (3) G 3 (6) fσ B (2) Mσ B (3) f G 1 (9) M G 1 (11) f σ B (1) M σ B = 1563(141) MeV (64)
where the error comes mainly from Λ. We have included the contributions of the σ B , σ B and G 1 .

t @GeV -2 D f G' 1 @GeVD Ú • 8 ‡ 6 5 Ï 4 L= 0.34, <aS G 2 >=6.35x10 -2 GeV 4 tc@GeV 2 D, <g 3 G 3 >=8.2 GeV 2 <aS G 2 > Figure 17: f G 1
as function of τ at N2LO from a three resonances parametrization.

• The decay constant f G 1

To extract the decay constant we shall work with the sum rules L c 2 or/and L c 3 used previously to get f G1 where L c 2 presents a τ -minimum (Fig. 17) while L c 3 an inflexion point. From L c 2 , one deduces:

f G 1 = 1000(65) tc (11) τ (150) Λ (1) G 2 (2) G 3 (15) fσ B (58) Mσ B (53) f G 1 (63) M G 1 (1) f σ B (8) M σ B (125) M G 1 = 1000(230) MeV. ( 65
)
15. The second radial excitation G 2

• The mass M G2

In so doing, we work with some judicious sum rules which optimize at smaller values of τ and are more sensitive to the high-mass meson contributions. This criterion is satisfied by R c 21 and R c 31 which one can see in Figs 18 and19 where the contributions of the σ B , σ and G 1 , G 1 have been substracted. R c 21 has a (τ, t c ) minimas at (0.18,12) in units of (GeV -2 , GeV 2 ) while R c 31 has the minimas at (0.24,10). We obtain :

M G2 = 3079(15) tc (10) τ (63) Λ (1) G 2 (2) G 3 (10) Mσ B (7) fσ B (20) M G 1 (36) f G 1 (6) M σ (1) f σ (193) M G 1 (207) f G 1 = 3079(294) MeV from R c 21 . (66)
and:

M G2 = 2879(9) tc (3) τ (64) Λ (1) G 2 (2) G 3 (10) Mσ B (8) fσ B (25) M G 1 (40) f G 1 (6) M σ (0) f σ (226) M G 1 (233) f G 1 = 2879(334) MeV from R c 31 , (67)
from which we deduce the mean value:

M G2 = 2992(221) MeV. ( 68
)
This second radial excitation is far above the ground states and first radial ecitations and may mix with the trigluonium scalar ground state [2,12]:

M 3g 3.1 GeV. ( 69
)
The decay constant f G2 We shall work with L c 3 (used to get f G1 ) for extracting the decay constant f G2 . The analysis is shown in Fig. 20 where one obtains:

f G2 = 797(27) tc (57) Λ (1) M σ B (3)M G1 (3) f G 1 (36) M G 1 (11) f G 1 (11) M G 2 = 797(74) MeV. ( 70
)
The errors due to the other terms not quoted are zero. One should note that, like in the case of f eff σ which is comparable in size, the previous coupling of the G 2 to the gluonic correlator can be interpreted as the sum of effective couplings of all higher states contributing to the spectral function.

• Truncation of the PT series

One can notice from the expression of the two-point correlator in Eq. 21 that the radiative corrections to the unit operator are huge.

M σ B and M G1

We compare explictily the effect of PT radiative corrections for the predictions of M σ B and M G1 and R c 42 in the case of "one resonance ⊕ QCD continuum" parametrization of the spectral function. For e.g t c = 4

GeV 2 , one obtains:

M σ B 768 (LO) → 1000 (NLO) → 1076 (N2LO) from R c 20 , M G1 1492 (LO) → 1522 (NLO) → 1530 (N2LO) from R c 42 , (71)
indicating that the α s corrections increase the LO masses by about 30/% for the σ B while the increase from NLO to N2LO is only about 8%. For M G1 the effect of radiative corrections is much smaller. This relatively small effect is due to the fact that radiative PT corrections tend to compensate in the ratios of sum rules.

The decay constant f σ B Fixing e.g t c 4 GeV 2 , we study the effect of the truncation of the PT series on the estimate of the coupling from L c -1 . Taking τ 0.5 GeV -2 at the inflexion point, we obtain the result in units of MeV for a "one resonance" parametrization :

f σ B = 202 (LO) → 406 (NLO) → 696 (N2LO) . ( 72
)
The effect of radiative corrections is huge as expected.

The decay constant f G1

Here we fix t c = 3.75 GeV 2 as obtained in the previous analysis . There is a large plateau in τ for different PT series. We obtain in units of MeV from L c 3 within a "one resonance ⊕ QCD continuum":

f G1 = 360 (LO) → 470 (NLO) → 580 (N2LO) , (73)
which shows relatively moderate PT corrections.

• Gluon condensates and the OPE Though the presence of the gluon condensates in the OPE are crucial for the stability points, their effects are relatively small ensuring the convergence of the OPE.

• The tachyonic gluon mass λ 2

As emphasized in [78] this effect is expected to explain the large scale of the scalar gluonium channel compared to the ordinary ρ meson channel. It is also expected [74] to give a phenomenological estimate of the uncalculated terms of the PT series.

M σ B and M G1

Including the λ 2 contribution in the OPE, we check that its effect on the mass determination is negligible which can be due to the fact that the effects tend to compensate in the ratio of moments.

The decay constants f σ B and f G1

Within a two resonances parametrization, we find for a given value of t c and τ that the tachyonic gluon mass decreases the value of f σ B by about 30-50 % while it increases only slightly f G1 by about 5%.

• Comparison with some other determinations

The σ B

The existence of the σ B with a mass about 1 GeV is compatible with the idea that it is the dilaton associated to conformal U (1) V invariance analogous to the case where the η -meson is associated to the U (1) A symmetry. The dilaton has been extensively discussed in an effective Lagrangian approach (see e.g. [46,47,89]) including the U (1) V symmetry.

-The mass in Eq. 52 is in a good agreement with the earlier sum rule result [7]:

M σ B (0.95 ∼ 1.10) GeV, (74)
from a least square fit of the USR ratio of sum rules. A similar result has been obtained from the analysis [40] based on Gaussian sum rules [41] including instanton effects.

-A full lattice calculation [48,49] finds a relative suppressed mass of the lightest glueball which moves from the quenched result about 1.6 GeV to the full QCD one around 1 GeV. However, a more approriate comparison with the lattice results requires a simulation beyond the one resonance contribution to the scalar gluonium correlator.

-A strong coupling analysis of the gluon propagator [90] finds a gluonium pole mass around M σ and M f0(980) from the mass gap [90].

The higher mass gluonium G 1 -Its mass is given in Eqs. 54 and 53 in the case of one or two resonances ⊕ QCD continuum parametrization of the spectral function. These results agree with the earlier ones in Ref. [2] from ratio of USR and can be compared to some other sum rules results [81,91,92]. We conclude fro the paper of Ref. [93] that their result for M G1 1.8 GeV is obtained outside the τ minima of their sum rule which leads to an overestimate.

-Our result for M G1 is slightly lower than the one from quenched and full lattice calculations [48][49][50][START_REF] Mathieu | [END_REF].

However, a precise comparison of the results from the two approaches is delicate as the lattice situation for the 1 GeV result remains unclear [48,49], while the lattice analysis is done by only retaining a single resonance contribution to the two-point correlator.

The gluonium radial excitations σ and G 2 and

We have also attempted to extract the masses and decay constants of the gluonium first radial excitations with the result in Eqs. 57, 58, 68 and 70. In the next phenomenological analysis, we are tempted to indentify the σ (1.32) with the f 0 (1.37). The G 2 (2.5) is (within the large errors) in the range of the observed f 0 (2.02) and f 0 (2.2) [START_REF] Zyla | [END_REF].

Comments

-In the example of the scalar gluonium channel, we have seen that (at least) NLO corrections are important in the determination of the scalar gluonium masses from the sum rules, while N2LO corrections are needed for a good determination of the decay constants. Therefore, one should make a great care for quoting the LO results given in the literature.

-The smallness of the tachyonic gluon mass contribution to the gluonium mass determinations may indicate that the PT series reaches its asymptotic form after the inclusion of the N2LO contribution which can be sufficient for extracting with a good accuracy these observables from the sum rules. However, we have also noticed that the decay constant of the σ B is largely affected (-30 to -40%) by the tachyonic gluon mass which tends to decouple the σ B from the two-point correlator.

-Some other analysis using "two resonances ⊕ QCD continuum parametrization from sum rules with instantons [40] and an inverse problem dispersive approach without instantons [94] reproduce our results

for M σ B and M G1 . The apparent discrepancies with the lattice results requests a clarification of the full QCD results and a two resonance parametrization of the two-point correlator.

-Some phenomenological implications of our results will be discussed in the next section.

17. The conformal charge ψ G (0)

• Test of the Low Energy Theorem (LET)

LET suggests that the value of ψ G (0) is given by Eq. 43. We test the accuracy of this relation by saturating the two-point function by the two lowest ground state resonances σ B and G 1 . In this way, we obtain :

ψ G (0) = i=σ,G,••• 2 f 2 G M 2 G ≈ 1.95 GeV 4 , (75)
which is slightly higher than the LET estimate in Eq.43 :

ψ LET G (0) = (1.46 ± 0.08) GeV 2 , (76)
without any appeal to some contributions beyond the OPE such as the instantons ones.

• Sum rule extraction of conformal charge ψ G (0)

Here, we shall, instead, determine ψ G (0) from the sum rule within the standard SVZ expansion. A similar analysis has been already done in [7] which we shall update here. In so doing, we work with the combination of sum rules5 : • The result of the analysis is shown in Fig. 21, where the curves present nice τ at (1.0 ± 0.04) GeV -2 .

ψ LSR G (0) = tc 0 dt t e -tτ 1 - tτ 2 1 π Imψ G (t) -L c -1 - τ 2 L c 0 QCD , (77)
A t c -minimum appears at (2 ∼ 4)GeV 2 and a stability from 8 GeV 2 . We obtain:

ψ LSR G (0) = 2.09(2) tc (0) τ (18) Λ (11) G 2 (10) G 3 (1) Mσ B ( 14
) fσ B (2)M G (7) f G (0) M σ B (2) f σ B = 2.09(29) GeV 4 , (78)
where the central value is 1.4 times the LET estimate in Eq. 43 and the one in [7] obtained within a one resonance and to LO.

One can notice that adding the two radial excitation contributions changes this value to 2.23 G 2 . The contribution of σ B is about 70% which is crucial for recovering the LET result.

• We analyze the effect of this result in the estimate of f σ B and f σ B obtained previously from the moment L c -1 where ψ G (0) contributes. We find e.g for f σ B that the inflexion point in τ moves from τ = 0.56 GeV -2 (see Fig. 13) to 0.46 for e.g t c = 3.75 GeV 2 to which corresponds to an increase of about 4 MeV, which is negligible compared to the large error of f σ B .

For f σ B , the LSR value of ψ G (0) increases the decay constant by about 100 MeV but this value is inside the range spanned by the large error of 342 MeV obtained in its determination. Therefore, we keep the result obtained in the previous section.

18. The slope ψ G (0) of the conformal charge

The slope ψ G (0) obeys the twice subtracted sum rule:

ψ G | LSR (0) = tc 0 dt t 2 e -tτ 1 π Imψ G (t) -L c -2 | QCD + τ ψ(0) G , (79)
where :

L c -2 (τ ) = β 2 (α s ) 2 π 2 τ -1 n=0,2,••• D -2 n , (80)
with:

D -2 0 = C 00 -2C 01 γ E + C 02 γ E - π 2 6 (1 -ρ 0 ) D -2 2 = τ C 21 (1 -γ E )λ 2 D -2 4 = -τ 2 40 -C 41 1 2 -γ E α s G 2 , D -2 6 = - τ 3 2 C 60 -C 61 1 3 -γ E g 3 f abc G a G b G c , D -2 8 = - τ 4 6 C 80 (1.1 ± 0.5) α s G 2 2 . ( 81
)
We show the analysis in Fig. 22 from which we deduce from t c = 2 ∼ 8 GeV 2 at the τ -minima ±0.04 GeV ψ G | LSR (0) = 0.95 (5)

tc (0) τ (2)Λ(0) G 2 (0.1) G 3 (1.5) Mσ B ( 25
) fσ B (1.1)M G1 (11.4) f G 1 (0.7) M σ (1.7) f σ (11) ψ G (0) = 0.95(30) GeV 2 , (82)
where the σ B contributes to about 80%. The present result updates the earlier one in Ref. [7].

Hadronic couplings of the σ B meson

• The σ B coupling to ππ One can use the previous result of the decay constant in Eq. 62 to predict the hadronic width of the σ B .

For this purpose, we consider the vertex (see e.g [1]):

V π (q 2 ) = π 1 |θ µ µ |π 2 : q = p 1 -p 2 with V (0) = O(m 2 π ) . (83)
In the chiral limit (m 2 π = 0), one has :

V (q 2 ) = q 2 ∞ 4m 2 π dt t(t -q 2 -i ) 1 π ImV (t) . (84)
Using the fact that V (0) = 1 [67,68], one can deduce the Low Energy Vertex sum rule (LEV-SR) [1]:

1 4 G=σ B ,••• g Gππ √ 2f G M 2 G = 1 , (85)
where the contribution of the light σ B is enhanced. Then, assuming (to a first approximation) that the σ B dominates the sum rule, one can deduce :

g σ B π + π -≈ (7.1 ± 2.5) GeV , (86)
to which corresponds the width:

Γ[σ B → π + π -+ 2π 0 ] = 3 2 |g σ π + π -| 2 16πM σ B 1 - 4m 2 π M 2 σ B 1/2 ≈ (1.4 ± 0.8) GeV. ( 87
)
To check the validity of this approximation, we use the average of the recent data from 2005 [START_REF] Zyla | [END_REF] which is :

Γ[f 0 (1.37) → ππ] 215(14) MeV =⇒ g σ B ππ (3.2 ± 0.1) GeV , (88)
while the ππ width of the f 0 (1.5) of 39 MeV [START_REF] Zyla | [END_REF] can be neglected. Including the f 0 (1.37) contribution, the previous sum rule gives:

g σ B π + π -5.7 GeV =⇒ Γ[σ B → π + π -+ 2π 0 ] 873 MeV, (89)
in agreement with previous estimates [1,2].

• The σ B coupling to K + K - Up to SU (3) breaking, the vertex sum rule in Eq.85 indicates that the σ B couples (almost) universally to Goldstone boson pairs, such that in this gluonium picture, one also expects a coupling of the σ B to K + K -. This feature has been observed from the analysis of ππ → K + K -scattering data [18][19][20].

• σ B mass and width confronted to the data It is informative to compare the previous mass and width to the one of the observed σ from scatterings data analysis :

The values of the mass and width obtained in Eqs. 52 and 87 for σ B are comparable with the onshell/Breit-Wigner mass and width obtained in Eq. 4.

These results together with the ones from data analysis reviewed in the first sections can indicate the presence of an important gluon component inside the observed σ/f 0 (500) meson wave function.

Hadronic couplings of the higher mass gluonium G

• G 1 coupling to ηη , ηη

In order to compute the g G1η1η1 coupling of the G 1 to the singlet η 1 , one starts from the [1,42]:

Ṽµν (q 1 , q 2 ) = i d 4 x 1 d 4 x 2 e i(q1x1+q2x2) 0|T Q(x 1 )Q(x 2 )θ µν (0)|0 , (90)
where θ µν is the energy-momentum tensor with 3 light quarks and :

Q(x) = α s 16π µνρσ G a µν G a ρσ , (91)
is the topological charge density, where one reminds that, in the large N c limit (no quark loops), the solution to the U (1) A problem is [1,42]:

Γ 2 (q) ≡ i d 4 x 0|T Q(x)Q(0)|0 → Γ YM 2 (q Λ) → Γ YM 2 (0) (180 MeV) 4 . (92)
Using the expression of the topological charge in the presence of quark loops [44] and the consistent one for the vertex Ṽµν (q 1 , q 2 ), one can deduce [1] the constraint:

η 1 |θ µ µ | n f =3 |η 1 9 11 12 f 2 π Γ YM 2 (0) = 2M 2 η1 = 1.15 GeV 2 , (93)
where the factor (9/11) comes from the ratio of the β-function β 1 for n f = 3 and n f = 0 and M η1 0.76

GeV [1,42]. Then, we deduce the LEV-SR [1]:

1 4 i=σ B ,G1 g Gη1η1 √ 2f i 1.15 GeV 2 . ( 94
)
Picking the singlet component of the physical η, η via the pseudoscalar mixing angle θ P -(18 ± 2) 0 [96] :

η = sin θ P η 1 + • • • , η = cos θ P η 1 + • • • , one deduces: G=σ B ,G1 g Gηη f G 4 √ 2 (sin θ P )1.15 GeV 2 -1.0 GeV 2 and g G1ηη = (sin θ P )g G1ηη . (95)
Following Ref. [1] by using the constraint from η → ηππ (allowing a 1% deviation from the precision of the data [START_REF] Zyla | [END_REF]) , and a Gell-Mann, Sharp, Wagner-type model, for the intermediate state σ B , one can deduce the upper bound for M σ B = 1.08 GeV :

|g σ B ηη | ≤ 0.34 GeV , (96)
its total width, Γ 4π 48.9 % and Γ 2(ππ) S /Γ 4π 26 %, we deduce from g f0σ B σ B and the LEV-SR in Eq. 101 the more accurate estimate :

f σ f σ B = (329 ± 30) MeV. ( 106
)
This value is in line with the intuitive expectations from ordianry ρ and ρ mesons.

Summary and conclusions

The different results obtained in this paper are summarized in Table 1 which we shall briefly comments below.

• The σ/f 0 (500) meson

In this paper, we have started to review our present knowledge on the nature of the σ meson from π + π -, γγ → π + π -, K + K -data and from J/ψ gluon rich channels (γππ, ωππ, φππ) and φ radiative decays,D (s) semileptonic decays and pp, pp production processes. The data (a priori) favour / does not exclude a large gluon component in the σ wave function where the oberved state may emerge from a maximal mixing with a qq state [4] rather than a pure four-quark ( ūu)(dd) or a ππ molecule state. The four-quark or ππ molecule state is not favoured by the too small predicted direct γγ width of the σ [17,19,20] and by the large K + K - coupling found from scattering data [18,19]. The assumption that the σ emerges from ππ rescattering is not also favoured from its large coupling to K + K -.

• The f 0 (980) meson from the data

From the scattering data analysis, it is unlikely that the f 0 (980) is an ss state because its has a non negligible / large coupling to ππ [18][19][20], while the strength of its γγ coupling is relatively small for a qq state [10,11]. Moreover, it is too light compared to the ss mass of about 1.4 GeV [11].

• Conclusions from the data

The previous observations may favour the interpretation that the σ and f 0 (980) emerge from a maximal mixing between a gluonium and a (ūu + dd) state like we have proposed in [2,4] and what had also been expected [68] in their pioneering paper.

• The σ B meson from QCD spectral sum rules (QSSR)

A QSSR analysis of the lowest ratios of sum rule L 20 favours a σ B (the subindex B refers to a pure unmixed gluonium state) having a mass M σ B =1.07 GeV (see Eq. 60) and a ππ width Eq. ( 89) compatible with the ones from scattering data for a Breit-Wigner / On-shell σ-mass.

Its decay constant can be deduced from L c 1 which is f σ B = 456(157) MeV (Eq. 62). One should also note that the set (Γ σ B , M σ B ) obtained from the sum rule is fairly compatible with the on-shell mass and width from π + π -→ π + π -, K + K -data .

• The "scalar gluonium" G 1 from QSSR An analysis of the higher weight USR L 21 and L 31 leads to a gluonium mass: M G1 = (1548 ± 121) MeV (Eq. 53) for a "3 resonances" ⊕ QCD continuum parametrization of the spectral function while its becomes (1515 ± 123) MeV (Eq.54) if one uses a "one resonance" ⊕ QCD continuum indicating that the presence of the σ B does not affect notably the value of M G1 . It is clear that higher moments are more appropriate in the QSSR approach to pick up the higher mass gluonium state due to the power mass suppression of the lowest mass σ B contribution in these sum rules.

The corresponding gluonium decay constant is f G1 = (365 ± 110) MeV (Eq. 63) which is comparable with the one of f σ B .

We have also seen that its small width into 2(ππ) S can be explained from the LET-V sum rule in Eq. 101 indicating that it is not a σ-like particle but instead decays into η η, ηη through th gluon vertex U (1) anomaly.

• The radial excitations σ , G 1 and G 2 from QSSR

We have attempted to extract the mass and coupling of the radial excitations σ of σ B , G 1 and G 2 of G 1 .

We found that M σ (1121 ± 117) MeV is rather low compared to the radial excitation of ordinary mesons where we notice that a similar feature has been observed in the analysis of the D * D molecule described by high-dimension quark currents.

f σ B 646(216) MeV can be large which we have interpreted as an effective sum of all radial excitations that the sum rule cannot disentangle. For a more phenomenoogical use, we identify the σ B with the observed f 0 (1.37) and extract the decay constant from its decay to 2(ππ) S ⊕ LEV-SR analysis from which we obtain:

f σ (329 ± 30) MeV (Eq. 101).

The mass of the radial excitation G 1 : M G 1 = 1563(141) MeV (Eq. 64) of G 1 is also found to be almost degenerated with the ground state mass M G1 = 1548(121) MeV (Eq. 61). Like in the case of σ B , we intrepret the large coupling f G 1 = 1000(230) MeV (Eq. 65) as a sum of the effective couplings of all higher radial excitations entering into the spectral function.

The 2nd radial excitation gluonium G 2 is found to have a relatively high mass M G2 2.99 GeV.

• Phenomenology

The global picture which emerges from the approach indicates that the observed gluonia candiates can de classified into two groups :

-The σ-like gluonia [σ/f 0 (0.5), f 0 (1.37)] which emerge from (σ B , σ B ) and which decay into ππ via OZI violating process while the f 0 (1.37) decays to 2(ππ) S via virutal σσ.

-The G-like gluonia [f 0 (1.5), f 0 (1.7)] which come from (G 1 , G 1 ) and which decay into η η and ηη through the gluon U (1) anomaly vertex but not copiously to ππ and 2(ππ) S .

-We may expect that high-mass 2nd radial excitation G 2 (2.99) has similar properties and brings some gluon component to the gluonia cnadidates above 2 GeV which contributes to their η η and ηη decays.

We remind that the masses of the S 2 (ūu + dd), S 3 (ss) quarkonia states and their radial excitations are predicted from the sum rules to be (in units of GeV) [2] :

M S2 1 , M S 2 (1.1 ∼ 1.4), M S3 (1.47) , M S 3 (1.7 ∼ 2.4), (107)
which can mix with the previous gluonia states.

A more complete analysis of the data can be done within some qq meson-gluonium mixing scheme. The case for the σ(500)-f 0 (980) has been studied in [4] but may/should be extended to the higher meson masses. Some attempts in this direction have been also done in [2] 6 but needs to be updated.

We have not updated our analysis of the γγ and radiative widths obtained in [1,2]. We expect that these results are still valid within the errors.

We emphasize that a more appropriate comparison of the results obtained in this work (e.g with lattice calculations) requests an analysis beyond the "one resonance" contribution to the two-point correlator from these alternative approaches. Indeed, a simple one resonance for parametrizing the spectral function is (obviously) insufficient for describing the complex and overpopulated spectra of the I = o scalar mesons.

• The conformal charge ψ G (0) from QSSR

We have completed our analysis by the estimate of the conformal charge ψ G (0) and its slope ψ G (0).

Using the previous values of the σ and gluonia G masses and couplings, we have obtained ψ LSR G (0)=2.09 (29) GeV 4 (Eq. 78) which can be compared with the LET prediction of about 1.46 GeV 4 (Eq.43).

The result indicates that the σ B meson contributes predominantly in the spectral function for reproducing the LET prediction of ψ G (0). It also shows that the instanton effect in the estimate of ψ G (0) is (a priori) not necessary for a correct estimate of the conformal charge ψ G (0). A similar feature has been observed in the estimate of the topological charge and its slope in the U (1) A sector [13] .

The value of the slope is ψ G (0)| LSR = 0.95(30) GeV 2 .

• Prospects

We plan to study in more details the mixing of the qq states with the higher gluonia masses σ (1.15), G 1 (1.53) and eventually with the G 2 (2.48).

We plan to extend the analysis to the case of (pseudo)scalar C = ± trigluonia channels [98] which has been studied in [12] and [99].

Figure 1 :

 1 Figure 1: L c -1 as a function of τ at N2LO for different values of tc.

Figure 2 : 2 >Figure 3 : 3 L= 2 >Figure 4 : 3 L=Figure 5 : 2 >Figure 6 : 2 >Figure 7 :

 223324352627 Figure 2: L c 0 as a function of τ at N2LO for different values of tc.

<aS G 2 Figure 8 :

 28 Figure 8: Mσ B from R c 20 as a function of τ at N2LO for different values of tc including the effect of the 2nd resonance with (M G 1 , f G 1 ) = (1.50, 0.39) GeV [2].

  the τ behaviour of R 42 for different values of t c .

. 34 ,Figure 9 :

 349 Figure 9: M G 1 as a function of τ at N2LO for different values of tc from a two resonances parametrization.

Figure 10 : 3 <aS G 2 >=6.35x10 - 2 Figure 11 :

 1032211 Figure 10: f G 1 from L c 2 as a function of τ at N2LO for different values of tc from a two resonances parametrization and for given values of (Mσ B , (fσ B ).

4 <g 3 G 3 >=0. 8

 4338 GeV 2 <aS G 2 >, fs B =561 MeV

Figure 12 : 3 G 3 >=8. 2 Figure 13 :

 1233213 Figure 12: f G 1 from L c 2 and L c 3 from the τ -minimum τ 0.5 GeV -2 from previous figures as a function of tc for given values of (Mσ B , fσ B ).

<g 3 G 3 >=8. 2 2 >Figure 14 :

 332214 Figure 14: M σ B as a function of τ at N2LO for different values of tc within a three resonances parametrization. from R c 20

Figure 15 :

 15 Figure 15: f σ B as a function of τ at N2LO for different values of tc within a three resonances parametrization. from L c -1

14 . 1 • The G 1 mass<g 3 G 3 >=8. 2 2 >Figure 16 :

 1411332216 Figure 16: M G 1 as a function of τ at N2LO for from a four resonances parametrization.

9 L=Figure 18 : 8 L= 0 . 34 ,Figure 19 :

 918803419 Figure 18: M G 2 as a function of τ at N2LO for a four resonances parametrization. from R c 21 .

15 L= 0 . 34 , 2 4 tc@GeV 2 D, <g 3 G 3 >=8. 2 Figure 20 :

 150342433220 Figure 20: f G 2 as a function of τ at for different values of tc for a four resonances parametrization from L c 3

2 L= 0 . 3 G 3 >=8. 2 Figure 21 :

 2033221 Figure 21: ψ(0) G from LSR as function of τ for different values of tc.

- 2 : 2 L= 0 .Figure 22 :

 22022 Figure 22: ψ (0) G from LSR as a function of τ for different values of tc.

An estimate of this mixing angle using Gaussian sum rules favours a maximal mixing[40].

For reviews, see e.g.[72,73].

The name inverse Laplace transform has been attributed due to the fact that perturbative radiative corrections have this property.

We remind that this combination of sum rules has been firstly (and successfully) introduced to quantify the deviation from kaon PCAC[95] due to SU (3) breakings and later to estimate the topological charge of the U (1) A channel[7] and its slope[7,13] within the standard SVZ expansion without any additional instanton contributions.

For a review prior 1996, see e.g[97].

which is a relatively small contribution. Neglecting this contribution and the f 0 (1.37) one, the LEV-SR in Eq. 95 leads to : g G1ηη 2.74 (78) GeV and g G1ηη 0.85 [START_REF] Zyla | [END_REF] GeV ,

We deduce the width for M G1 = 1.55 GeV, p ηη = 20 MeV and p ηη = 517 MeV :

G1

(2.5 ± 1.4) [(2.6 ± 0.9) Data] MeV, (98) which implies:

The results agree perfectly within the errors with the data.

One can also write a low-energy vertex sum rule (LEV-SR) for this coupling [1]:

Using the dispersion relation of the vertex, one can write the sum rule:

where here it is more appropriate to use the pole mass M σ 0.5 GeV for the virtual σ B .

Assuming that the σ dominates the sum rule which is also the only one which decays copiously into 2(ππ) S [START_REF] Zyla | [END_REF], one can deduce :

where we have used the previous estimate of the coupling in Eq. 58.

We confront this prediction with the data of f 0 (1.37.1, 5) → 2(ππ). In so doing, we use the six most recent measurements of the f 0 (1.37) 4π width since 1995 compiled by PDG [START_REF] Zyla | [END_REF] from which we deduce the mean:

Using the ratio [START_REF] Zyla | [END_REF] :

we deduce:

where the estimate is off by about a factor two indicating that the sum rule may have oversetimated the decay constant of the σ -meson which is (presumably) the sum of effective couplings of higher state resonances.

Declaration of interests

☒ The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.


Paolo Gianessi

Audrey Cerqueus

Damien Lamy

Xavier Delorme

Using Reconfigurable Manufacturing Systems to minimize energy cost: a two-phase algorithm

Keywords: Energy-efficiency, Flexible and reconfigurable manufacturing systems, Process Planning/Equipment Selection

Energy consumption has become a major concern for society, and since the industrial sector is the largest consumer, companies are urged to improve energy efficiency of their production systems. This paper investigates how Reconfigurable Manufacturing Systems (RMS), and particularly their scalability feature, can be exploited in order to minimize the energy cost in production systems w.r.t. a Time-Of-Use pricing scheme. In the case of RMS, the resulting energy cost optimization problem is a Bilevel Optimization problem, as it jointly considers both the line balancing and the production planning. After introducing the problem and its features, a solving approach based on a simulated annealing algorithm and a linear program is proposed. The approach is then validated on designed instances based on classical test problems taken from the literature. Results show that considerable savings in terms of energy cost can be achieved w.r.t. dedicated lines, even when optimally designed, ultimately showing the great potential of RMS towards energy efficiency.

INTRODUCTION

Energy consumption worldwide is supposed to rise by nearly 50% by 2050 and the industrial sector is the largest energy consumer, with an electricity usage that is expected to grow by 40% in the projected period (U.S. [START_REF] Eia | International energy outlook 2019[END_REF]. However, the concern of society for the environment is driving companies to adopt more energy-efficient and sustainable production systems. This also implies the design and management of energy-efficient manufacturing systems (MS) to increasingly consider renewable energy sources, which usage is expected to spread in all sectors during the next decades [START_REF] Battaïa | Sustainable and energy efficient reconfigurable manufacturing systems[END_REF]. Consequently, energy efficiency must be considered at all levels, be it strategic, tactical or operational [START_REF] Giret | Sustainability in manufacturing operations scheduling: A state of the art review[END_REF]. When considering energy efficiency of production systems, three measures are usually referred to (see e.g. [START_REF] Masmoudi | Job-shop scheduling problem with energy consideration[END_REF]): total energy consumption; total energy cost w.r.t. to a given pricing scheme -more often Time-Of-Use (TOU); power peak limitation. Moreover, world gross output is expected to double worldwide in the next decades, with high volatility of the market which requires production systems to quickly adapt their throughput. The notion of Reconfigurable Manufacturing System (RMS), introduced by [START_REF] Koren | Reconfigurable manufacturing systems[END_REF], aims to achieve such reactivity by means of efficiently reconfigurable production systems. Typical RMS are composed by several serial workstations with multiple parallel identical machines used at each stage. The machines on each workstation are generally computer numerical control machines or reconfigurable machine tools, but can also consist in other types of resources (e.g., workers with cobot). RMS can be an interesting lever to deal with variable energy availability or pricings, as it is the case with TOU, since it allows to change the system structure to better fit energy requirements. In costlier periods, a less consuming configuration, even with lower productivity, can be used, before switching to a higher-throughput higher-consumption configuration in periods with lower energy prices and meet the demand. In this paper, we deal with the problem of balancing an RMS in order to define a most fitting set of configurations so as to minimize the energy cost associated with a given demand and w.r.t. a energy cost profile. The problem has been shown to be a Bilevel Optimization problem in Cerqueus et al. (2020b) which first studied it. The authors have proposed a two-phase decomposition approach to tackle it but have only shown its effectiveness on a case study with realistic energy cost parameters. In this work, we extend the two-phase approach and perform an extensive computational experience so as to validate it and better assess the potential gains in terms of energy cost that could derive by the use of RMS. The plan of the paper is as follows. Section 2 reviews similar works in the literature. Section 3 outlines the problem at study. Section 4 delves into the details of the decomposition approach. Section 5 describes the computational sessions. Finally, Section 6 shows the conclusions and perspectives of this work.

LITERATURE REVIEW

RMS have been introduced in [START_REF] Koren | Reconfigurable manufacturing systems[END_REF]. They aim at reaching as much flexibility as flexible MS while keeping as high a production rate as dedicated lines. Specific features of RMS are modularity, integrability, convertibility, diagnosability, customization and scalability [START_REF] Koren | Value creation through design for scalability of reconfigurable manufacturing systems[END_REF]. Scalability is enabled by the capacity to integrate or remove parallel machines on workstations and motivated by the need to adapt the production level to fluctuations of demand or energy price. The literature on RMS mainly deals with system design, layout, process planning, and reconfigurable control [START_REF] Bortolini | Reconfigurable manufacturing systems: Literature review and research trend[END_REF], and it seems that few papers consider scalability as a lever for productivity. [START_REF] Deif | Assessing capacity scalability policies in RMS using system dynamics[END_REF] investigates a model for assessing the scalability capacity of a make-to-order RMS according to different demand scenarios and performance measures. In [START_REF] Wang | Scalability planning for reconfigurable manufacturing systems[END_REF] a scalability planning methodology for reconfigurable manufacturing is explored. The approach consists in changing the capacity of an existing system by successive reconfigurations, in order to minimize the number of machines required to respect a new throughput. In [START_REF] Hees | A production planning system to continuously integrate the characteristics of reconfigurable manufacturing systems[END_REF], a production planning system that focuses on the scalability of the production capacities and the adaptation of functionalities is investigated. In [START_REF] Moghaddam | Configuration design of scalable reconfigurable manufacturing systems for part family[END_REF], two different approaches are considered to design multi-product and scalable RMS for multiple production periods. The objective is to minimize design and reconfiguration costs while fulfilling a given demand. Both approaches are based on the estimated demand of different parts: (i) up-or downgrade the RMS depending on estimates for each production period, or (ii) select and reconfigure the reconfigurable manufacturing tools for all periods based on longer-term estimates. [START_REF] Putnik | Scalability in manufacturing systems design and operation: State-of-the-art and future developments roadmap[END_REF] states that scalability can improve the optimization of manufacturing system design and management, and help to develop new paradigms for sustainability, but few papers actually consider this last criterion. [START_REF] Zhang | Modeling and verification of reconfigurable and energy-efficient manufacturing systems[END_REF] introduces the concept of energyefficient RMS and investigates a discrete event simulation model to evaluate the systems energy efficiency. For RMS, [START_REF] Choi | A holistic production planning approach in a reconfigurable manufacturing system with energy consumption and environmental effects[END_REF] investigates a multi-objective production planning problem that considers energy consumption, throughput, and inventory holding costs to assess the performance of the planning. A configuration corresponds to a production plan which is adjoined by a total energy consumption. In [START_REF] Touzout | Multi-objective sustainable process plan generation in a reconfigurable manufacturing environment: exact and adapted evolutionary approaches[END_REF] a multi-objective problem for sustainable process planning in RMS is addressed. Three optimization criteria are considered: (i) the total production cost, (ii) the total completion time, (iii) the amount of greenhouse gases emitted. [START_REF] Gianessi | Simple assembly line balancing problem with power peak minimization[END_REF] is one of the first papers dealing with energy at the design stage for dedicated lines, more precisely by minimizing peak power. The recent survey of [START_REF] Battaïa | Sustainable and energy efficient reconfigurable manufacturing systems[END_REF] shows that few research projects exploit the capacity of RMS to improve energy efficiency and sustainability in production. To the best of our knowledge, the first study that has considered scalability as a lever to adapt energy consumption in the context of RMS is the one of Cerqueus et al. (2020b).

PROBLEM DESCRIPTION

In this work, we consider the optimization problem of determining both a set of configurations of an RMS at the design stage, and the planning of their use. The set of configurations must be large enough so that the RMS can adapt the production planning according to the varying energy costs of a TOU pricing scheme, while fulfilling a demand ∆ over a given period T . The planning then consists in selecting which configuration should be used at each moment of the timespan T : the objective is to minimize the total energy cost of production. This problem, which has first been studied in Cerqueus et al. (2020b), is actually a Bilevel Optimization problem (see e.g. [START_REF] Colson | Bilevel programming: A survey[END_REF], or [START_REF] Moore | Extensions to the multilevel programming problem[END_REF] for an application to production planning). Indeed balancing and planning problems are usually considered separately by manufacturers with different decision-makers and time horizons -the former is a strategic problem, the latter a tactical/operational issue. In the case of RMS, however, it is useful to design the production system from the start in such a way that its scalability can be exploited as much as possible w.r.t. to a given objective -here energy cost. As a consequence, and regardless of the specific objective, the two decision layers become strongly interconnected: one can easily see that on one hand the performance of a configuration set cannot be assessed without constructing the planning of how configurations are used over the period T , while on the other hand the planning problem cannot be solved if the configuration set is not known in advance. Moreover, the design of an RMS is more complex than that of a dedicated MS. The design decisions of a dedicated MS mainly concern how the n production tasks are assigned to the m ≤ n available workstations, so as to comply with technological constraints and optimize production criteria. Even the Simple Assembly Line Balancing Problem (SALBP), that studies a straight paced assembly line in which design decisions only concern the takt time and/or the number of workstations, is NP-Hard and has given rise to a large amount of scientific works (see e.g. [START_REF] Scholl | State-of-the-art exact and heuristic solution procedures for simple assembly line balancing[END_REF]. The consideration of a different criterion can only make the problem even harder, as shown for instance in [START_REF] Gianessi | Simple assembly line balancing problem with power peak minimization[END_REF] for power peak. With RMS, the task is much harder, as the goal of the design stage is to determine a set of line configurations that differ in the assignment of tasks to workstations and in the number of identical resources assigned to each workstation. In this work, as it has been done in Cerqueus et al. (2020b), we focus on scalability and consider the particular case in which all the configurations of the set share the same balancing, and only differ in the number r k of resources of each workstation k. As a direct consequence of this assumption, system reconfigurations can be seen as switching on/off resources, and hence reconfiguration times can be neglected since a reconfiguration does not imply a reassignment of tasks and materials to workstations. Under this assumption, and by considering the minimization of the energy cost of production as the main objective, the goal of the design stage is to determine the balancing x that gives rise to the most fitting set of descending configurations C(x) in terms of per-time-unit energy consumption Q i and takt time c i , i ∈ C(x). The Bilevel Optimization problem at hand can then be expressed as:

min x F(x, y ); max x S(x) (1) s.t. B(x) ≤ 0 (2) y = argmin y E(x, y) (3)
s.t. P(y) ≤ 0 (4)

In the upper-level problem (1)-( 3), F is an overall cost function and S a scalability measure of x;

(2) enforce the balancing constraints;

(3) certifies that planning y has minimum energy cost E, given balancing x.

Binary upper-level variables x jk take value 1 if task j ∈ {1,...,n} is assigned to workstation k ∈ {1,...,m} w.r.t. an upper bound n max on the number of tasks assigned to the same workstation. Each task features a processing time t j and a total energy consumption e j (i.e. the integral of the power consumption profile of task j over its processing time). The basic configuration that descends from balancing x has r k = 1 ∀k ∈ {1,...,m}, and workstation k could have takt time c k = W k r k equal to its workload W k = j t j x jk . The takt time of the line, c, is then set according to the bottleneck workstation k = arg max c k . Other configurations can be derived from balancing x by increasing the number r k of resources of workstations, but only those having a higher production rate, i.e. a lower c, are worth being considered. By definition, c can only be reduced by increasing by one the number r k of resources of the bottleneck station. This can result in a possible new bottleneck station. This process can be repeated until the bottleneck station would exceed the maximum number of resources per workstation r max . By following this reasoning, one can obtain the full set C(x) of configurations originating from x, with descending takt time values. As for energy consumption, we assume that during idle time I k = c -c k stations have a residual energy consumption which is proportional by a factor α to I k , the number of resources r k and the per-resource, per-time-unit energy consumption of tasks assigned η k W k , with η k = j e j x jk . Hence, the energy consumption of workstation k during a takt is

E k = η k (1 + α r k I k W k )
. This allows to associate each configuration i ∈ C(x) with a per-time-unit energy consumption Q i :

Q i = 1 c i k E i k = 1 c i k η k (1 + α r i k I i k W k ) (5)
in which η k and W k are invariants of the configuration i, while c i , r i k , I i k and thus E i k depend on i. The lower-level problem (3)-( 4) represents the planning decisions, i.e. how configurations of C(x) should be used w.r.t. a given pricing scheme to minimize the production energy cost while ensuring fulfilment of demand ∆. Since we are considering a TOU pricing scheme, the time horizon T is divided in several cost periods and lower-level variables y ip represent the portion of period p spent producing according to configuration i ∈ C(x). Constraints (4) enforce planning feasibility and demand fulfilment, while term E(x, y) of ( 3) is the total energy cost of producing according to configurations C(x) and planning y.

DECOMPOSITION APPROACH

As in Cerqueus et al. (2020b), the algorithmic approach presented in this paper for the Bilevel Optimization problem is a two stage decomposition approach. Phase 1 addresses the balancing problem through a metaheuristic and yields a balancing x, and a set C(x) of configurations generated from it, which seek to optimize both the scalability of the production system and its capacity to fulfill a large range of economical demands. Starting from set C(x), Phase 2 solves the planning problem to minimize the energy cost of the production process needed to satisfy demand ∆ over time horizon T . In the following, we describe the two phases in detail.

Phase 1: Generating a Configuration Set

Phase 1 of the metaheuristic consists in a Simulated Annealing (SA) algorithm (see e.g. [START_REF] Aarts | Simulated annealing[END_REF]. In the two-phase approach for the problem at hand, the SA-based algorithm begins by randomly building a feasible balancing x, i.e. that complies with the other previously mentioned constraints, as well as with precedence constraints (here denoted by symbol ≺) among tasks:

(∀ j 1 , j 2 ∈ {1,...,n} : j 1 ≺ j 2 ) k k(x j2,k -x j1,k ) ≥ 0
This also defines the basic configuration, i.e. with r k = 1 ∀k, of set C(x). Other configurations i ∈ C(x) are then derived from it as previously described in Section 3. The neighbors of a balancing x are those which can be obtained by moving one task j from its current workstation

k to another k = k, k 1 ≤ k ≤ k 2 , where (∀ j 1 ≺ j) k k x j1,k ≤ k 1 ; (∀ j ≺ j 2 ) k 2 ≤ k k x j2,k
i.e. such that no precedence constraint involving j is violated. Moreover, a candidate destination workstation cannot be chosen if n max tasks have already been assigned to it. New workstations can be opened to be assigned tasks, up to a maximum number m. A balancing x is evaluated by means of a weighted sum Φ(x) of two terms, namely a hypervolume metric H(x) and the production rate P(x) of the most productive among the configurations derived from x. Hypervolume H(x) is computed based on the takt time c i and per-time-unit energy consumption Q i of configurations i ∈ C(x) that are non-dominated configurations w.r.t. c and Q (i.e. those for which no other configuration i ∈ C(x) exists s.t. c i ≤ c i and Q i ≤ Q i ). We suppose without loss of generality that configurations i ∈ C(x) are sorted by decreasing c values. The values of c i and Q i are normalized using suitable upper (c U , Q U ) and lower (c

L , Q L ) bounds: c = c-c L c U -c L and Q = Q-Q L Q U -Q L .
The hypervolume of balancing x is then:

H(x) = (1 -c1 )(1 -Q1 ) + i≥2 (c i-1 -ci )•(1 -Qi ) (6)
To have the same order of magnitude, the second part of the weighted sum is also normalized:

P(x) = c U c L c U -c L ( 1 c |C(x)| - 1 c U ) (7)
The fitness Φ(x) of x (with weight factor λ) is then:

Φ(x) = λ H(x) + (1-λ) P(x) (8)
Hypervolume H allows to use an aggregated criterion to evaluate the scalability of a set C(x) (Cerqueus et al., 2020a). The maximization of H is meant to hopefully lead to configuration sets whose c and Q values are at the same time as diversified and as low as possible. The latter aspect, in particular, is expected to yield configurations which are both highly producing and energy efficient. Figure 1 depicts the calculation of H for an example set.

The second term of fitness function Φ in ( 8) is proportional to P to further stress the fact to seek for a configuration with a consistent highest production level.

0 c Q H c1 c|C(x)| Q1 Q|C(x)| × × × × × × × × × 1 1
Fig. 1. Hypervolume calculation for a set of 8 configurations (1 dominated, in red).

Phase 2: Planning Configurations over the Timespan

The planning problem of Phase 2 consists in optimally deploying the configurations of C(x) over a period T to fulfill demand ∆ and is solved as a Linear Program (LP).

We suppose that the considered pricing scheme is of type TOU, i.e. with a set P of time periods each featuring a energy unit cost U p ; durations of periods D p are s.t. T = p∈P D p . Lower-level variables y ip are then [0, 1] real variables modeling the portion of time period p spent producing with configuration i ∈ C(x). The LP model is: min

i∈C(x),p∈P D p • U p • Q i • y ip (9) s.t. i∈C(x),p∈P Dp c i • y ip ≥ ∆ (10) i∈C(x) y ip ≤ 1 ∀p ∈ P (11)
D p y ip is the time spent using configuration i ∈ C(x) to produce during period p. Hence:

• (9) represents the energy cost to minimize, since Q i D p y ip is the associated energy consumption, and • (10) ensures that demand ∆ is satisfied, since 1 c i D p y ip is the corresponding produced quantity, resulting from the production rate 1 c i . Finally, (11) impose that production in a period p cannot last more than D p , i y ip < 1 meaning that production is suspended for a time (1i y ip )D p .

COMPUTATIONAL EXPERIENCE

We tested our approach on benchmark SALBP instances derived from [START_REF] Scholl | Data of assembly line balancing problems[END_REF] and [START_REF] Otto | Systematic data generation and test design for solution algorithms on the example of SALBPGen for assembly line balancing[END_REF]. For the instances of [START_REF] Scholl | Data of assembly line balancing problems[END_REF], since the author provides sets of instances with same tasks and different c values, we derived one instance per set, with median takt time value. The obtained SALBP instances are enhanced with energy, demand, time horizon and TOU pricing features. The demand is directly derived from the takt time of the SALBP instance. The energy consumption e j of each task j is computed as the product of its processing time t j and a randomly generated average power consumption (uniform distribution in [5;50]). We considered a energy cost profile in which a time horizon T of 24h (hours) is subdivided in 6 periods: 1 low-cost period of 8h with a (energy) unit cost of 18, 3 medium-cost periods of total 10h with unit cost 65, and 2 high-cost periods of total 6h with unit cost 108. In order to have (for comparison purposes) a reference dedicated line we run a SALBP-1 model on IBM CPLEX with a 3 hours time limit. We discarded Scholl instances for which no SALBP-1 feasible solution could be found within time limit (namely two instances, called Scholl and Barthol2). We only considered 36 instances of [START_REF] Otto | Systematic data generation and test design for solution algorithms on the example of SALBPGen for assembly line balancing[END_REF], 18 of which are taken from the set of instances with n = 20 tasks and 18 from the set with n = 50 tasks, evenly distributed in order to have a representative sample of each of the parameter values that the authors of [START_REF] Otto | Systematic data generation and test design for solution algorithms on the example of SALBPGen for assembly line balancing[END_REF] used for the generation of the instances. All of such Otto instances have a SALBP-1 feasible solution within time limit and are thus considered here. The determination of the reference dedicated line also allows to complete the instances with the upper limit on the number of stations, which is fixed to one third higher than the number of stations obtained when solving the SALBP-1 version. As for the simulated annealing, the descent factor for the temperature is 0.98 and the initial temperature is set to 10. The weight factor in the fitness function Φ (see ( 8)) is λ = 0.5. The total number of iterations and the length of the steps are proportional to the number of tasks in the instance (resp. 1250n and 5n). The simulated annealing having a random component, for each instance the whole of SA and LP is replicated 10 times. Table 1 compares our approach with the performance in terms of energy cost saving that could be achieved w.r.t. the reference dedicated line. Two levels of demands are considered, "lower demand" is the demand directly derived from the takt time of the instances and "higher demand" is a demand 25% higher. Phases 1 and 2 of the metaheuristic use the same set of configurations for both lower and higher demands, as demand is not a parameter of the metaheuristic (contrary to SALBP-1 which considers the demand as known). The value N/A for the higher demand for Bowman denotes that the demand is higher than what could be achieve by the reference dedicated line (the desired takt time is smaller than the higher task processing time). However, our approach finds a solution, since it can use parallel resources on the stations. Instances with name marked by a '*' are those for which the solution of the SALBP-1 with higher demand is not optimal (due to exceeded time limit). For lower demands, using our planning instead of a dedicated line results in an overall 65% average reduction in cost. This does not seem to depend on the dataset, nor on the instance size (69% for Otto with n = 20, 64% for Otto with n = 50, 63% for Scholl). For higher demands, the cost reduction is lower (56% on average). For all instances but six, the reduction of cost by our approach is significantly lower for higher demands than for lower demands. The difference between categories and size of instances is more visible than for lower demands (64% for Otto with n = 20, 54% for Otto with n = 50, 51% for Scholl). As for the occupation rates of the different time periods of the TOU scheme, without surprise the dedicated lines use almost all of the time horizon T . Only high-cost periods are used less than 100% (95% on average for lower demands, 97% for higher demands). For the RMS, low-cost periods are densely used (96% on average for lower demands, almost 100% for higher demands), medium-cost periods are scarcely used (average 3% for lower demands, 14% for Table 1. Energy cost of the planning of an RMS yielded by our method (SA and LP), compared with that of a dedicated line. Term t SA is the solving time of our approach. Columns 4, 5 and 7 are resp. the number of configurations output by the SA and the average number (over 10 replications) of those used during the planning to satisfy lower and higher demands. Cost columns are energy cost percentage decrease values, i.e. the energy cost of the RMS, minus that of the dedicated line, divided by the latter.

higher demands), high-cost periods are never used. The average number of used configurations is always below 2. The typical planning generated for the RMS has the following structure. Let j, k ∈ C(x) be the configurations with lowest per-takt energy consumption and highest productivity, respectively. Then:

• if the demand can entirely be fulfilled during the low cost period, then the configuration j and k are used so that the per-produced-item energy consumption is minimized and the demand is answered; • if the demand is greater than such limit, then configuration k is used for the full low-cost period and the exceeding demand is fulfilled in priority during the medium-cost periods, according to the same scheme.

This planning profile is a consequence of the fact that there is only a 5% difference in terms of per-produceditem energy consumption between the most and the least consuming configurations, which is largely smaller than the difference of per-time-unit energy cost between low-, medium-and high-cost periods (40%). Few replications (8 out of 610) use a configuration that is neither the most productive, nor that with the lowest pertakt energy consumption, but instead offers an interesting tradeoff between productivity and energy cost. Since this solution structure is recurring, simple heuristics could be derived to easily build the production planning of an RMS and achieve significant cost economies.

CONCLUSION

In this article, we address the Bilevel Optimization problem of balancing a Reconfigurable Manufacturing System (RMS) and derive a set of system configurations that would allow to obtain a production plan of minimum energy cost w.r.t. a Time-Of-Use (TOU) pricing scheme and capable of meeting a given demand over a given time horizon. The purpose of this study is to show the relevance of RMS when dealing with energy-efficient production, an issue nowadays more and more sensitive.

The problem has been first studied in Cerqueus et al. (2020a), which investigated whether the scalability of RMS can be a lever for energy efficiency. We deal here with the same problem and propose an extended version of the two-phase metaheuristic. A thorough computational experience is performed to better assess both the problem interest and the approach effectiveness. In Phase 1, a Simulated Annealing seeks for a balancing of the production line. A balancing is evaluated by assessing the quality of the derived configuration set via a weighted sum of the highest attainable production rate and of the scalability of the production system (defined as the diversity and the performance of its configuration set). Phase 2 solves a Linear Program to optimally use the configurations derived from such balancing over the considered time horizon.

Tests have been run on benchmark instances of the SALBP, enhanced with problem-specific features, and results have been compared with dedicated lines composed of optimal or near-optimal SALBP solutions. The comparison shows that substantial savings can be achieved in terms of energy costs, even with tight levels of demand, by using a low number of configurations on average. Many future research paths can be outlined. First, Phase 1 could benefit from the use of other neighborhood operators or metaheuristic schemes, e.g. Variable Neighborhood Search. Moreover, the proposed metaheuristic can be improved by a feedback mechanism so as to include in the evaluation of a balancing the optimal cost of a planning made up of the configurations derived from it. A feedback mechanism would also allow to better exploit the Bilevel nature of the problem, as opposed to the two-steps hierarchization considered here, and ultimately lead to algorithms capable of achieving globally better solutions.

To this end, the most promising direction seems the design of an exact Bilevel Optimization algorithm. The multiobjective nature of the upper-level problem also deserves further investigations. The structure of solutions could be studied from a theoretical point of view to seek for properties (e.g. in terms of dominance w.r.t. takt time and energy parameters) that could help design even more performing algorithms to obtain energy-efficient RMS. By removing the same-balancing assumption of Cerqueus et al. (2020a), a far wider range of industrial cases could be considered, as well as potentially better (since less constrained) configuration sets. In this case, reconfiguration times could not be neglected anymore in the planning problem but more energy-efficient configurations could be available and ultimately allow larger economic savings.

Lastly, future research works could also consider some other industrial constraints, such as a power peak limit.