Georges Nassopoulos

Patricia Serrano-Alvarado

Pascal Molli

Emmanuel Desmontils

Extracting Basic Graph Patterns from Triple Pattern Fragment Logs

Keywords: Linked Data, Triple Pattern Fragments, log analysis, Basic Graph Patterns

The Triple Pattern Fragment (TPF) approach is de-facto a new way to publish Linked Data at low cost and with high server availability. However, data providers hosting TPF servers are not able to analyze the SPARQL queries they execute because they only receive and evaluate queries with one triple pattern. In this paper, we propose LIFT: an algorithm to extract Basic Graph Patterns (BGPs) of executed queries from TPF server logs. Experiments show that LIFT extracts BGPs with good precision and good recall generating limited noise.

Introduction

The Triple Pattern Fragment (TPF) approach is de-facto a new way to publish Linked Data at low cost and with high availability for data providers [15]. WarDrobe [1] provides more than 38 billions of triples distributed over 65 datasets. Following the TPF approach, most of the SPARQL query processing is now executed on the client-side, TPF servers only receive and evaluate queries with one triple pattern. Consequently, data providers of TPF servers do not know the executed SPARQL queries and cannot analyze them as data providers do with queries of SPARQL endpoints.

Knowing executed SPARQL queries is fundamental for data providers. Mining logs of SPARQL endpoints allows to detect recurrent patterns in queries for prefetching [START_REF] Lorey | Detecting SPARQL Query Templates for Data Prefetching[END_REF] or for benchmarking [8]. It provides the type of queries issued, the complexity and the used resources/predicates [START_REF] Möller | Learning from Linked Open Data Usage: Patterns & Metrics[END_REF][START_REF] Picalausa | What are Real SPARQL Queries Like?[END_REF]. It allows also to distinguish between man or machine made queries [11,12]. Currently, such analysis cannot be done on logs of TPF servers because they only contain information about single triple patterns.

In this paper, we propose LIFT (LInked data Fragment Tracking): an algorithm to extract Basic Graph Patterns (BGPs) from logs of TPF servers. Compared to the state of art, [START_REF] Verborgh | Initial Usage Analysis of DBpedia's Triple Pattern Fragments[END_REF] reported statistics from the logs of the DBpedia's TPF server. In previous work [9], we proposed an algorithm to extract BGPs of federated SPARQL queries from logs of a federation of SPARQL endpoints.

Here, we address a similar scientific problem but in the context of a single TPF server.

The main challenge to extract BGPs is the concurrent execution of SPARQL queries on one TPF server. If we find a function f , to extract BGPs from isolated traces of one SPARQL query, is f able to extract the same BGP from traces of concurrent SPARQL queries? LIFT faces this problem by tracking the bindings among different triple pattern queries to detect joins. We experimented LIFT with different levels of concurrency. We demonstrate in which conditions, it extracts BGPs with good precision and good recall generating limited noise. Thanks to LIFT, we were able to extract the frequent BGPs from the TPF log published in the USEWOD 2016 dataset [5].

Next section introduces a motivating example and our problem statement. Section 3 presents LIFT and Section 4 shows our experiments. Section 5 presents related work. Finally, conclusions and future work are outlined in Section 6.

Motivating example and problem statement

In Figure 1, two clients, c 1 and c 2 , execute concurrently queries Q 1 and Q 2 over the DBpedia's TPF server. Q 1 asks for movies starring Brad Pitt and Q 2 for movies starring Natalie Portman. 1 Both queries have one BGP composed of several triple patterns (tp n ).

SELECT ?movie ?title ?name WHERE { ?movie dbpedia-owl:starring ?actor .

(tp1) ?actor rdfs:label "Brad Pitt"@en .

(tp2) ?movie rdfs:label ?title .

(tp3) ?movie dbpedia-owl:director ?director .

(tp4) ?director rdfs:label ?name (tp5) FILTER LANGMATCHES(LANG(?title), "EN") FILTER LANGMATCHES(LANG(?name), "EN") } c1(173.28. 19.114) : Query Q1 SELECT ?titleEng ?title WHERE { ?movie dbpprop : starring ?actor .

(tp 1 ) ?actor rdf s : label "N atalie P ortman"@en .

(tp 2 ) ?movie rdf s : label ?titleEng .

(tp 3 ) ?movie rdf s : label ?title (tp 4 ) FILTER LANGMATCHES(LANG(?titleEng), "EN") FILTER (!LANGMATCHES(LANG(?title), "EN")) } c2(173.28. [START_REF]-deduced 7 times BGP[END_REF].114) : Query Q2 DBpedia's TPF server ?predicate = rdf s : label & ?object = "Brad P itt"@en . . . ?predicate = rdf s : label & ?object = "N atalie P ortman"@en . . . TPF clients decompose SPARQL queries into a sequence of triple pattern queries. Table 1 presents some traces of the TPF server for query Q 1 . Oddnumbered lines represent received triple pattern queries and even-numbered ones represent sent triples after evaluation on the RDF graph.

Lines 1 and 3, correspond to triple pattern queries for tp 2 and tp 1 of Q 1 .2 We can observe that the object in Line 3, comes from a mapping seen in Line 2. This injection of a mapping obtained from a previous triple pattern query, is clearly a join implemented in a nested-loop from tp 2 towards tp 1 .

As the TPF server only sees triple pattern queries, the original queries e.g., Q 1 and Q 2 are unknown to the data provider. In this work, we address the following research question: Can we extract the BGPs from a TPF server log?

In our example, we aim to extract two BGPs from the TPF server log, one corresponding to Q 1 , BGP [1]= {tp 1 .tp 2 .tp 3 .tp 4 .tp 5 } and another corresponding to Q 2 , BGP [START_REF] Gallego | An Empirical Study of Real-World SPARQL Queries[END_REF]= {tp 1 .tp 2 .tp 3 .tp 4 }. Before presenting our scientific problem, we introduce the following definitions.

Definition 1 (BGP). A BGP (Basic Graph Pattern) is a set of triple patterns. Any tuple

∈ (I ∪ L ∪ V ) × (I ∪ V ) × (I ∪ L ∪ V ) is a triple pattern s, p, o
, where I is the set of all IRIs, L the set of all literals and V the set of all variables disjoint from L and I. 3Definition 2 (TPF server log). A TPF server log is a totally ordered sequence of execution traces structured in tuples ip, ts, tp, µ o where ip is the IP address of the client, ts is the timestamp of the http request, tp is a triple pattern, and µ o is the set of RDF triples returned by the TPF server transformed in mappings.

We denote by E(Q i ), the log produced by a TPF server when evaluating the SPARQL query Q i and by E(Q 1 ... Q n ) the log of n concurrent queries.

Definition 3 (Approximation ≈ of BGPs). A BGP approximates another (≈) if, to some extent, both contain same triple patterns and same joins.

To measure such approximation we can use precision and recall of triple patterns and joins of one BGP against another. The average of precision and recall can be used as a measure of global quality of the approximation.

Definition 4 (Problem of BGP reversing). Given a log corresponding to the execution of one query, E(

Q i ), find a function f (E(Q i )) producing a set of BGPs {BGP 1 , ..., BGP n }, such that: Property 1. f (E(Q i ))

approximates (≈) the BGPs existing in the original query.

If we consider that BGP (Q i ) returns the set of BGPs of

Q i then f (E(Q i )) ≈ BGP (Q i ). Property 2. f (E(Q i )
) guarantees resistance to concurrency, i.e., BGPs obtained from the log of isolated queries, approximate (≈) results obtained from the log of concurrent queries:

f (E(Q 1 )) ∪ ... ∪ f (E(Q n )) ≈ f (E(Q 1 ... Q n )).
We evaluate the BGPs extracted by f with the precision, recall and quality of triple patterns and joins returned by f against those existing in original queries. If f (E(Q 1 )) extracts the BGP= {tp 1 .tp 2 .tp 3 .tp 4 .tp 5 }, then precision, recall and quality of triple patterns and joins are perfect according to the BGP present in Q 1 , even if variables have different names. But if f (E(Q 1 )) misses one triple pattern (e.g., tp 5 ), then precision is 4/4, recall is 4/5 and quality is (1 + 4/5)/2.

In Figure 1, if c 1 and c 2 have different IP addresses, it is possible to separate

E(Q 1 Q 2 ) into E(Q 1 ), E(Q 2 )
and apply the reversing function to each trace. However, in the worst case, c 1 and c 2 have the same IP address, i.e., a web application running on the cloud that executes queries Q 1 and Q 2 in parallel. Thus, we expect that

f (E(Q 1 Q 2 )) ≈ f (E(Q 1 )) ∪ f (E(Q 2 )).

LIFT: a reversing function

We propose LIFT as a system of heuristics to implement f . The idea is to detect nested-loop joins. In Table 1, the mappings returned in Line 2 are reused in the next triple pattern query at Line 3. We track such bindings in order to link variables of different triple pattern queries. In this paper, we make the following hypothesis: (i) we consider only bound predicates, (ii) we do not consider the server's web cache (this information can be easily obtained by data providers), and (iii) we do not consider the client's TPF cache. The first hypothesis can be omitted, but we kept it because analysis of SPARQL queries show that predicates are frequently bound [START_REF] Gallego | An Empirical Study of Real-World SPARQL Queries[END_REF]. Figure 2 presents a simplified log of E(Q 3 ), E(Q 4 ) and E(Q 3 Q 4 ) where: Q 3 = SELECT * W HERE {?x p2 toto . ?x p1 ?y} and Q 4 = SELECT * W HERE {?x p3 titi . ?x p1 ?y . ?x p4 tata}.

For the sake of simplicity, timestamps are transformed into integers. The IP address of the TPF client is the same for Q 3 and Q 4 , so we removed the ip column. Variables are named ?s or ?o. µ o represents the mappings of variables resulting from the evaluation of tp. We call them output-mappings. Observe the client first requests more selective triple patterns, i.e., ?x p2 toto for Q 3 and ?x p3 titi for Q 4 , leaving at the end less selective ones i.e., ?x p1 ?y for both queries. Then mappings returned by selective patterns are bound into less selective ones producing a nested-loop. See that mappings c1 and c2 are bound in the variable ?x of the second triple pattern of Q 3 . Similarly, mappings c3 and c4 are bound in the variable ?x of the other triple patterns of Q 4 . We call inputmappings these injected mappings. Modified triple patterns are the inner part of the nested-loop that we call inner loop. We call outer loop the triple patterns whose mappings are used to bound variables, e.g., ?x p2 toto .

The basic intuition of LIFT is to detect if mappings obtained in a request are bound in next requests. This can be challenging because mappings can be : (i) bound several times (e.g., in star queries), (ii) bound partially as a side-effect of LIMIT and FILTER clauses, (iii) or bound into a different concurrent query.

As a real log can be huge, LIFT analyzes the log on a sliding window defined by a gap, i.e., a time interval. When LIFT reads an entry e in the log with a timestamp ts, it considers only entries reachable within the gap i.e., ts ± gap. Algorithm 1 shows the three phases of LIFT. This algorithm can aggregate triple patterns that do not belong to the same nested-loop as it is the case in our example of Figure 3, where CTP [3] aggregates triple patterns of Q 3 and Q 4 . We suppose that this case is not likely, especially when the gap is small but if it is the case, next algorithm splits ctps to separate these nested-loops. If the µ i of a ctp is a subset of the µ o of a previous ctp, then we consider that the 2 corresponding variables can be linked. This happens in the example described in Figure 4a with CTP [START_REF] Gallego | An Empirical Study of Real-World SPARQL Queries[END_REF] and CTP [START_REF] Lorey | Detecting SPARQL Query Templates for Data Prefetching[END_REF]. We consider that ?σ of CTP [START_REF] Lorey | Detecting SPARQL Query Templates for Data Prefetching[END_REF] is linked to ?s of CTP [START_REF] Gallego | An Empirical Study of Real-World SPARQL Queries[END_REF]. We formalize this behavior at Lines 6 to 7 of Algorithm 3. A direct inclusion does not occur if Algorithm 2 aggregated too many log entries as it is the case with CTP [3]. Q 3 and Q 4 have a common triple pattern ?x p1 ?y and Algorithm 2 aggregates them. We solve this problem by splitting a ctp. The idea is to produce a dtp from a ctp c k if it exists an intersection between the µ o of a ctp c l and the µ i of c k . In the example described in Figure 4a, CTP [3] produces two splits: one when analyzing CTP [1] (DTP [3] is produced) and another when analyzing CTP [START_REF] Gallego | An Empirical Study of Real-World SPARQL Queries[END_REF] (DTP [START_REF] Lorey | Detecting SPARQL Query Templates for Data Prefetching[END_REF] is produced) because both µ o intersect the µ i of CTP [3]. Splitting affects input-mappings, but also timestamps and output-mappings. After splitting, we obtain input-mappings that are subsets of previous output-mappings. Intersection and splitting is shown in Lines 5 and 8 of Algorithm 3.

Algorithm 1: Global algorithm of LIFT

Nested-loop join detection

foreach (c k , vi) ∈ { (c k , vi) | c k ∈ CTP ∧ c k .id > c.id ∧ ingap(c k ,c,gap) ∧ ∃ vi ∈ vars(c k .µi) | c k .µi(vi) ∩ c.µo(vo) = ∅ } do 6 if c k .µi(vi) ⊆ c.µo(vo) then
For lack of space, we do not detail the function split. It basically groups, from the TPF log, values that belong to the intersection, this generates correct timestamps, output-mappings and input-mappings. We register the split relationship with a split predicate that links a ctp with its produced dtps. In our example, we have 2 split relations; split (3,3) and split (3,[START_REF] Lorey | Detecting SPARQL Query Templates for Data Prefetching[END_REF]. Splitting has an effect on CTP traversal that we see in Line 3. Output-mappings of produced dtps must be analyzed so when the nested-loop detection analyzes a splitted ctp, it inserts in the CTP List, produced dtps. split(c) returns the set of dtps produced by splitting c.

BGP Extraction

Figure 5 represents the connected components of the DTP Graph shown in Figure 4b. From this representation, it is easy to compute the final BGPs with a variable renaming and restitution of an IRI/literal in place of ω when there is only one input mapping, e.g., toto, titi and tata. Detected joins are underlined.

In our example, LIFT rebuilds perfectly BGPs of queries Q 3 and Q 4 . This example is executed with gap = 8. If we reduce the gap, then some joins are not detected and recall decreases. If we execute concurrently more queries having same triple patterns, then LIFT can deduce joins that do no exist in original queries and consequently precision will decrease. In Section 4, we measure experimentally the precision and recall of LIFT in different situations.

Experiments

The goals of the experiments are twofold, (i) to evaluate precision and recall of LIFT's results and (ii) to show that LIFT extracts meaningful BGPs from a real TPF log.

First, we analyze a small set of queries taken from the TPF web site. We evaluate precision and recall of LIFT deductions from logs of queries executed in isolation in Section 4.1, and from logs of queries executed concurrently in Section 4.2. These two sections evaluate to which extent LIFT guarantees Properties 1 and 2 of our problem statement (cf. Definition 4), correspondingly.

Then, we analyze LIFT with an important number of real user queries taken from logs of USEWOD 2016 [5]. In Section 4.3, we evaluate precision and recall of 14,259 queries coming from logs of the DBpedia's SPARQL endpoint. In Section 4.4, we analyze 4,720,874 single triple pattern queries coming from logs of the DBpedia's TPF server.

Evaluation of LIFT with queries in isolation

We used 29 over 30 queries of the TPF web site 4 . Concerned datasets are DBpedia 2015-04, UGhent, LOV or VIAF. We captured http requests and answers of queries using the webInspector 1.2 tool5 . Source code of LIFT is available on GitHub 6 For each query Q i , we run LIFT(E(Q i ), ∞). Figure 6 presents precision and recall of LIFT deductions against original queries 7 , they show to which extent LIFT(E(Q i )) ≈ BGP (Q i ) (cf. Definition 4, Property 1). In average, LIFT obtained 97% of recall and 75% of precision of joins. That gives a quality of 84.5%. LIFT deduces perfectly 15 of 30 BGPs:

Q 1 -Q 6 , Q 9 , Q 11 , Q 15 -Q 18 ,
Q 22 , and Q 29 -Q 30 . Q 28 was not analyzed because it needs two TPF servers. LIFT does not detect UNION queries because they are processed on the client side. Q 9 is a query with a BGP like {(tp1 UNION tp2) . tp3}. In this case, LIFT detects 2 BGPs, {tp1 . tp3} and {tp2 . tp3}. Q 29 is also a UNION query but without join, as there is no intersection between mappings, LIFT detects two separate BGPs containing one triple pattern. We consider this behaviour correct.

Figure 7 describes Q 7 and its deduced BGPs. BGP [1] is correct, while BGP [START_REF] Gallego | An Empirical Study of Real-World SPARQL Queries[END_REF] is not. This is due to traces of the first request made by TPF clients for each triple pattern of each query to decide the join ordering (TPF servers send cardinalities of concerned triples in their answers). Thus, when processing Q 7 , the client makes a first request of each triple pattern and then decides to begin with the first triple pattern. Then it binds resulting mappings into the ?book variable of the second triple pattern to retrieve corresponding authors. This nested-loop is deduced in BGP [1]. But as output mappings of the first request of the second triple pattern intersects with the values of the inner loop, LIFT deduces BGP [START_REF] Gallego | An Empirical Study of Real-World SPARQL Queries[END_REF] with a self-join that is very unlikely and that can be easily filtered in a post-processing. Such situation appears in 6 of 29 queries:

Q 7 , Q 12-14 , Q 21 and Q 25 .
In some cases, LIFT deduces additional triple patterns and thus false joins with well deduced triple patterns because of the intersection of mappings. This is more challenging to filter. In Figure 7, Q 8 has an additional triple pattern, the last one, and a join with the second triple pattern. This is the case for

Q 8 , Q 10 , Q 14 , Q 20 , Q 23-27 .
In addition, LIFT merges triple patterns when they are very similar as it is the case in Q 19 and Q 20 where some triple patterns have same predicate and variables in the same position (subject/object).

ID Original query

Deduced BGPs

Q7

SELECT DISTINCT

?book ?author WHERE { ?book rdf:type dbpo:Book . ?book dbpo:author ?author } LIMIT 100 BGP [1]: {?s1 rdf :type dbpo:Book . ?s1 dbpo:author ?o2} BGP [START_REF] Gallego | An Empirical Study of Real-World SPARQL Queries[END_REF]: {?s3 dbpo:author ?o3 . ?s3 dbpo:author ?o4}

Q8

{SELECT ?award WHERE { ?award a dbpedia-owl:Award . ?award dbpprop:country ?language . ?language dbpedia-owl:language dbpedia:Dutch_language} {?s1 dbpedia-owl:language dbpedia:Dutch_language . ?s2 dbpprop:country ?s1 . ?s2 rdf :type dbpedia-owl:Award . ?s1 rdf :type dbpedia-owl:Award} Fig. 7: LIFT deductions for Q 7 and Q 8 . Prefix dbpo corresponds to dbpedia-owl.

Does LIFT results resist to concurrency?

We implemented a tool to shuffle several TPF logs according to different parameters. Thus, given E(Q 1 ), ..., E(Q n ), we are able to produce different significant representations of E(Q 1 ... Q n ). We grouped queries targeting the same dataset into a set of randomly chosen queries as shown in Table 2. For each query set, we evaluate how

LIFT(E(Q 1 ), gap) ∪ ... ∪ LIFT(E(Q n ), gap) ≈ LIFT(E(Q 1 ... Q n ),
gap) in terms of precision and recall for different gap values (cf. Definition 4, Property 2). Gap varies from 1% to 100% of the log duration. Each query set was shuffled 4 times and we calculate the average of results by gap. In Figure 8 we report results by join, (a) and (b) show precision whereas (c) and (d) show recall.

Dataset

Query sets

DBpedia 2015

DB1 = {Q 1 , Q 8 , Q 14 , Q 22 } DB4 = {Q 4 , Q 12 , Q 24 } DB2 = {Q 3 , Q 11 , Q 15 , Q 20 } DB5 = {Q 7 , Q 16 , Q 21 , Q 5 } DB3 = {Q 6 , Q 13 , Q 19 , Q 27 } DB6 = {Q 9 , Q 10 , Q 29 , Q 30 } UGhent U G = {Q 2 , Q 23 , Q 25 , Q 29 , Q 30 } LOV LV = {Q 17 , Q 18 , Q 26 , Q 29 , Q 30 } VIAF V F = {Q 29 , Q 30 }
Table 2: Query sets executed concurrently over a TPF server.

Precision and recall globally improve when gap increases. When gap is small (less than 50%) precision decreases significantly. A small gap leads LIFT to split values of an inner loop i.e., the ctpExtraction algorithm cannot aggregate in one ctp all triple patterns of the inner loop.

Concerning recall, LIFT is moderately impacted by concurrency. In some cases, LIFT favours recall by producing all possible joins in the nested-loop detection.

Concerning precision, LIFT is more impacted by concurrency and results depend on concurrent executed queries. When executed queries have triple patterns that are semantically or syntactically similar, then LIFT generates many false joins that impact precision. A post-processing could filter these false joins.

Evaluation of LIFT with real SPARQL queries

The goal of this experiment is to evaluate to which extent LIFT is able to deduce BGPs of an important number of real user queries. We analyzed 10 hours of one day (2015-10-30) of the log of the DBpedia SPARQL endpoint of USEWOD 2016 dataset [5]. From 380,834 http requests containing SPARQL queries, we analyzed 14,259 queries that represent our ground truth. We kept only SELECT queries having BGPs with more than one triple pattern and that return not null results. We made an aggregation of same queries and kept only one, by hour and by user. That is because, if in a TPF log, the same query appears several times, LIFT deduces only one query, i.e., LIFT(E(Q 1 ) E(Q 1 ), gap) ≈ BGP (Q 1 ). Then, we constraint to 200 the number of queries by user for each hour. The OPTIONAL operator was transformed in a JOIN and we filter queries having a join over predicates. 8We installed locally a TPF server with DBpedia 3.9 to execute and collect the TPF server logs that were given as input to LIFT. Gap was set to one hour. We calculate precision, recall and quality of deduced BGPs against the ground truth. LIFT returned BGPs in xml files which were loaded in a MySQL database. This allows us to make several analysis.

Globally, we obtained 69% of precision, 64% of recall and 66% of quality. Results are better for recurrent queries (in particular for precision) for which we obtained 85% of precision, 55% of recall and 70% of quality. That means that queries executed frequently under different execution contexts (i.e., concurrent queries) are better deduced than queries requested few times. We analyzed results by user and calculate their averages, we obtained 84% of precision, 74% of recall and 78% of quality. From this analysis, we observed that there are particular users whose queries are less well deduced, but in general LIFT deductions by user are good.9

Analysis of a real TPF log

The goal of this experiment is to extract BGPs from the log of the DBpedia's LDF server of USEWOD 2016 dataset [5]. This log contains http requests from October 2014 to November 2015. We analyzed the first quarter of the log representing 4,720,874 single triple pattern queries (until 27th February 2015). We pruned 1% of the log corresponding to entries that are not well formed TPF requests. LIFT needs not only received triple pattern queries but also corresponding answers. To obtain answers, we executed the triple pattern queries of the log using a TPF client 10 . Then, we run LIFT with log slices of one hour with a maximum gap (one hour).

We obtained 1236 BGPs containing more than one triple pattern. Table 9 describes the most recurrent queries. Unsurprisingly, some of them are the queries available on the TPF web site. BGP [1] corresponds to Q 1 , BGP [START_REF] Gallego | An Empirical Study of Real-World SPARQL Queries[END_REF] is like BGP [1] except that dbpedia-owl:starring is replaced by dbpprop:starring. BGP [3] corresponds to the query used as the motivation example of [13], BGP [START_REF] Lorey | Detecting SPARQL Query Templates for Data Prefetching[END_REF] corresponds to Q 3 , BGP [5] to Q 6 , etc. In this top 20 list, only few queries were unknown: BGP [START_REF] Möller | Learning from Linked Open Data Usage: Patterns & Metrics[END_REF], BGP [7], BGP [8], BGP [12] and BGP [13]. We can analyze deduced BGPs to obtain some statistics, for instance, the most common type of join is subjectsubject (2693), followed by subject-object ( 1063), what is coherent with analysis shown in [START_REF] Gallego | An Empirical Study of Real-World SPARQL Queries[END_REF]. LIFT's deductions from the TPF log of USEWOD 2016 are available in an xml document 11 .

Compared to the state of art, [START_REF] Verborgh | Initial Usage Analysis of DBpedia's Triple Pattern Fragments[END_REF] reports statistics on TPF logs that allow to verify server availability, number of evaluated requests and cache usage. These statistics are useful but can be obtained at triple pattern granularity. LIFT, is able to extract BGPs from TPF logs, that can be subsequently analyzed to retrieve other information, e.g., frequently executed BGPs as done in Sections 4.3 and 4.4.

Extracting information from logs is traditionally a data mining process [3]. Sequential pattern mining [7] focuses on discovering frequent subsequences from an ordered sequence of events. However, searching for BGPs cannot be reduced to searching for frequent episodes in a TPF log. Suppose, two queries Q 1 : {?x p1 o1 . ?x p2 ?y} and Q 2 : {?x p1 ?y . ?y p3 ?z}. The TPF query engine executes the joins with a nested-loop. So, ?x p1 o1 and ?x p1 ?y will appear once in the log, while patterns with IRIs p2 ?y and IRIs p3 ?z will appear many times according to the selectivity of the triple patterns on p1. Searching for frequent episodes will raise up episodes with p2 and p3 but joins were between p1, p2 and p1, p3. Clearly, TPF logs need to be pre-processed before analysis.

In previous work [9], we focused on analyzing SPARQL logs. We proposed FETA, an algorithm to reverse BGPs of federated SPARQL queries from logs collected from a set of SPARQL endpoints. LIFT and FETA share the problem statement, but have several major differences. First, FETA analyzes a SPARQL log gathered from many data providers in order to infer BGPs of federated queries. LIFT analyzes a single TPF log of a single data provider to infer BGPs of SPARQL queries. This makes possible to run LIFT on real available logs as those published by USEWOD 2016 [5]. Second, the input log of FETA is composed of SPARQL queries that are more complex than TPF queries. On the other hand, SPARQL queries include variable names that can be used to disambiguate some situations. In TPF logs, such variables are not available. Next, TPFs contain metadata that helps clients to decide a join ordering. To obtain the cardinality of a triple pattern, a client actually requests the triple pattern making indistinguishable these requests from regular query processing. This makes BGP extraction from TPF logs more challenging.

Conclusions and future work

LIFT extracts BGPs from TPF server logs. It tracks joins executed with nestedloops following a strategy of merging and splitting compatible triple pattern queries. LIFT deductions have good precision and good recall mainly when evaluating frequent queries. The main challenge for LIFT is the concurrent execution of queries from the same client targeting exactly the same triples (queries sharing IRIs).

LIFT opens several perspectives. First, LIFT focused on nested-loops. Clients can also rely on symmetric hash to perform joins that can also be detected. Detecting nested-loop and symmetric hash together from the server-side is challenging. Second, in LIFT, we favoured recall, mainly by splitting input mappings for every intersection with output mappings. We can introduce more constraints for splitting that should improve precision against recall, e.g., to separate repeated executions of identical queries which traces appear in the same TPF. Third, LIFT is currently an algorithm that post-processes logs. We can transform LIFT into a streaming algorithm able to extract BGPs in real-time. This will also allow to process very large TPF logs. Finally, LIFT is able to process logs on different datasets by just merging them. In this case, one data provider is able to detect cross-join against different datasets hosted on the same server. If different data providers are ready to collaborate by just sharing their logs, they should be able to infer federated running across multiple servers.

Fig. 1 :

 1 Fig. 1: Concurrent execution of queries Q 1 and Q 2 .

Fig. 2 :

 2 Fig. 2: Examples of simplified TPF logs.

Fig. 3 :

 3 Fig. 3: TPF log and CTP List produced by Algorithm 2 with E(Q 3 Q 4 ) and gap = 8.

Algorithm 3

 3 describes how to link variables of different ctps produced by Algorithm 2. It builds a DTP Graph of deduced triples patterns (dtp) where nodes have the same structure as ctps and edges represent a relation of inclusion between input-mappings (µ i ) and output-mappings (µ o ) of 2 different dtps. Figure 4b presents the DTP Graph produced by Algorithm 3 with the CTP List of Figure 4a. Dashed links represent linked variables deduced by Algorithm 3. id ts tp µo µi 1 1,1 ?s p2 ?ω ?s c1 c2 ?ω toto 2 2,2 ?s p3 ?ω ?s c3 c4 ?ω titi 3 3,8 ?σ p1 ?o ?o a b c d ?σ c1 c2 c3 CTP List of candidate triple patterns obtained by Algorithm 2. id ts tp µo µi 1 1,1 ?s p2 ?ω ?s c1 c2 ?ω toto 2 2,2 ?s p3 ?ω ?s c3 c4 ?ω titi 3 3,5 ?σ p1 ?o ?o a b ?DTP Graph of deduced triple patterns obtained by Algorithm 3.

Fig. 4 :

 4 Fig. 4: CTP List and DTP Graph produced by Algorithm 3 with gap = 8.

Algorithm 3 : 2 foreach c ∈ CTP do 3 if 4 else

 3234 Detection of nested-loop joins 1 Function nestedLoopDetection (CT P, gap) is input : a list CT P of ctps; a gap in time units (seconds) output: An edge-labelled DTP Graph of dtps split(c) = ∅ then CTP.insertAfter(c.id,split(c)); DTP.addnode(c); foreach vo ∈ vars(c.µo) do 5

7 DTP 8 else

 78 .addnode(c k ) ; DTP.addEdge(c,c k ,(vo,vi)); DTP.addnode(s=split(c k , vi, c, vo)) ; DTP.addEdge(c,s,(vo,vi));

Fig. 5 :

 5 Fig. 5: Connected components of the DTP Graph produced by the execution of Algorithm 3 for gap = 8.

Fig. 6 :

 6 Fig. 6: Precision and recall by query of the TPF web site.

Fig. 8 :

 8 Fig. 8: Precision and recall.

Table 1 :

 1 Excerpt of log of the DBpedia's TPF server for query Q 1 .

  This allows to gather triple pattern queries that seem to be part of the same inner loop. 2. Next, LIFT looks for an inclusion relationship among output-mappings and input-mappings of ctps. If such an inclusion exists a deduced triple pattern (dtp) is created. If instead of inclusion an intersection exists, LIFT splits ctps to obtain a dtp with inclusion. If neither inclusion nor intersection exists an isolated dtp is created. This produces a DT P Graph where nodes are dtps and edges are inclusion relationships between dtps. 3. Finally, LIFT extracts BGPs from the DTP graph. Ideally, LIFT(E(Q 3 Q 4 ), gap) should compute the 2 BPGs of Q 3 and Q The objective of ctpExtraction is to aggregate together log entries that seem to participate in the same outer or inner loop. Aggregated entries are ctps. A ctp is a tuple ip, ts, tp, µ o , µ i where ip is an IP address, ts is a pair of timestamps (s.min, ts.max) representing a range; when creating a ctp both timestamps are identical and correspond to the timestamp of the corresponding entry in the log. tp is a triple pattern query, µAlgorithm 2 outlines the extraction of a CTP List from a TPF log with a particular gap. Figure3illustrates the effect of executing Algorithm 2 on log E(Q 3 Q 4 ) with gap=8. The log is processed in sequential order. Lines 5 to 7 initialize input-mappings by replacing constants by variables σ or ω. Next, Lines 9 to 10 merge (i.e., aggregate) current ctp with an existing and compatible one. An existing ctp is compatible if it has the same tp, it is produced by the same ip address, and fits in the gap. The ingap(c, c k , gap) function returns true if c.ts.min -c k .ts.max ≤ gap. If the current ctp is compatible with an existing one, output/input-mappings and timestamps are merged. When updating timestamps, the lower timestamp remains always the same, only the upper timestamp can grow. A variable of tp cannot belong to µ o and µ i simultaneously.

	ts tp	µo
	Section 3.1 details the CTP extraction. Section 3.2 describes the nested-loop detection. Finally, Section 3.3 presents the phase of extraction of BGPs. 1 ?s p2 toto ?s c1 c2 id ts tp µo µi 2 ?s p3 titi ?s c3 c4 1 1,1 ?s p2 ?ω ?s c1 c2 ?ω toto 3 c1 p1 ?o ?o a 2 2,2 ?s p3 ?ω ?s c3 c4 ?ω titi 4 c3 p4 tata 3 3,8 ?σ p1 ?o ?o a b c d ?σ c1 c2 c3 c4 5 c2 p1 ?o ?o b 6 c3 p1 ?o ?o c ?σ c3 c4 4 4,7 ?σ p4 ?ω ?ω tata 7 c4 p4 tata 3.1 Extraction of candidate triple patterns 8 c4 p1 ?o ?o d (b) CTP List of candidate triple patterns ob-
	(a) Log of a TPF server.	tained by Algorithm 2.

1 Function LIFT(log, gap) is input : a TPF server log; a gap in time units (seconds) output: a set of BGP s data : CT P a list of ctps, DT P a graph of dtps 2 CT P ← ctpExtraction (log, gap) 3 DT P ← nestedLoopDetection (CT P, gap) 4 return BGP ← bgpExtraction(DT P ) 1. First, LIFT merges triple pattern queries having same characteristics into candidate triple patterns (ctp). 4 : {?s p2 toto . ?s p1 ?o}, {?s p3 titi . ?s p1 ?o . ?s p4 tata}. o (output-mappings) is the list of solution mappings for variables of tp. µ i (input-mappings) is a set of mappings built during the ctpExtraction. Basically, we replace any constant of tp by a variable, we use σ for subject and ω for object. Replaced constants are regrouped in µ i . Algorithm 2: Extraction of Candidate Triple Patterns 1 Function ctpExtraction (log, gap) is input : a TPF server log; a gap in time units (seconds) output: a list CT P of ctp 2 CT P ←[ ] 3 foreach e ∈ log do 4 c ← read(e) as (ip,(ts,ts), tp, µo, µi) switch c.tp do 5 case ?s p o: do c.tp ← ?s p ?oin ; c.µi ← ?ω|o 6 case s p ?o: do c.tp ← ?sin p ?o ; c.µi ← ?σ|s 7 case s p o: do c.tp ← ?sin p ?oin ; c.µi ← ?σ|s, ?ω|o 8 if ∃ c k ∈ CT P | ingap(c,c k ,gap) ∧ c k .ip = c.ip ∧ c.tp = c k .tp then 9 c k . µo ∪ c. µo ; c k . µi ∪ c. µi ; c k .ts.max=c.ts.max; 10 else CTP.add(c) 11 return CT P

These queries come from http://client.linkeddatafragments.org/.

TPF clients only request bound parts of a triple patterns, variables are omitted.

We do not consider blank nodes in this paper.

http://client.linkeddatafragments.org/

https://sourceforge.net/p/webinspector/wiki/Home/

https://github.com/coumbaya/lift

Queries, TPF logs and LIFT results are available at: https://github.com/coumbaya/lift/blob/master/experiments.md

All scripts to process logs of the DBpedia SPARQL endpoint are in Python and available at https://github.com/edesmontils/BE4DBPedia/

These results are available at: http://documents.ls2n.fr/lift/dbpedia-sparql-results.tar

https://github.com/LinkedDataFragments/Client.js


Alexis Bitaillou

Benoît Parrein

Guillaume Andrieux

How to bring machine learning in industrial networks?

Keywords: Industrie 4.0, industrial networks, machine learning, cognitive networks, wireless network

Industry 4.0 opens new opportunities and challenges. One of these challenges is to provide a set of open and standardized technologies. The long life IT technologies such Ethernet and Wi-Fi can be considered as potential candidate. In this article, we present the industry 4.0 and related technologies like historical and candidate network technologies for the industry, IIoT and cognitive networks. We propose a concept using cognitive networks combined with machine learning to make IT networking technologies compliant with industrial requirements.

I. INTRODUCTION

To stay competitive, industries have to create value added. As more and more industries master industry 3.0, a new kind of industry is emerging: the industry 4.0. Industry 4.0 has several ambitions, for example to optimize resources. It based on Internet of Things (IoT) and Cyber-Physical Systems (CPS). Technically, industry 4.0 has opened new challenges. One of them is to adapt IoT to the industrial environment. Industry has particular constraint such as very low latency, high reliability and security. There are a lot of industrial network technologies. Even if operational technology (OT) and information technology (IT) are merging, the constraints and requirements are not the same. We will start by the technical background of our proposition.

II. BACKGROUND A. The fourth industrial revolution

The industry of future has been a subject of research since its first formalization in 2013 in Germany [START_REF] Kagermann | Recommendations for Implementing the Strategic Initiative INDUSTRIE 4.0: Securing the Future of German Manufacturing Industry ; Final Report of the Industrie 4.0 Working Group[END_REF]. Its name is different following the country. For examples, the industry of future is called Industrie 4.0 in Germany, Made in China 2025 in China and Nouvelle France Industrielle in France [START_REF] Kagermann | Recommendations for Implementing the Strategic Initiative INDUSTRIE 4.0: Securing the Future of German Manufacturing Industry ; Final Report of the Industrie 4.0 Working Group[END_REF], [START_REF] Li | China's manufacturing locus in 2025: With a comparison of "Made-in-China 2025[END_REF]. The main objective stills the same: make the industries more competitive by increasing the value-added and optimizing resources and processes. The common causes are the high cost of labour, more generally the high operating cost, and the competition of others countries' industries. The fourth industrial revolution is based on Industrial Internet of Things (IIoT), cyber-physical systems (CPS), manufacturing digitalization and integration of trends of IT (e.g. big data). We will take the case of German Industrie 4.0. Industrie 4.0 has several simple ambitions. The first step is to optimize efficiency of manufacturing processes and resources. For example, electricity could be saved by switching it off during breaks and week-ends. One of the ambitions is to use a unique set of standardized and open technologies. This recommendation is important because it lets companies choose the technologies they need. The others ambitions are about cyber-security, safety and in-house training for workers. Industry 4.0 opens new challenges such as smart manufacturing, the construction of network environment or CPS integration [START_REF] Zhou | Industry 4.0: Towards future industrial opportunities and challenges[END_REF], [START_REF] Lu | Industry 4.0: A survey on technologies, applications and open research issues[END_REF]. There are already researches on some of those topics. For example, the cloud manufacturing can be considered as a prototype of smart manufacturing [START_REF] Li | Cloud manufacturing: A new service-oriented networked manufacturing model[END_REF].

B. Industrial Internet of Things

According to the ITU-T Y2060/Y.4000, the Internet of Things (IoT) is "a global infrastructure for the information society, enabling advanced services by interconnecting (physical and virtual) things based on existing and evolving interoperable information and communication technologies" [START_REF] Itu | Overview of the Internet of Things[END_REF]. The Industrial Internet of Things (IIoT) is industrial version of IoT. The IIoT is composed of sensors, actuators and other cyber-physical systems. The IIoT has several purposes the industry's processes and management. Shrouf et al. [START_REF] Shrouf | Smart factories in Industry 4.0: A review of the concept and of energy management approached in production based on the Internet of Things paradigm[END_REF] proposed using IIoT to increase energy efficient of manufacturing. This approach consists in monitoring and analysing energy consumption. After the analysis, a strategy data-driven can be applied. Another element of Industry 4.0 is to limit down-time. Xu et al. [START_REF] Xu | Intelligent fault prediction system based on internet of things[END_REF] proposed using things like sensors to monitor production equipment and machines. Then they used tools like big data in order to predict fault. A preventive maintenance can be done before the predicted fault and potentially reduce the down-time. IIoT can help in a key point of Industry 4.0: the automation. The cyber-physical production system could be seen as the coupling of the robotics et IIoT [START_REF] Jazdi | Cyber physical systems in the context of Industry 4.0[END_REF]. The "IT world" has the challenges to connect the Industrial Internet of Things with the rest of the IT structure.

C. Networks

Historically, industrial networks are composed of wired elements. Technologies such as RS232 or CAN have been well tested for years. However, cables involve a low mobility and high deployment cost. These requirements are not fully compatibles with the industry of future and its needs of high flexibility [START_REF] Varghese | Wireless requirements and challenges in Industry 4.0[END_REF]. They still present in industries because they are considered as reliable, effective and safe. Alternatively, several wireless technologies have appeared. They are more and more reliable and faster and faster, IEEE 802.11 for example. There are a lot of potentially industrial technologies [START_REF] Wollschlaeger | The Future of Industrial Communication: Automation Networks in the Era of the Internet of Things and Industry 4[END_REF].

We start with wired technologies. In wired technologies, there are two main categories: the field buses and the Ethernetbased technologies. The field buses gather at least a dozen of standards [START_REF] Felser | The fieldbus war: History or short break between battles?[END_REF]. FIP, PROFIBUS, CAN and many others belong to this category. Since 2000, several field bus providers have incorporated Ethernet to their solution. Concretely, their solution consists of using Ethernet with their add-on. As Ethernet had not supported real time, the add-on provides real time and sometimes isochronous support. For example, SERCOS III, PROFINET and EtherCAT are Ethernet-based technologies. Since 2011, Ethernet has been improved to support real time but isochronous stills not supported.

We will start by low range wireless technologies. The main technologies in this category are RFID, IEEE 802.15.4 and Bluetooth. IEEE 802.15.4 is the base of WIA-PA (Wireless Networks for Industrial Automation-Process Automation), WirelessHART and ISA100.11a, three industrial wireless technologies [START_REF] Zurawski | Industrial communication technology handbook[END_REF]. RFID gathers several ISO standards, from electronic article surveillance to contactless smart card. For example, RFID can be used logistics for object tracking. With a range between 100 m and 250 m range, IEEE 802.11 is the only middle range wireless technology. However, IEEE 802.11 is not design to industrial networks. Its non-determinist nature and its "high" power consumption are its principal default for industrial usage [START_REF] Cheng | Det-WiFi: A Multihop TDMA MAC Implementation for Industrial Deterministic Applications Based on Commodity 802[END_REF]. IEEE 802.11ah is designed to machineto-machine (M2M) communications [START_REF] Adame | IEEE 802.11ah: The WiFi approach for M2M communications[END_REF]. It overcomes power consumption problem but the bandwidth is low (8 Mbit/s max.).

Historically, cellular networks are not design for industrial purpose. The power consumption is quite high, so devices on battery have short life. LTE-M, Narrowband-IoT (NB-IoT) and EC-GSM-IoT have been included in the LTE standard [START_REF] Liberg | Cellular Internet of things: technologies, standards, and performance[END_REF]. All of them are compatible with M2M. Their bandwidth is lower than "classical" LTE (e.g. 3 Mbit/s for LTE-M2), but they extend battery life. LTE-M, NB-IoT and EC-GSM-IoT can fully reuse already deployed 4G networks. This kind of technology is use in smart grid for example.

With the Internet of Things, a new kind of long range has appeared: the Low Power Wide Area Networks (LP-WAN) [START_REF] Mekki | A comparative study of LPWAN technologies for largescale IoT deployment[END_REF]. Their range comes from few kilometres up to 40 km. Their bandwidth rarely exceeds 1 Mbit/s. For example, LoRa uses sub-GHz ISM band (e.g. 868 MHz in EU) and 2.4 GHz ISM band. In LoRa, the maximum bandwidth is 50 kbit/s in sub-GHz ISM band and 2 Mbit/s at 2.4 GHz. The range is up to 5 km in urban area and up to 20 km in rural area. The other preponderant technology is Sigfox. The maximum range is 10 km in urban area and 40 km in rural area. Nevertheless, the maximal bandwidth is limited to 600 b/s in upstream and 100 b/s in downstream (respectively 140 and 4 messages per day). IEEE 802.15.4g and 802.15.4k are the only IEEE standards in this category. Other proprietary solutions exist such as Weightless, DASH7 and Ingenu RPMA.LPWAN can be used to check cows' health in agriculture for example.

D. Interferences

Several wireless technologies use the same ISM bands (sub-GHz and 2.4 GHz mainly). In some situations, it creates interferences. There are two kinds of interferences: i) multipath interferences and ii) multi-user interferences. Multipath interferences appear when a wave takes different paths from the sender to the receiver [START_REF] Saleh | A Statistical Model for Indoor Multipath Propagation[END_REF]. To solve this problem, Bottomley et al. [START_REF] Bottomley | A generalized RAKE receiver for interference suppression[END_REF] proposed to use RAKE antenna. Tsai and Chang [START_REF] Tsai | The feasibility of combating multipath interference by chirp spread spectrum techniques over Rayleigh and Rician fading channels[END_REF] proposed to use Chirp Spread Spectrum (CSS). For example, LoRa uses mainly CSS. These solutions get round the problem without cancelling multipath interferences.

The second kind of interferences is multi-user interferences (MUI). They happen when two sources broadcast in the same time. According Nobilet [START_REF] Nobilet | Etude et optimisation des techniques MC-CDMA pour les futures générations de systèmes de communications hertziennes[END_REF], orthogonal codes can avoid multi-user interferences. However, orthogonal codes can be broken by reflection. So, in case of reflection, orthogonal codes don't protect from MUI and multipath interferences appear. Another solution is simply to use time division multiplexing such as TDMA (Time Division Multiplexing Access). In theory, devices can broadcast at precise moment and not at the same time. In practices, TDMA is not efficient and not scalable.

Interferences are important because they are the causes of collision. Veijlgaard et al. [START_REF] Vejlgaard | Interference Impact on Coverage and Capacity for Low Power Wide Area IoT Networks[END_REF] shown in their experiment that interferences can increase lost packet rate up to 50% for LoRa and 60% for Sigfox. The coverage of LoRa and Sigfox also decreases inside and outside. Petrova et al. [START_REF] Petrova | Interference Measurements on Performance Degradation between Colocated IEEE 802.11g/n and IEEE 802.15.4 Networks[END_REF] shown that up to 90% of IEEE 802.15.4 can be discarded when IEEE 802.15.4 network and IEEE 802.11n network are concurrent. They also showed that changing the channel of IEEE 802.15.4 or IEEE 802.11n can solve the problem. This task can be achieved with cognitive radio.

E. Cognitive radio and cognitive network

The cognitive radio (CR) appeared in 1998 created mainly by Mitola [START_REF] Mitola | Cognitive radio[END_REF]. It builds on the top of Software-Defined Radio (SDR). The cognitive radio has been created to add smartness in SDR. The cognitive radio transforms radio nodes from simple protocol users into smart radio domain-aware agents. In cognitive radio, nodes can reconfigure their radio settings such as power delivered, frequencies used or modulation automatically depending on the situation. In order to make this objective, Mitola created RKRL (Radio Knowledge Representation Language), a description language. With RKRL, users can define a set of settings and conditions. The nodes follow the configuration according to the tested conditions. This previous definition has been completed by Haykin in 2005 [START_REF] Haykin | Cognitive radio: Brain-empowered wireless communications[END_REF]. According to Haykin, CR has to make communications reliable when needed, and make efficient radio spectrum usage. Haykin also proposed using cognitive radio to solve spectrum scarcity problem. This opens new challenges about sensing, detection of users, MAC strategies and routing. Cognitive networks (CN) extent the concept of radio cognitive to support wired networks and heterogeneous networks [START_REF] Thomas | Cognitive networks: Adaptation and learning to achieve end-to-end performance objectives[END_REF]. CR focuses on physical and link layer, whereas cognitive network has the OSI stack as scope. CN can perform multi-objectives optimization, unlike cross-layer design [START_REF] Dietterich | Machine Learning for Cognitive Networks: Technology Assessment and Research Challenges[END_REF]. In 2008, Manoj et al. [START_REF] Manoj | CogNet: A cognitive complete knowledge network system[END_REF] proposed CogNet, a cross-layer design cognitive network architecture.

III. PROPOSITION

Our objective is to help to build the network of industry of the future. We see that the ambitions of the industry 4.0 define some constraints such as using open and interoperable technologies or the need of flexibility. As industry 4.0 is based on (I)IoT and CPS, so we need consider their particularities. The choice of the network technologies is finally oriented by all these constraints. Our solution is to use cognitive network in order to provide low end-to-end delay, low latency and a high reliability. In a cognitive network, the parameters of the OSI stack can be finely tuned to provide better performance without change the infrastructure. Optionally, our solution can consider out-ofstack parameters such as energy and security requirements. For example, it is possible to adjust the power transmission and other parameters according a policy that maximizes performance and power saving. In order to maximize the efficiency and the effectiveness of the solution, we may adopt a crosslayer design. The project MobileMAN shown the cross-layer approach can increase performance [START_REF] Borgia | Mobileman: Design, integration, and experimentation of cross-layer mobile multihop ad hoc networks[END_REF]. Nevertheless, the cross-layer design has some drawbacks, especially on the longevity. The number of parameters and the number of value taken by the parameters could be very high. For example, there are about 6720 possible combinations just for LoRa [START_REF] Bor | LoRa Transmission Parameter Selection[END_REF]. The choice of optimal values of the cognitive network is so not trivial. Machine learning (e.g. Deep Neural Networks, SVR) can help to guess the best parameters according the objectives. The figure 1 illustrates our proposition.

IV. CONCLUSION

In this paper, we presented the industry 4.0 and related technologies such as IoT, industrial networks and cognitive network. We also proposed an architecture compliant with the objectives of the industry 4.0. Our concept can be mainly composed of IT technologies such Wi-Fi, Ethernet or machine learning. Our proposition is also independent of specific or proprietary technology. We plan to validate our architecture by simulations and experimentations. The simulations will be done with software like ns2. For example, FIT IoT-LAB could be a good testbed for this kind of experimentations.

Figure 1 .

 1 Figure 1. Proposed architecture inspired from CogNet [28].


Loris Di Cairano
email: l.di.cairano@fz-juelich.de

Matteo Gori
email: gori6matteo@gmail.com

Marco Pettini
email: pettini@cpt.univ-mrs.fr

Coherent Riemannian-geometric description of Hamiltonian order and chaos with Jacobi metric

Keywords: numbers: 05.20.Gg, 02.40.Vh, 05.20.-y, 05.70.-a Hamiltonian Chaos, Differential Geometry

By identifying Hamiltonian flows with geodesic flows of suitably chosen Riemannian manifolds, it is possible to explain the origin of chaos in classical Newtonian dynamics and to quantify its strength. There are several possibilities to geometrize Newtonian dynamics under the action of conservative potentials and the hitherto investigated ones provide consistent results. However, it has been recently argued that endowing configuration space with the Jacobi metric is inappropriate to consistently describe the stability/instability properties of Newtonian dynamics because of the non-affine parametrization of the arc length with physical time. To the contrary, in the present paper it is shown that there is no such inconsistency and that the observed instabilities in the case of integrable systems using the Jacobi metric are artefacts.

I. INTRODUCTION

Elementary tools of Riemannian differential geometry can be successfully used to explain the origin of chaos in Hamiltonian flows or, equivalently, in Newtonian dynamics. Natural motions of Hamiltonian systems can be viewed as geodesics of the configuration-space manifold M equipped with the Riemannian metric g J , known as the Jacobi metric (or kinetic energy metric). The stability/instability properties of such geodesics can be investigated by means of the Jacobi-Levi-Civita (JLC) equation for geodesic spread. It has been shown that chaos in physical geodesic flows does not stem from hyperbolicity of M : phase space trajectories/geodesics are destabilized both by regions of negative curvature and by parametric instability caused by positive curvature varying along the geodesics [START_REF] Pettini | Geometrical hints for a nonperturbative approach to Hamiltonian dynamics[END_REF][START_REF] Casetti | Geometric approach to Hamiltonian dynamics and statistical mechanics[END_REF][START_REF] Pettini | Geometry and Topology in Hamiltonian Dynamics and Statistical Mechanics[END_REF]. Another remarkable fact is that the JLC equation written for a geometrization of Hamiltonian systems in an enlarged configuration space-time endowed with a metric due to Eisenhart [START_REF] Eisenhart | Dynamical Trajectories and Geodesics[END_REF] yields the standard tangent dynamics equation commonly used in numerical computations of the largest Lyapunov exponent (LLE). These two different geometric framework have been proven to give the same information about order and chaos and about the strength of chaos as well. This has been checked in the case of two-degrees of freedom Hamiltonian systems, for the Hénon-Heiles model and for two coupled quartic oscillators, respectively [START_REF] Cerruti-Sola | Geometric description of chaos in two-degrees-of-freedom Hamiltonian systems[END_REF][START_REF] Pettini | On the Riemannian description of chaotic instability in Hamiltonian dynamics[END_REF], and in the case of a large number of degrees of freedom (from 150 up to 1000) [START_REF] Cerruti-Sola | Lyapunov exponents from geodesic spread in configuration space[END_REF]. It has been found that the JLC equation stemming from the Jacobi metric gives exactly the same quantitative results of the tangent dynamics equation. This notwithstanding, in Ref. [START_REF] Cuervo-Reyes | Non-affine geometrization can lead to non-physical instabilities[END_REF] it has been argued that the non-affine parametrization of the arc-length with time in configuration space endowed with the Jacobi metric leads to nonphysical instabilities. More precisely, the JLC equation written for the Jacobi metric seems to give chaos also for a system of harmonic oscillators. Moreover, such alleged nonphysical instabilities are found to be stronger in systems with few degrees of freedom. Such an argument seems to radically exclude the use of Jacobi metric in configuration space to consistently investigate Hamiltonian chaos; in fact, the non-affine parametrization of the arc length with respect to physical time is an unavoidable consequence of the way Jacobi metric is derived from Maupertuis' least action principle. Some mathematical works [START_REF] Montgomery | Who's afraid of the Hill boundary?[END_REF][START_REF] Giambo | Morse theory for geodesics in singular conformal metrics[END_REF][START_REF] Giambo | On the normal exponential map in singular conformal metrics[END_REF], partly motivated by the results reported in Ref. [START_REF] Cuervo-Reyes | Non-affine geometrization can lead to non-physical instabilities[END_REF], have investigated the behaviour of the geodesics in configuration space endowed with the Jacobi metric near the so called Hill's boundaries, i.e. the regions in configuration space where E = V (q) and where the Jacobi metric is singular (g J = 0). In particular, all these works emphasized the phenomenon of geodesic reflection near Hill's boundaries and the relevance of these reflections to characterize periodic orbits, with, among the others, an original proposal dating back to Ref. [START_REF] Seifert | Periodische bewegungen mechanischer systeme[END_REF]. Finally, other authors [START_REF] Yamaguchi | Geometric approach to Lyapunov analysis in Hamiltonian dynamics[END_REF] have tackled the geometrization of Hamiltonian dynamical systems by lifting the Jacobi metric from configuration space M to its cotangent bundle T * M , i.e.

to the phase space. In this framework JLC equations are rewritten as a system of first order linear differential equations on the tangent bundle T T * M of phase space: this allows to identify the adequate degrees of freedom to compute the Geometric Largest Lyapunov Exponent(GLLE). Within this framework, the GLLE (JLC equations) and LLE (tangent dynamics) are found in very good agreement to characterize the chaotic regime (strong/weak chaos) of the Hénon-Heiles model. All these studies, together with the manifest contradiction among the outcomes of Ref. [START_REF] Cuervo-Reyes | Non-affine geometrization can lead to non-physical instabilities[END_REF] and those of Refs. [START_REF] Cerruti-Sola | Geometric description of chaos in two-degrees-of-freedom Hamiltonian systems[END_REF][START_REF] Pettini | On the Riemannian description of chaotic instability in Hamiltonian dynamics[END_REF] have motivated the present work.

The present paper is organized as follows. In Section II we briefly discuss some aspects of the construction of the Jacobi metric, with emphasis on the consequences of non-affine parametrization of the arc length with time for the JLC equation, and on the presence of boundaries where the metric is singular. In Section III, the JLC equation is rewritten by introducing a parallel transported frame for which explicit expressions are then given in Subsections III.1 and III.3 for N = 2 and N = 3, respectively. Then in Subsections III.2 and III.4 the results of the corresponding numerical simulations are reported for two and three (resonant) harmonic oscillators, respectively, using the most critical values of the parameters which, according to Ref. [START_REF] Cuervo-Reyes | Non-affine geometrization can lead to non-physical instabilities[END_REF], should yield non physical instabilities. It is shown that this is not the case. In Section IV, we discuss how the measure concentration phenomenon, which takes place at a large number of degrees of freedom, completely removes any problem making unnecessary even resorting to a parallel transported frame. Some conclusions are drawn in Section V.

II. EFFECTS OF NON-AFFINE PARAMETRIZATION OF THE ARC LENGTH WITH JACOBI METRIC

Among the different possibilities of rephrasing Newtonian dynamics in geometric terms, as reported in [START_REF] Casetti | Geometric approach to Hamiltonian dynamics and statistical mechanics[END_REF][START_REF] Pettini | Geometry and Topology in Hamiltonian Dynamics and Statistical Mechanics[END_REF], the Jacobi metric in configuration space leads to the mathematically richest structure [(M, g J ) is geodesically complete in the sense of the Hopf-Rinow theorem].

Let us consider a thorough investigation of the geometrization of Newtonian dynamics by means of (M, g J ) for systems described by a Lagrangian of the form

L( q, q) = 1 2 g ij (q) qi qj -V (q) , (1)
where

qi (t) = dq i dt (t) (2)
g ij are the components of the kinetic energy metric on configuration space M with the associated Levi-Civita connection ∇ specified by the Christoffel coefficients Γ i jk , and V (q) is the potential energy. It is well known that, given a chart (U, φ) on the configuration space M and a curve γ : I ⊂ R -→ M , the natural motions φ(γ(t)) = γ(t) = q(t) are the class of curves that make stationary the action functional

S[q(t)] = t 1 t 0 L( q(t), q(t))dt (3)
on the class of curves with q(t 0 ) = a and q(t 1 ) = b fixed, i.e.

δS[γ] = 0 with γ(t 0 ) = a and γ(t 1 ) = b . (4)
The Newton equations are derived from the Euler-Lagrange equations, i.e.

∇ γ γ = -gradV (5)
where grad i f (q) = g ik ∂ k f is the gradient. Natural motions γ(t) belong to a class of curves of configuration space satisfying the "physical" variational principle (4), therefore these curves can be identified with geodesics of configuration space which also satisfy a variational principle but in this case of a "geometrical" kind. In fact, geodesics γ(s) are curves of a Riemannian manifold (M, g) endowed with a metric g that makes stationary the length functional between two fixed points, i.e.

δl [γ] = 0 with l [γ(s)] = s 1 s 0 gij dq i ds dq j ds ds γ(s 0 ) = a γ(s 1 ) = b . (6)
One possible way to provide an identification of natural motion with geodesics is provided by the introduction of the Jacobi metric on a subspace of configuration space. Let us consider the total energy function

H( q, q) = qi ∂L ∂ qi -L = 1 2 g ij (q) qi qj + V (q) (7)
which is obviously a conserved quantity along the natural motions, i.e. dH( q(t), q(t))/dt = 0.

Since the Lagrangian is a homogeneous function in qi it follows that H( q, q) = 2W -L where

W = 1 2 g ij qi qj = H( q, q) -V (q) (8)
is the kinetic energy. If we consider a class of isoenergetic trajectories q(t; E) in configuration space, that is having the same total energy value H( q(t; E), q(t; E)) = E, as the kinetic energy is non negative, the trajectories of the system in configuration space are confined in the region M V <E = {q ∈ M |V (q) < E}. Moreover, for isoenergetic trajectories the kinetic energy can be expressed as a function of the coordinates, i.e.

W (q(t; E)) = E -V (q) (9)
thus the action functional of Eq. ( 3) can be rewritten in the form

S E [q(t; E)] = t 1 t 0 L( q(t; E), q(t; E)) dt = t 1 t 0 [E + 2W (q(t; E))] dt = (t 1 -t 0 )E + t 1 t 0 2W (q(t; E)) dt . (10)
as we are interested in the variational principle and t 0 , t 1 , E are fixed quantities, the first term in the last equality can be neglected. The integral in Eq.( 10) can be interpreted as a length integral in configuration space, in fact

S E [q(t; E)] = t 1 t 0 2W (q(t; E)) dt = t 1 t 0 2W (q(t; E)) 2W (q(t; E)) dt = t 1 t 0 2 [E -V (q)] ḡij (q) qi qj dt = = t 1 t 0 2 [E -V (q)] ḡij qi qj dt = t 1 t 0 g ij qi qj dt = t 1 t 0 ds dt 2 dt = s 1 s 0 ds ( 11
)
where the new metric

g ij := 2W (q)g ij , W (q) := E -V (q) ( 12
)
called also Jacobi metric, has been introduced with the associated the arc-length element

ds 2 = g ij dq i dq j = 2 [E -V (q)] g ij dq i dq j = 4[W (q)] 2 dt 2 . ( 13
)
A central point of the following discussion is related to the non-affine parametrization of the arc-length s with respect to the physical time t, which is clearly a necessary consequence of the construction of Jacobi metric. Moreover, we observe that the Jacobi metric is related to kinetic energy metric g ij through a conformal rescaling via a factor proportional to the kinetic energy [E -V (q)] preserving the signature of the metric only in the interior of

M V ≤E , i.e. M V <E = M E = {q ∈ M |V (q)
< E}, the so called Hill's region. By endowing the region M E with the Jacobi metric g the natural motions with fixed energy E are the same as geodesics γ(s) of the manifold (M E , g).

This approach has remarkable consequences, as it is discussed in the following. In fact, the geometric description of Newtonian dynamics, identifying the solutions of Newton equations with the geodesics of suitable Riemannian manifolds, provides a powerful conceptual and mathematical framework to study the stability/instability of dynamics in terms of the stability properties of a geodesic flow, described by the geodesic spread equation relating stability/instability with geometry.

The standard observable to define the presence of dynamical chaos and to measure its strength is the largest Lyapunov exponent, and, as we will see in the following, in the geometrical framework a Geometrical Lyapunov exponent can be defined. Now, a key point in derivation of Jacobi metric from Maupertuis' principle is to set the relation between the physical time t and the arc-length s as

ds 2 = 4W 2 (q)dt 2 . ( 14
)
It follows that a generic geodesic parametrized by the arc-length s, i.e. q(s) = { qi (s)} i∈ [1,n] has a unit velocity

g J ( q(s), q(s)) = 2W dq i ds δ ij dq j ds = 1
while if parametrized with respect to the physical time

dq i dt g ij dq j dt = 2W . (15)
So, in general, the physical time is a non-affine parametrization of the geodesics of Jacobi metric. We derive in what follows the effect of the reparametrization (14) on the geodesic equation. Let us introduce the vector fields

Y i := dq i ds (16)
and

X i := dq i dt = ds dt dq i ds = 2W Y i (17)
defined along the geodesic q i (s) = q i (s(t)) = qi (t). Using the definition of geodesic for the

vector field Y ∇ Y Y = 0 (18)
it is possible to derive the equations for the vector field X = αY (with α = (2W ) -1 )

∇ αX (αX) = α (∇ X α) + α 2 X = 0 (19) that implies ∇ X X = -∇ X (log α)X = X∇ X (log 2W ) . (20)
In a natural coordinate system {q i } i=1,...,N the equation (20) reads

d 2 q i dt 2 + Γ i jk dq j dt dq k dt = δ i k dq k dt dq j dt ∂ j log(2W ) . (21)
The Christoffel symbols in Jacobi metric take the form

Γ i jk =Γ i jk + 1 2 δ i j ∂ k log(2W ) + δ i k ∂ j log(2W ) -g jk grad i log(2W ) (22)
where Γ i jk are Christoffel symbols of the kinetic energy metric g ij and grad i f := g ik ∂ i f is the gradient with respect to kinetic energy metric. Substituting ( 22) in (21) and using (15)

we obtain

∇ X X = d 2 q i dt 2 + Γ i jk dq j dt dq k dt = 1 2 grad i (2W ) = -grad i V (23)
i.e. the Newton's equations of the dynamical system. As already mentioned above, dynamical chaos can now be investigated by means of the equation for the geodesic spread describing the stability of a geodesic flow of a Riemannian manifold. The geodesic spread is measured by a vector field J which locally gives the distance between nearby geodesics. This vector field evolves along a reference geodesic according to the Jacobi-Levi Civita equation which, in components, reads

∇ 2 J i ds 2 + R i jkl dq j ds J k dq l ds = 0 ( 24
)
where

R i jkl = ∂ k Γ i jl -∂ l Γ i jk + Γ m lj Γ i km -Γ m kj Γ i ml ( 25
)
is the Riemann curvature tensor. In Jacobi metric we have

∇ 2 X i ds 2 = 1 2W ∇ dt 1 2W ∇X i dt = 1 4W 2 ∇ 2 X i dt 2 + 1 4W ∇X i dt d dt 1 W (26)
that substituted into (24) yields to

∇ 2 J i ds 2 + R i jkl dq j ds J k dq l ds = 0 → ∇ 2 J i dt 2 + R i jkl dq j dt J k dq l dt = ∇J i dt d dt ln W . (27)
Every non-affine parametrization in Jacobi metric can be reduced to a relation between the physical time and the arc-length parameter of the form t → ds = f (q)dt, where f (q) is a function of the point and generates a term proportional to the derivative of f (q) into each equations. Thus, the JLC equation written for the Jacobi metric is not invariant for time reparametrization. We can rephrase the main point raised by the work in Ref. [START_REF] Cuervo-Reyes | Non-affine geometrization can lead to non-physical instabilities[END_REF] as attributing to the right hand-side of Eq.( 27) the origin of chaotic-like instabilities even for integrable systems, that is, the origin of non-physical artifacts. In Ref. [START_REF] Cuervo-Reyes | Non-affine geometrization can lead to non-physical instabilities[END_REF] it has been surmised that the larger the fluctuations of kinetic energy W the more dramatic the occurrence of non-physical instabilities stemming from Eq.( 27). Remarkably, the geometrization of Newtonian dynamics in an enlarged configuration space-time equipped with the Eisenhart metric tensor [START_REF] Pettini | Geometrical hints for a nonperturbative approach to Hamiltonian dynamics[END_REF] yields the JLC equation in the form of the Tangent Dynamics Equation which is commonly used to compute the Largest Lyapunov Exponent (LLE) [START_REF] Pettini | Geometrical hints for a nonperturbative approach to Hamiltonian dynamics[END_REF][START_REF] Pettini | Geometry and Topology in Hamiltonian Dynamics and Statistical Mechanics[END_REF], therefore the authors of Ref. [START_REF] Cuervo-Reyes | Non-affine geometrization can lead to non-physical instabilities[END_REF] claim that this is the only consistent geometrization of Newtonian dynamics to investigate chaos, while the Jacobi metric would be unsuitable for the same task.

II.1. Geometrical Lyapunov exponent

We conclude the present Section by giving a definition of a Geometrical Lyapunov exponent. Of course, the starting point is the JLC which, in intrinsic notation, is

∇ 2 ξ J + R(J, ξ)ξ = 0 . ( 28
)
By defining Y = ∇ ξ J and R ξ J := R(ξ, J)ξ the above equation becomes

∇ ξ J = Y ∇ ξ Y = R ξ J . (29)
Then, by putting

A :=   0 1 1 R ξ 0   , J :=   J Y  
the JLC equation reads

∇ ξ J = AJ . (30)
A Geometrical Lyapunov exponent, in analogy with the definition of the standard Lyapunov exponent, can be defined after having expressed as a function of physical time the solution of Eq. ( 30) as

λ G := lim t-→∞ 1 t J (t) g J (0) g . ( 31
)
where the norm of J is

J (t) 2 g = 2 W (q)δ ij {J i (t)J j (t) + (∇ qJ ) i (t)(∇ qJ ) j (t)} .
Let us note that throughout the literature, [START_REF] Pettini | Geometrical hints for a nonperturbative approach to Hamiltonian dynamics[END_REF], [START_REF] Pettini | Geometry and Topology in Hamiltonian Dynamics and Statistical Mechanics[END_REF], [START_REF] Cerruti-Sola | Geometric description of chaos in two-degrees-of-freedom Hamiltonian systems[END_REF] and [START_REF] Cerruti-Sola | Lyapunov exponents from geodesic spread in configuration space[END_REF], the geometrical Lyapunov exponent has been expressed as a function of physical time t, whereas in [START_REF] Cuervo-Reyes | Non-affine geometrization can lead to non-physical instabilities[END_REF] the Geometrical Lyapunov exponent is defined as a function of the arc-length. As a final comment, let us remark that the existence of many different frameworks to rephrase Newtonian dynamics in geometric terms [START_REF] Pettini | Geometry and Topology in Hamiltonian Dynamics and Statistical Mechanics[END_REF] can lead, a-priori, to different quantitative evaluations of the strength of chaos, to the contrary, the use of different geometric frameworks must lead to the same qualitative description of the stability/instability properties of the dynamics. For example, the transition from weak to strong chaos must be at least qualitatively and coherently reproduced in any geometric framework, as is for example actually shown for high dimensional Hamiltonian flows in Ref. [START_REF] Cerruti-Sola | Lyapunov exponents from geodesic spread in configuration space[END_REF].

III. PARALLEL TRANSPORTED FRAME FOR A SYSTEM OF N PARTICLES IN JACOBI MANIFOLD

Let us now work out the JLC equation for a parallel transported orthonormal frame along a reference geodesic. The advantage of this representation with respect to the use of standard local coordinates is that by making parallel transported frames to "incorporate" the geodesics reflection, when they approach the Hill's boundaries in configuration space, eliminates a source of artefacts in the numerical solution of the JLC equation. In fact, the sharp reflection of geodesics close to the Hill's boundary of a mechanical manifold would require a prohibitively high numerical precision to avoid the introduction of an error amplification mimicking chaos even for integrable systems. To the contrary, with respect to a parallel transported frame -"incorporating" the geodesic reflection -nearby geodesics are no longer affected by the fake error amplification due to the reflection. As a matter of fact, we will show that the solutions of the JLC equation -written for the Jacobi metric -have the correct physical meaning also in the "pathological" cases where unphysical instabilities were found in Ref. [START_REF] Cuervo-Reyes | Non-affine geometrization can lead to non-physical instabilities[END_REF].

The parallel transported frame is built by requiring that the covariant derivative ∇ ξ X of all vectors X with respect to the geodesic flow ξ is orthogonal to ξ. Let us introduce a reference frame e i := ∂/∂q i on the tangent bundle T M and the corresponding dual frame

θ i := dq i such that θ i (e j ) = δ i j .
The Jacobi metric tensor is the multi-linear map g J : T M × T M -→ R such that, if written with respect to the natural basis, is given by

g J := g ij θ i ⊗ θ j (32)
the components of which are

g ij := g(e i , e j ) = 2(E -V (q))δ ij
with the differential arc-length

ds 2 = g ij θ i θ j . (33)
To build a parallel transported frame we consider a geodesic γ : I ⊂ R -→ M , with tangent vector field ξ : I -→ T M and the Jacobi vector field J : I -→ T M such that respect to the natural basis {e i } i∈ [1,n] they are written as

J = J i e i ξ = ξ i e i
where the Jacobi field, by definition, verifies

∇ 2 ξ J + R(J, ξ)ξ = 0 (34)
We shall show that there exist reference frames parallel transported along any geodesic. Such systems exist if an orthogonal tensor field Ω : I -→ SO N (γ(I)) is defined at each point along the geodesic flow. This tensor field allows to define the required frame where the basis on T M is given by {E i } i=i,••• ,N and dual frame on T * M {Θ i } i=i,••• ,N and, by definition of dual frame we have Θ i (E j ) = δ i j . Therefore, the parallel transported frame is defined by

E i = Ω • e i , ∇ ξ E i = 0, g(E i , E j ) = δ ij Θ i = θ i • Ω -1 , ∇ ξ Θ i = 0, g * (Θ i , Θ j ) = δ ij (35)
In the Jacobi-Levi Civita equation, written in (34), we can define the Ricci tensor along the geodesic flow Ric(ξ) : γ(I) -→ T 1 1 (γ(I)) by

Ric(ξ) := R(•, ξ)ξ . (36)
A canonical isomorphism exists between the rank-two tensor Ric(ξ) and a symmetric matrix N × N , that we shall denote with R. For every ω ∈ T * M and X ∈ T M , such a matrix is given by

Ric(ξ)(X, ω) := Ric(ξ) i j X j ω i = g(Y ω , R(X)) (37)
where Y ω := g -1 ω. We note the following symmetry properties of Ric(ξ)

Ric(ξ) = R jkl i ξ k ξ l e i ⊗ θ j = R jkli ξ k ξ l θ i ⊗ θ j = R ilkj ξ k ξ l θ i ⊗ θ j = R ilk j ξ k ξ j e j ⊗ θ i (38)
and by redefining the indices j i and l k the symmetry of Ric(ξ) is evident. Hence the associated matrix can be diagonalized at each point along the flow. Let note that the symmetry of Riemann tensor entails

Ric(ξ)ξ = R(ξ, ξ)ξ = 0 . (39)
Then, the tangent vector to the geodesic ξ is an eigenvector of the matrix R associated to Ric(ξ) with vanishing eigenvalue. A posteriori, one observes that for harmonic oscillatorsgeometrized through the Jacobi metric -the eigenbasis of the matrix R coincides with the parallel transported frame. In general, this is not true and we have two different orthogonal matrices, that one which diagonalises R and another one which transforms the natural basis into the parallel transported frame. Fortunately, in the present case, given the natural frame

{e i } i∈[1,N ]
on the tangent bundle T M , there exists an orthogonal transformation Ω : I -→ SO N (γ(I)) such that

E i = Ω • e i (40)
namely, such that it transforms the natural basis into the parallel transported frame

{E i } i∈[1,N ]
and, moreover, such that it diagonalises the matrix R:

R d = Ω -1 • R • Ω (41)
This allows to write the Jacobi field with respect to such a basis and thus

J = J i e i = J i (Ω -1 ) k i Ω l k e l = [(Ω -1 ) k i J i ][Ω l k e l ] = Jk E k . ( 42
)
Let us consider the JLC equation ( 34) and proceed to substitute J → J then

∇ 2 ξ J + Ric(ξ) J = 0 d 2 Ji ds 2 E i + Jk Ric(ξ) l k E l = 0 d 2 Jk ds 2 + Jk Λ k E k = 0 ( 43
)
where now the repeated indices do not stand for summation. In this way, we obtain N -1 second order differential equations in the unknown functions

{ Jk } i∈[1,N -1]
, and the N -th equation is that for the vector tangent to the reference geodesic, equation which corresponds to the geodesic equation thus, by denoting with J1 the function for this equation, we have

d 2 J1 ds 2 = 0 d 2 Jk ds 2 + Jk Λ k = 0 (44)
It is interesting to remark that {Λ k } are just sectional curvatures, i.e the principal directions of curvature identified by the vectors

{E k } k∈[1,N ] , i.e Λ k := E k , R(E k , E 1 )E 1
Now, passing from the arc-length parameter s to the physical time t the equations (44

) become d 2 Jk dt 2 (t) - d dt log(W (q(t))) d Jk dt (t) + 4W (q(t)) 2 Jk (t)Λ k (t) = 0 (45) whence Jk (t) -a(t) Jk (t) + b(t)J k (t) = 0 k ∈ [2, N ] J1 (t) -a(t) J1 (t) = 0 a(t) := d dt log(W (q(t))) b k (t) := 4W (q(t)) 2 Λ k . (46)
The second equation in 46 can be immediately solved by setting h(t) := J1 (t) and then tackling the differential equation

ḣ(t) + a(t)h(t) = 0 =⇒ h(t) = W (q(t))
giving

J 1 (t) = t t 0 W (q(η)) dη . (47)
The other equations can be written in canonical form, namely, by redefining the function J as follows

J = A(t)f (t) (48)
and substituting it into (46) we obtain

f (t) + ḟ (t) -a(t) + 2 Ȧ(t) A(t) + f (t) -a(t) Ȧ(t) A(t) + Ä(t) A(t) + b(t) = 0 . ( 49
)
By choosing A(t) such that the coefficient of ḟ (t) vanishes, we get the following conditions for A(t)

Ȧ A = a(t) 2 Ä A = a(t) 2 4 + ȧ(t) 2
(50) thus giving

A(t) = exp 1 2 a(t) dt = exp 1 2 d dt log(W (q(t))) dt = exp 1 2 log(W (q(t))) = W (q(t)) . (51)
The equation (49) becomes

f (t) + f (t) - a(t) 2 4 + ä(t) 2 + b(t) = 0 f (t) + f (t)   - 3 4 Ẇ W 2 + 1 2 Ẅ W + b(t)   = 0 (52)
To apply this procedure to the equations (44), we set

Jk := A(t)f k (t) (53)
and for every k ∈ [2, N ] we obtain the final form for the components of the Jacobi-Levi

Civita equation f k + ω (k) (t)f k (t) = 0 (54)
where the functions in the second term of the l.h.s. can be interpreted as time dependent frequencies

ω (k) (t) = 4W (q(t)) 2 Λ k (t) - 3 4 Ẇ W 2 + 1 2 Ẅ W (55)
where

Ẇ = -V = -ξ, ∇ R n V R n Ẅ = -V = -HessV (ξ, ξ) + ∇ R n V 2 R n (56)
III.1. Parallel Transported Frame for a system of 2 harmonic oscillators

In Refs. [START_REF] Cerruti-Sola | Geometric description of chaos in two-degrees-of-freedom Hamiltonian systems[END_REF][START_REF] Pettini | On the Riemannian description of chaotic instability in Hamiltonian dynamics[END_REF] it has been shown that for the Hnon-Heiles model and for two coupled quartic oscillators, respectively, the geometrization through the Jacobi metric perfectly discriminates between ordered and chaotic motions by investigating the stability/instability of geodesics through the JLC equation expressed in a parallel transported frame. In this section, we are going to show that for two harmonic oscillators (of course an integrable system) the norm of the geodesic separation vector remains bounded, in spite of the fluctuations of kinetic energy which, according to the claim of Ref. [START_REF] Cuervo-Reyes | Non-affine geometrization can lead to non-physical instabilities[END_REF], should have entailed apparent instability of the regular motions of this system. The Hamiltonian of this system is

H(p 1 , p 2 , q 1 , q 2 ) = 1 2 (p 2 1 + p 2 2 ) + κ 2 (q 2 1 + q 2 2 ) (57)
and the associated Jacobi metric, having set W (q) := E -V (q 1 , q 2 ), is

g J = 2W (q)(dq 1 ⊗ dq 1 + dq 2 ⊗ dq 2 ) . ( 58
)
The Ricci tensor along the geodesic flow is

[Ric(ξ)] i j = Eκ W 2   ξ 2 ξ 2 -ξ 1 ξ 2 -ξ 1 ξ 2 ξ 1 ξ 1 .   (59)
The eigenvectors of this matrix are

E 1 = (ξ 1 , ξ 2 ) t E 2 = (-ξ 2 , ξ 1 ) t (60)
with the corresponding eigenvalues:

Λ 1 = 0 Λ 2 = Eκ 2W 3 . (61)
With these eigenvectors the matrix for the basis transformation is simply obtained in the form

Ω =   ξ 1 -ξ 2 ξ 2 ξ 1   (62)
so that the JLC equation for the components of the parallel transported Jacobi vector field J is written as

d 2 J1 ds 2 = 0 d 2 J2 ds 2 + Eκ 2W 3 J2 = 0 (63)
which, written for the physical time t and with the notations of the preceding Section, read

d 2 J 1 dt 2 - Ẇ W dJ 1 dt = 0 d 2 f 2 dt 2 + ω2 f 2 = 0 (64)
where

ω2 (t) = 2Eκ W - 3 4 Ẇ W 2 + 1 2 Ẅ W (65)

III.2. Numerical results

In Ref. [START_REF] Cuervo-Reyes | Non-affine geometrization can lead to non-physical instabilities[END_REF], it has been claimed that there exists a set of initial conditions for which the JLC equations written for the Jacobi metric lead to unstable solutions even in the case of two harmonic oscillators because of the affine parametrization of the arc-length with time and the consequent fluctuations of the kinetic energy. By starting from the general solutions for a system of two harmonic oscillators given by

q 1 (t) = A 1 cos(ωt + φ 1 ) q 2 (t) = A 2 cos(ωt + φ 2 ) ( 66
)
where A 1 , A 2 and φ 1 , φ 2 are determined by the initial conditions, the authors of Ref. [START_REF] Cuervo-Reyes | Non-affine geometrization can lead to non-physical instabilities[END_REF] used polar coordinates to rewrite the previous equations in the compact form

r(t) = R 2 + ∆ 2 cos(2ω t)
where

R 2 = A 2 1 + A 2 2 2 , ∆ 2 = A 2 1 -A 2 2 2 , ( 67
)
then they have reported that every initial condition fulfilling the condition

I = ∆ R 4 > 4 7
(68) yields unstable solutions. We have adopted the same initial conditions for two harmonic oscillators and representing the solutions of the JLC equations with respect to a parallel transported frame, already used in [START_REF] Cerruti-Sola | Geometric description of chaos in two-degrees-of-freedom Hamiltonian systems[END_REF], and we have found that the equations (64) with the frequencies (65) display stable solutions. We have solved the equations by using two different conditions both fulfilling Eq. ( 68). These are

I 0 = ∆ R 4 = 56 81 (69)
I 1 = ∆ R 4 = 21 25 ( 70
)
Let us consider the condition (69). This is obtained by plugging through Eqs.(66) the initial condition

A 1 → 1 2 , A 2 → 1 10 9 + 2 √ 14 , φ 1 → 0, φ 2 → π 2 (71)
into the frequency (65). The time variation of ω2 (t) is reported in Figure 1. Correspondingly, the time dependence of the Geometric Lyapunov Exponent, λ(t), reported in Figure 2, clearly displays the typical time dependence found for regular trajectories, that is, λ(t) ∼ 1/t.

Coming now to the second condition given in (70), again this is implemented by plugging through Eqs.(66) the initial condition

A 1 → 1 2 , A 2 → 1 4 -5 + √ 21 , φ 1 → 0, φ 2 → π 3 ( 72
)
into the frequency (65).

The time variation of ω2 (t) is now reported in Figure 3. The corresponding time variation of the Geometric Lyapunov Exponent, λ(t), is now reported in Figure 4. Again it is found that λ(t) decays as 1/t, as expected for regular motions.

III.3. Parallel Transported Frame for a system of 3 harmonic oscillators

A relevant step forward is now obtained by considering three degrees of freedom because now the geodesic separation vector has two nontrivial components (since the component parallel to the velocity vector does not accelerate). Therefore, we consider 3 harmonic oscillators described by the Hamiltonian

H(p 1 , p 2 , p 3 , q 1 , q 2 , q 3 ) = 1 2 (p 2 1 + p 2 2 + p 2 3 ) + κ 2 (q 2 1 + q 2 2 + q 2 3 ) (73)
and the corresponding Jacobi metric, again with W (q) := E -V (q 1 , q 2 , q 3 ), is

g J = 2W (q)(dq 1 ⊗ dq 1 + dq 2 ⊗ dq 2 + dq 3 ⊗ dq 3 ) (74)
The eigenvectors of the Ricci tensor Ric(ξ, ξ) are

Ẽ1 = ξ =      ξ 1 ξ 2 ξ 3      Ẽ2 =      ξ 1 (q 2 ξ 2 + q 3 ξ 3 ) -q 1 (ξ 2 ξ 2 + ξ 3 ξ 3 ) ξ 2 (q 1 ξ 1 + q 3 ξ 3 ) -q 2 (ξ 1 ξ 1 + ξ 3 ξ 3 ) ξ 3 (q 1 ξ 1 + q 2 ξ 2 ) -q 3 (ξ 1 ξ 1 + ξ 2 ξ 2 )      Ẽ3 =      q 3 ξ 2 -q 2 ξ 3 q 1 ξ 3 -q 3 ξ 1 q 2 ξ 1 -q 1 ξ 2      (75)
with the associated eigenvalues

Λ 1 = 0 Λ 2 = Eκ 2W 3 Λ 3 = - κ 4W 2 ξ 1 ξ 1 (-4E + 3κq 2 2 + 3κq 2 3 ) + ξ 2 ξ 2 (-4E + 3κq 2 1 + 3κq 2 3 ) -6κq 2 q 3 ξ 2 ξ 3 + ξ 3 ξ 3 (-4E + 3κq 2 1 + 3κq 2 2 ) -6κq 1 ξ 1 (q 2 ξ 2 + q 3 ξ 3 ) , (76) respectively.
With these eigenvectors the matrix for the basis transformation now is

Ω =      ξ 1 ξ 1 (q 2 ξ 2 + q 3 ξ 3 ) -q 1 (ξ 2 ξ 2 + ξ 3 ξ 3 ) q 3 ξ 2 -q 2 ξ 3 ξ 2 ξ 2 (q 1 ξ 1 + q 3 ξ 3 ) -q 2 (ξ 1 ξ 1 + ξ 3 ξ 3 ) q 1 ξ 3 -q 3 ξ 1 ξ 3 ξ 3 (q 1 ξ 1 + q 2 ξ 2 ) -q 3 (ξ 1 ξ 1 + ξ 2 ξ 2 ) q 2 ξ 1 -q 1 ξ 2      . (77)
In order to make the parallel transported reference frame orthonormal, we have to orthogonalise the above eigenvectors with respect to the Jacobi metric, thus for i = 1, 2, 3 we have

n 2 i g J ( Ẽi , Ẽi ) = 1 (78)
namely

n i = 1 2W δ αβ E α i E β i (79)
These factors are

n 2 1 = 1 2W (q) n 2 2 = 1 (q 2 2 + q 2 3 )(ξ 1 ) 2 -2q 1 q 2 ξ 1 ξ 2 + (q 2 1 + q 2 3 )(ξ 2 ) 2 -2q 1 q 3 ξ 1 ξ 3 -2q 2 q 3 ξ 2 ξ 3 + (q 2 1 + q 2 2 )(ξ 3 ) 2 n 2 3 = (ξ 1 ) 2 + (ξ 2 ) 2 + (ξ 3 ) 2 ((ξ 1 ) 2 + (ξ 2 ) 2 )q 2 3 -2q 1 q 3 ξ 1 ξ 3 -2q 2 ξ 2 (q 1 ξ 1 q 3 ξ 3 ) + q 2 2 ((ξ 1 ) 2 + (ξ 3 ) 2 ) + q 2 1 ((ξ 2 ) 2 + (ξ 3 ) 2 ) (80)
The reference frame is composed by the normalized vectors

E i = n i Ẽi (81)
and by representing the Jacobi vector field as J = J i E i , we write the Jacobi-Levi Civita equations for the three components as

d 2 J 1 ds 2 = 0 d 2 J 2 ds 2 + Λ 2 J 2 = 0 d 2 J 3 ds 2 + Λ 3 J 3 = 0 (82)
Finally, passing to the physical time and by using Eq. ( 53) we get

f k + ω (k) (t)f k (t) = 0 with, again, ω (k) (t) = 4W (q(t)) 2 Λ k (t) - 3 4 Ẇ W 2 + 1 2 Ẅ W .
The results reported in the following are worked out by numerically integrating the following equations

d 2 J 1 dt 2 - Ẇ W dJ 1 dt = 0 d 2 f 2 dt 2 + ω (2) f 2 = 0 d 2 f 3 dt 2 + ω (3) f 3 = 0 (83)
with the following expressions for ω (2) (t) and ω (3) (t) :

ω (2) (t) := 2Eκ W - 3 4 Ẇ W 2 + 1 2 Ẅ W ω (3) (t) := 4W 2 Λ 3 - 3 4 Ẇ
W 2 + 1 2 Ẅ W (84)

III.4. Numerical results

In Ref. [START_REF] Cuervo-Reyes | Non-affine geometrization can lead to non-physical instabilities[END_REF], the alleged definitive argument to rule out the use of Jacobi metric to consistently describe the stability/instability of Hamiltonian dynamics was given by considering N decoupled harmonic oscillators. The claim was that non-vanishing fluctuations of kinetic energy (due to the non affine parametrization of the arc length with time) entail parametric resonance in the JLC equation mimicking chaos for an integrable system. The authors considered the solutions of this system in the form

q k (t) = cos (ωt + θ k ) , k = [1, N ] , (85)
where

θ k = k 2πf N
, with the phases θ k distributed on a fraction f of the interval 2π. It has been reported that the smaller f the larger fluctuation of kinetic energy and the larger the Lyapunov exponent. The fluctuation of kinetic energy √ σ is given in [START_REF] Cuervo-Reyes | Non-affine geometrization can lead to non-physical instabilities[END_REF] by

√ σ = sin(2πf ) √ 2N sin(2πf /N ) . (86)
Notice that in the N → ∞, the kinetic energy fluctuation magnitude √ σ has a non-vanishing value, so that the authors claim that any dimension this basic integrable system would display non-physical instabilities. In the next Section we will argue against this claim on the basis of an argument related with the concentration of measure at high dimension.

As shown in the preceding Section, and before in Refs. [START_REF] Cerruti-Sola | Geometric description of chaos in two-degrees-of-freedom Hamiltonian systems[END_REF][START_REF] Pettini | On the Riemannian description of chaotic instability in Hamiltonian dynamics[END_REF], a consistent description of order and chaos is obtained using the Jacobi metric and writing the JLC equation for a parallel transported frame which is quite simple to be found for N = 2, while the first non trivial extension is given by the N = 3 case. The JLC equation for the three dimensional case is given by Eqs.(83). There are three principal directions of curvature, namely, the sectional curvatures [Eqs. ( 76)], identified by the planes generated by the velocity vector along a geodesic, E 1 = ξ, and the parallel transported basis vectors E k with k = 2, 3.

These sectional curvatures coincide with the eigenvalues of the operator Ric(ξ, ξ) = R(•, ξ)ξ;

one of these is obviously zero because Ric(ξ, ξ)ξ = 0 while the other two are given by

Λ k = g(E k , Ric(ξ, ξ)E k ).
In Figure 6 Equations ( 83) have been integrated with a fourth-order Runge-Kutta algorithm along the q k (t) given by (85) and setting the phases θ k uniformly distributed on a fraction f of the interval 2π. It is well evident that the norm of the Jacobi geodesic separation vector is always bounded, coherently with the decay with 1/t of the running value of λ(t). The JLC equation written for the Jacobi metric and with a parallel transported frame provides the correct result: no instability of the trajectories of an integrable system is found, contrary to the claim of Ref. [START_REF] Cuervo-Reyes | Non-affine geometrization can lead to non-physical instabilities[END_REF].

IV. CONCENTRATION OF MEASURE OF THE VOLUME OCCUPIED BY AC-CESSIBLE CONFIGURATIONS WITH JACOBI METRIC

We know from the work in Ref. [START_REF] Cerruti-Sola | Lyapunov exponents from geodesic spread in configuration space[END_REF] that the Jacobi-Levi Civita equation, written in the natural reference frame for a large number of degrees of freedom, appropriately works by Hamiltonian is

� �� ��� ���� �� � �� � �� -� �� -� ���� � � λ � =���� (�)
H(P , Q) = P 2 R N 2 + Q 2 R N 2 = N i=1 1 2 |p i | 2 + |q i | 2 (87)
and the Jacobi metric in the Hill's region M E is

g J = 2 E - Q 2 R N 2 δ ij dq i ⊗ dq j . (88)
The associated Riemannian volume form is

ν J = det g J dq 1 ∧ • • • ∧ dq N (89)
where the determinant of the metric is

det g J = 2 N E - Q 2 R N 2 N . (90)
Therefore, the total volume of M E is given by the following integral

V N (E) = 2 N/2 M E E - Q 2 R N 2 N/2 dq 1 • • • dq N (91)
It is worth noticing a remarkable property of the Jacobi metric, its associated volume is proportional (up to an N -dependent factor) to the microcanonical ensemble measure

Ω N,µ (E) = M E |dp 1 ∧ . . . ∧ dp N ∧ dq 1 ∧ . . . ∧ dq N | (92)
where

M E = {(P , Q) ∈ T * M E |H(P , Q) = E}.
Due to the quadratic form of the kinetic energy 2W = R 2 , the integral (92) can be rewritten as

Ω N,µ (E) = M E A(S N -1 ) R≤ √ 2(E-V (Q)) R (N -1) dR dq 1 . . . dq N = = A(S N -1 ) N M E [2 (E -V (Q))] N/2 dq 1 . . . dq N = A(S N -1 ) N V N (E) , (93)
where A(S N -1 ) is the area of the unitary sphere. Let us now consider the volume of M E given by integral V N (E) and show that it concentrates around an N -1-dimensional manifold

Σ V = Q ∈ M E |V (Q) = V ,
in other words, the overwhelming contribution to the volume integral is given by microscopic configurations far from the boundary of M E . Although it could be questionable to provide a statistical argument for an integrable system for which ergodic hypothesis does not hold, the statistical averaging is intended as an averaging over all the possible configurations compatible with the constraint H(P , Q) = E. It is convenient to rewrite the integral (91) with the volume element expressed in spherical coordinates

|dq 1 ∧ • • • ∧ dq N | = Q N -1 sin N -2 (φ 1 ) sin N -3 (φ 2 ) • • • sin(φ N -2 )dφ 1 • • • dφ N -2 dQ (94)
where Q = Q R N = 2V (Q), because integrating over the angular variables one obtains

V N (E) = 2 N/2 A(S N -1 ) √ 2E 0 exp N 2 ln E - Q 2 2 + (N -1) ln Q dQ = 2 (N -1)/2 A(S N -1 ) E 0 dV exp N 1 2 ln (E -V ) + 1 2 - 1 N ln V = 2 (N -1)/2 A(S N -1 )E N 1 0 dx exp [-N F (x)] (95)
where x = V /E is the relative value of the potential energy with respect to the total energy and

F (x) = - 1 2 ln (1 -x) - 1 2 - 1 N ln x . (96)
As we are interested in the limit of large N , we can apply the Laplace approximation to evaluate the previous integral, i.e. we consider the Taylor expansion around the minimum with respect to x of F in the interval (0, 1),

1 0 exp [-N F (x)] dx ≈ exp [-N F (x)] 1 0 exp - (x -x) 2 2σ 2 dx . (97)
where σ 2 = (F (x)N ) -1 > 0. For a generic value of N , the solution of

F (x) = 0 is x = 1 2 N -2 N -1 , (98)
which is actually a minimum since From (99) and σ 2 = (F (x)N ) -1 , it follows that the largest part of the volume (∼ 99.7%) is concentrated around Σ xE in the interval [x -3σ, x + 3σ] with

F (x) = 4(N -1) 3 N 2 (N -2) > 0 for N > 2 . ( 99
σ = 1 2 √ N     1 - 2 N 1 - 1 N     1/2 ≈ 1 2 √ N . ( 100
)
This statistical argument shows that in the case of a large number of degrees of freedom the volume of the manifold is concentrated around a submanifold constant energy hypersurface V = E/2, far from the boundary of the Hill region where Jacobi metric is singular.

For example, with N = 100 kinetic energy fluctuations of absolute value of 10% occur with a probability of 68% while the probability of configurations hitting the boundary (E -V = 0) is ∼ 5.52 × 10 -88 .

V. DISCUSSION

Even though the point raised in Ref. [START_REF] Cuervo-Reyes | Non-affine geometrization can lead to non-physical instabilities[END_REF] is interesting, the conclusion put forward by the authors is incorrect. The fluctuations of kinetic energy along a trajectory/geodesic of the Jacobi metric associated with an integrable system, like a collection of harmonic oscillators, are by no means responsible for the activation of parametric instability mimicking a chaotic behaviour. When the number of degrees of freedom of a Hamiltonian system is small, the associated geodesics can often approach the boundary ∂M E = {q ∈ M |V (q) = E} of the mechanical manifold M E , and, in so doing, the geodesics bounce on ∂M E . The sharp reflection of the geodesics on the so-called Hill's boundaries [START_REF] Montgomery | Who's afraid of the Hill boundary?[END_REF][START_REF] Giambo | Morse theory for geodesics in singular conformal metrics[END_REF][START_REF] Giambo | On the normal exponential map in singular conformal metrics[END_REF][START_REF] Seifert | Periodische bewegungen mechanischer systeme[END_REF] are at the origin of numerical instabilities which in principle could be perhaps avoided by a prohibitively high precision of the integration algorithm for the Jacobi-Levi-Civita equation describing the geodesic spread. However, throughout this paper we have shown that this problem can be fixed by choosing a parallel transported coordinate system. The stability/instability of geodesics is an intrinsic property thus in principle independent of the choice of the coordinate system, however, not all the coordinate systems are necessarily equivalent from the point of view of their numerical implementation and reliability of the corresponding outcomes. And, in fact, the sharp reflection of the geodesics by the boundaries ∂M E is accounted for by a sudden reflection of the coordinate axes of the parallel transported frames thus separating the true geometric origin of stability/instability of geodesics from the source of numerical artefacts related with their peculiar shape.

We have then shown that when the number of degrees of freedom increases, then the probability of approaching the boundary of the corresponding mechanical manifold M E gets lower and lower and, even if at finite N the kinetic energy fluctuates it does not affect the strength of chaos measured through the outcomes of the JLC equation written for both the Jacobi and Einsenhart metrics which are in perfect agreement, as shown in Ref. [START_REF] Cerruti-Sola | Lyapunov exponents from geodesic spread in configuration space[END_REF]. Thus already for a few tens of degrees of freedom the JLC equation for (M E , g J ) written in natural chart [START_REF] Pettini | Geometry and Topology in Hamiltonian Dynamics and Statistical Mechanics[END_REF][START_REF] Cerruti-Sola | Lyapunov exponents from geodesic spread in configuration space[END_REF]]

d 2 J k dt 2 + 1 E -V ∂ k V δ ij dq i dt -∂ j V dq k dt dJ j dt + [∂ 2 kj V ] J j + 1 E -V (∂ k V )(∂ j V ) -∂ 2 ij V + (∂ i V )(∂ j V ) E -V dq i dt dq k dt J j = 0 .
can be safely used, at most with the exclusion of a zero measure set of initial conditions.

For very weakly coupled harmonic oscillators and N = 128 these equations give λ as small as 2 × 10 -6 .

In conclusion, the study of order and chaos of Hamiltonian flows -identified as geodesic flows of the Jacobi metric in configuration space -is legitimate and coherent, although not

FIG. 1 .FIG. 2 .FIG. 3 .FIG. 4 .

 1234 FIG. 1. Time dependence of the function ω2 (t) defined in (65), for the condition (69).

  Figures5 and 6, in Figures7 and 8, and in Figures9 and 10, respectively.

FIG. 5 .FIG. 6 .FIG. 7 .FIG. 8 .FIG. 9 .FIG. 10 .

 5678910 FIG. 5. Comparison between the GLE (blue curve) and 1/t (red straight line). f = 0.05, σ(0.05) = 0.6968.

)

  This means that the largest part of the volume V N (E) is concentrated around the hypersurface at constant potential energy Σ xE , a result very close to what is expected from the virial theorem [14].

Boltzmann prescription for the microcanonical partition function

This fact stands on the edge of a long standing debate about the "correct" prescription for the microcanonical partition function. Such a debate is out of the scope of the present work: we just report that Boltzmann prescription gives a result more consistent with the virial theorem in the considered case.