Dalal Madakat
email: dalal.madakat@onera.fr

Jérôme Morio
email: jerome.morio@onera.fr

Daniel Vanderpooten

A biobjective Branch and Bound procedure for planning spatial missions

Keywords: Space debris, Biobjective time dependent traveling salesman problem, Pareto Frontier, Fathoming rules

More than 90% of the space objects orbiting around the earth are space debris. Since the orbits of these debris often overlap the trajectories of spacecraft, they create a potential collision risk. The problem of removing the most dangerous space debris can be modeled as a biobjective time dependent traveling salesman problem (BiTDTSP). In this paper, we study an approach based on a branch and bound procedure to determine the Pareto frontier of the BiTDTSP.

Introduction

Since the launch of the first satellite in 1957, an increasing number of objects has been put into orbit around the earth. The number of space objects with diameter above 10 cm is estimated, nowadays, at about 15,000. Only 10% of these objects are operational, the other objects constitute space debris. Even if future launches are suspended, the number of space debris will continue to increase, e.g. owing to collisions between debris, making the collision with an operational satellite more probable [START_REF] Liou | Collision activities in the future orbital debris environment[END_REF]. Since consequences of collisions with debris may prove dramatic, avoidance maneuvers or missions to remove debris are necessary [START_REF] Crowther | Orbital debris: a growing threat to space operations[END_REF]. As removing all space debris would be quite expensive, the idea is to determine the debris which are more likely to cause collisions and remove them at least cost. Removal is performed by 1 achieving a space rendezvous between a moving space vehicle and each debris followed by a soft capture using a robotic arm. The shuttle has to meet each debris on its orbit until all debris have been dealt with and then return to its initial orbit. This tour has to be achieved in the less expensive and the fastest possible way. Most collisions are not debris/satellite but rather debris/debris and result in an increasing number of space debris [START_REF] Kessler | Collisional cascading: The limits of population growth in low earth orbit[END_REF][START_REF] Anselmo | Updated results on the long-term evolution of the space debris environment[END_REF]. So the earlier the removal is finished, the less new debris are generated.

In this article, we study the problem of removing a list of space debris and extend the results obtained in [START_REF] Madakat | Biobjective planning of an active debris removal mission[END_REF]. We propose an exact algorithm based on a branch and bound procedure to compute the set of non dominated (cost, duration) vectors and give for each of these vectors a feasible solution. Two specificities that make our model more realistic are taken into account: between each pair of targets, several transfer possibilities, with different costs and times, are considered. These costs and times depend on the start time from the initial target in the pair. Therefore, the problem is modeled as a biobjective time dependent traveling salesman problem defined on a multigraph. To the best of our knowledge, BiTDTSP has not been studied before, although several articles study the time dependent traveling salesman problem (TDTSP) and the biobjective traveling salesman problem (BiTSP). TDTSP is a variant of TSP where distances depend on the arrival time to each vertex. Malandraki and Daskin [START_REF] Malandraki | A restricted dynamic programming heuristic algorithm for the time dependent traveling salesman problem[END_REF] described dynamic programming (DP) algorithms for TDTSP extended by Bellman [START_REF] Bellman | Dynamic programming treatment of the travelling salesman problem[END_REF] and finally Held and Karp [START_REF] Held | A dynamic programming approach to sequencing problems[END_REF]. Schneider [START_REF] Schneider | The time-dependent traveling salesman problem[END_REF] also proposed a simulated annealing heuristic to deal with TDTSP. BiTSP has been addressed by several authors. Gendreau et al. [START_REF] Bérubé | An exact-constraint method for bi-objective combinatorial optimization problems: Application to the traveling salesman problem with profits[END_REF] used the -constraint method to efficiently generate the Pareto front of the traveling salesman problem with profits. Schmitz and Niemann [START_REF] Schmitz | A bicriteria traveling salesman problem with sequence priorities[END_REF] were interested in a BiTSP problem motivated by various applications in the context of service delivery in which the second objective relates to priorities among locations to be visited. Paquette and Stutzle [START_REF] Paquete | Design and analysis of stochastic local search for the multiobjective traveling salesman problem[END_REF] analyzed algorithmic components of Stochastic local search (SLS) algorithms for the multiobjective traveling salesman problem. Based on the insights gained, they engineered SLS algorithms for this problem. Lust and Jaszkiewicz [START_REF] Lust | Speed-up techniques for solving largescale biobjective tsp[END_REF] proposed a heuristic resolution based on the two-phase local search method with speed-up techniques for BiTSP.

The paper is organized as follows. In Section 2, we model the problem of removing dangerous space debris. The proposed approach is presented in Section 3. Section 4 is devoted to implementation issues. Computational experiments and results are reported in Section 5 and conclusions are provided in a final section.

Problem description

Context and notations

Given n debris to be removed, the space shuttle has to move from its own orbit and visit the n debris in order to collect them and then return to its first orbit. The total quantity of fuel burned during each transfer from a debris i to a debris j, denoted by c ij , represents the transfer cost. The duration of the transfer is denoted by d ij . The quantity of fuel burned during the mission should not exceed the shuttle capacity, thus the cost cannot exceed a fixed cost c max . Moreover, the duration of the mission should not exceed a fixed duration d max . We assume here that the mission is not carried out by an unmanned space shuttle. To remove each debris, the shuttle has to perform a rendezvous with each debris on its orbit. Thus in the following we associate each debris with its orbit. Costs and durations depend on the way the rendezvous is achieved.

A space rendezvous between a debris and a space shuttle is an orbital maneuver where both arrive at the same orbit and approach to a very close distance. There are several ways to achieve a space rendezvous as shown in [START_REF] Petropoulos | Low-thrust transfers using primer vector theory and a second-order penalty method[END_REF][START_REF] Luo | Interactive optimization approach for optimal impulsive rendezvous using primer vector and evolutionary algorithms[END_REF][START_REF] Arutyunova | Pontryagin's maximum principle for constrained impulsive control problems[END_REF]. In our case, we have chosen to perform a rendezvous following the Lambert method where the shuttle moves between the two orbits undergoing exactly two pulses. Figure 1 shows how the transfer is performed in the Lambert elliptical case. The cost of this transfer is the quantity of fuel burned in order to perform the first and second pulses.

Transfer arc

Initial orbit

Final orbit

The first boost The second boost Figure 1: The transfer orbit of the space shuttle using the elliptical maneuver For each debris i, t i denotes the time at which the shuttle reaches the orbit of debris i. Once the shuttle has reached orbit i, it can immediately start the next transfer to reach another orbit j or wait before beginning the transfer. Indeed, waiting on an orbit may be cheaper and/or quicker to reach the next orbit. The duration of the transfer is the sum of the waiting time in the departure orbit and the travel time to the arrival orbit. In the following we assume that each elementary transfer requires a minimum duration δt trans and a minimum cost δc min . We assume as well that the service time on an orbit i takes a duration of δt serv i .

Formulation of the problem

The studied space objects and the possible transfers between their orbits are modeled by a complete valued digraph G = (V, A) where V represents the set of object's orbits numbered from 0 to n. Vertex 0 represents the shuttle and vertices from 1 to n represent the n debris. The set A corresponds to the set of feasible transfers between orbits. Several arcs may link each pair of orbits depending on the moment on when the shuttle reaches the departure orbit. To each 3-tuple (i, j, t i ), corresponds a set of feasible transfers A ij (t i ). Each element of A ij (t i ) induces a bivalued arc linking i to j whose value is a pair (c ij , d ij ) corresponding to the cost and duration of the travel. As the arcs representing possible transfers depend on the time at which the shuttle reaches the departure orbit, G is a dynamic multigraph. In the figure 2, we have shown the way the digraph is constituted as the shuttle performs its mission.

After it reaches an orbit i, the shuttle must serve it before going on. Thus, the shuttle can leave the orbit i at least at t i d = t i + δt serv i . Due to mission duration constraints, the shuttle should leave the orbit i before an instant limit t i l . When the shuttle leaves i at t i l , it has barely time to achieve the mission before t 0 + d max . Thus, if the shuttle reaches l after visiting V ⊂ V debris, one has

t i l = t 0 + d max -(n + 1 -|V |) × δt trans - i∈V \V δt serv i -δt trans
In the equation above, (n + 1 -|V |) × δt trans is the least time needed to travel to each unvisited debris, i∈V \V δt serv i is the time needed to serve debris and δt trans is required to go back to the first orbit. The duration of time corresponding to possible departure times is discretized using a time step δt. Hence, the transfer can start at any time t i d +w×δt, where w ∈ I(t i ) and I(t i ) determines the set of possible departure times if the shuttle reaches orbit i at t i that is I(t i ) = {0, 1, ...,

t i l -t i d δt }).
For each departure time several transfer durations are possible, the set of possible durations is denoted as D. A departure time t i , w ∈ I(t i ) and transfer duration p ∈ D define a new arc allowing the shuttle to reach j at t j = t i d + d ij (t i , w, p). The valuation of this arc is:

Preprocessing

The cardinality of A ij (t i ) depends on the number of possible departure times from i and on the number of possible durations to perform the transfer from i to j. Therefore A ij (t i ) may contain a very large number of possible transfers. The computation of the cost corresponding to each transfer possibility using Lambert algorithm is very time-demanding (see Section 5.2.2). The algorithm spends much more time computing transfer costs than optimizing the shuttle trajectory. Therefore, optimizing the algorithm performances requires limiting the number of calls to the cost computation function. This is achieved by reducing the cardinality of A ij (t i ) as will be seen in Section 5.2.1.

The branch and bound procedure

In order to enumerate the non dominated vectors for BiTDTSP, we propose a branch and bound enumeration scheme. In our case, each node at level k in the search tree corresponds to a partial solution visiting k debris starting from the shuttle orbit. Its successors are partial solutions visiting k + 1 debris. A partial solution s is characterized by: • V s the set of debris visited in that partial solution.

• l the last visited element.

• T (s) the time spent to visit the set of debris V s .

• C(s) the visit cost incurred to visit the set of debris V s .

Branching

The branching procedure consists in choosing the next node that will be separated and in generating its successors. This requires generating A lj (t l ) for each debris j / ∈ V s . A strategy of selection dictates which of the active nodes in the search tree will be separated next. Different search strategies are adopted in our case. Once a node is chosen, we generate its successors by computing A lj (t l ) for each j / ∈ V s .

Bounding

A bounding procedure allows deleting some partial solutions that may not lead to efficient complete solutions. We present some fathoming rules. Fathoming rules can lead to substantial improvements in algorithmic performance [START_REF] Bazgan | Solving efficiently the 0-1 multi-objective knapsack problem[END_REF]. Given a partial solution s visiting k debris, an extension of s is any feasible solution in which the k first visited debris are identical to s.

In the following, the set of possible extensions of a partial solution s is denoted by Ext(s). In our context we define the two following fathoming rules.

• First fathoming rule D 1

A pool of good candidate solutions is initially generated (see Section 3.3). If the performance of the most optimistic extension of a partial solution is dominated by one of the solutions pool, this partial solution can be discarded. Consider a partial solution s = (V s , l, T (s), C(s)), the performance of any extension of s is at least

T * (s) = T (s) + (|V | -|V s | + 1) × δt trans + i∈V \Vs δt serv i + δt serv l and C * (s) = C(s) + (|V | -|V s | + 1) × δc min .
Thus, if x is a solution from the pool corresponding to a criterion vector (C(x), T (x)) such that:

T (x) ≤ T * (s) C(x) ≤ C * (s)
then each extension of s is in the best case as good as x and can thus be discarded.

• Second fathoming rule D 2

Consider two partial solutions s = (V s , l, T (s), C(s)) and s = (V s , l , T (s ), C(s )) such that

V s ⊇ V s , l = l , T (s) ≤ T (s ) and C(s) ≤ C(s )
The last rule is valid thanks to two assumptions, whose validity will be discussed in Section 5.2.1.

Assumption 1 We assume that, for any orbits i and j, if t i and t i are two possible arrival times to i such that the time difference between these dates is greater than the period of orbit i then each transfer possibility corresponding to the arrival time t i is dominated by at least one transfer possibility corresponding to the arrival time t i .

More precisely:

∀w ∈ I(t i ), ∀p ∈ D t i ≥ t i + P i ⇒ ∃w ∈ I(t i ) ∃p ∈ D such that c ij (t i , w, p) ≤ c ij (t i , w , p ) d ij (t i , w, p) ≤ d ij (t i , w , p
) where P i refers to the period of orbit i. According to this assumption, waiting more than one period on any orbit i before traveling to orbit j does not lead to efficient transfers.

Assumption 2

The second assumption concerns the triangle inequality of travel costs and durations. Given three debris i, j and l, all transfers that allow reaching l from i through j when the shuttle reaches i at t i are dominated by at least one transfer that corresponds to an earlier possible arrival time t i to orbit i and that reaches l directly from i.

This assumption can be formalized as follows:

∀w 1 ∈ I(t i ) ∀(p 1 , p 2 ) ∈ D 2 ∀w 2 ∈ I(t i + δt serv i + d ij (t i , w 1 , p 1 )) t i ≤ t i ⇒ ∃w ∈ I(t i )∃p ∈ T t i + δt serv i + d ik (t i , w, p) ≤ t i + δt serv i + d ij (t i , w 1 , p 1 ) + δt serv j + d jk (t j , w 2 , p 2 ) c ik (t i , w, p) ≤ c ij (t i , w 1 , p 1 ) + c jk (t j , w 2 , p 2 )
Proposition 1 Assumptions 1 and 2 ensure the validity of the fathoming rule D 2 . Proof 1 Let s = (V s , l, T (s), C(s)) and s = (V s , l , T (s ), C(s )) be two partial solutions such that:

V s ⊇ V s , l = l , T (s) ≤ T (s ) et C(s) ≤ C(s ).
Assume that the travel is performed by a shuttle N 1 in the first solution while in the second solution it is achieved by another shuttle N 2.

Case where: V s = V s Let e be an extension of s . If, in the solution e , N 2 leaves l at t such that T (s) + δt serv l ≤ t ≤ T (s) + δt serv l + P l hence N 1 can wait on the orbit l during t -(T (s) + δt serv l )which is possible since this duration is smaller than the period of orbit l.

Knowing that N 1 and N 2 have achieved a rendezvous with debris l and that they are moving at the same speed, at t they will be in the same place.

Hence, it is possible for N 1 to perform the same sequence of transfers that N 2 has achieved and to finish the mission in the same time than N 2 at a lower cost. Thus e is a possible extension for s too.

If N 2 leaves the orbit l at t ≥ T (s) + δt serv l + P l then, according to Assumption 1, there exists a transfer that leaves l at t ∈ [T (s) + δt serv l , T (s) + δt serv l + P l ] and allows N 1 to visit the next debris visited in e earlier and at lower cost. This transfer brings us back to the same initial situation with the only difference that another orbit has been visited. Repeating the same procedure, we construct an extension of s that dominates e .

Case where: V s V s

Let A = V s -V s be the set of debris visited in the solution s and not visited in the solution s . Let us also define B the set of debris that have to be visited in extensions of both s and s . Hence, an extension of s has to visit the set of debris A ∪ B.

Let e be a completion of s that visits the set A∪B in a given order. Thus, the completion has the form {l,

A 1 , B 1 , A 2 , B 2 , .., A h , B h } where A = h i=1 A i and B = h i=1 B i .
The idea is to construct an extension of s from e that costs less and which allows to achieve the mission more quickly than e .

The shuttle N 1 reaches l at t l = t 0 + T (s) while N 2 reaches it at t l = t 0 + T (s ). t l ≤ t l so according to the triangle inequality, there is a sequence of transfers that allow the shuttle N 1 to reach the first element visited in B 1 faster than N 2. Elements of B 1 have to be visited in e and e as well, with a similar reasoning in the case 1, we can extends e in a way that N 1 reaches the last element in B 1 faster than N 2. This brings us back to the initial situation with the difference that we have visited A 1 ∪ B 1 in addition. By repeating this process, we construct an extension e of s that allows N 1 to reach the last element to visit in B h faster than N 2.

This validates the second fathoming rule.

The Branch and Bound algorithm

A solution pool is initially generated so as to use the fathoming rule D 1 as early as possible. For the procedure to be efficient, this pool must be generated quickly while containing reasonably good solutions. For this purpose, we compute costs and durations corresponding to solutions visiting debris in their orbit inclination order which will form the set initializing the solution pool. This order was chosen because it corresponds to some solutions that promote least cost transfers. In fact, transfers that require large changes in inclination cost more than transfers performed between near inclination orbits [START_REF] Wadsley | Optimal visitation order for spacecraft servicing missions[END_REF]. Besides, numerical results show that numbering debris according to their inclination order decreases computation time (see Section 5.2.2). The set of these solutions is stored in an ordered list N D which is updated each time a complete solution is generated.

A working list L contains current partial solutions. Initially L contains one element which is the initial partial solution (∅, 0, 0, 0). As long as there still exist elements in L, the algorithm computes successors of the first element in the list. The separated node is then removed and its successors are inserted in the working list L. Different orders were tested for list L. Once a solution is complete it is inserted into N D and then we remove from N D all dominated solutions. At the end of the algorithm, all elements stored in N D are non dominated. For each partial solution s, Γ s denotes the set of successors of s. if Dominated = false then Insert s i in L unless it is dominated with respect to D 2 Remove from L partial solutions that are dominated by s i Algorithm 1 returns exactly the set of non dominated vectors and provides for each non dominated vector one efficient solution. In procedure 'Update', we first use D 1 to test dominance as it is easier to verify and may cut many branches of the search tree. Notice that the way partial solutions are stored influences the computation time since it can limit the number of dominance tests which will be checked. Three different orders were tested in our case as will be seen in next section.

Implementation issues 4.1 Debris order

The order in which debris are numbered is a crucial implementation issue in branch and bound procedures. We would like to number debris in an order allowing least cost solutions to be generated first. Two elements influence the cost of the transfer between two orbits: the difference between the two orbits inclination and the difference between the two orbits altitudes. In the following, we test two ways of numbering and compare them with a random numbering:

• The orbits inclination numbering.

• The orbits altitude numbering.

Data storage

In the solution pool, complete solutions are inserted in the lexicographic order of the criteria vector: (cost, duration). 1 is defined as follows: Two solutions s 1 = (V 1 , l 1 , C(s 1 ), T (s 1 )) and s 2 = (V 2 , l 2 , C(s 2 ), T (s 2 )) are such that s 1 1 s 2 if and only if:

C(s 1 ) < C(s 2 ) or (C(s 1 ) = C(s 2 ) and T (s 1 ) ≤ T (s 2 ))
The interest of this order is that for each pair of solutions (s 1 , s 2 ) such that s 1 1 s 2 , one has s 2 cannot dominate s 1 . The way a new solution is inserted in the solution pool is shown in Figure 4. For the partial solutions, we have tested the three orders 1 , 2 and 3 . The order 2 is defined for each pair (s 1 , s 2 ) as follows: s 1 2 s 2 if and only if:

|V 1 | > |V 2 | or (|V 1 | = |V 2 | and C(s 1 ) < C(s 2 )) or (|V 1 | = |V 2 | and C(s 1 ) = C(s 2 ) and T (s 1 ) ≤ T (s 2 ))
The interest of this order is that it allows separating least cost partial solutions from those that visit more debris first and decreases the number of dominance tests. The insertion of a new element is done as explained in Figure 5. In the red part of the list, we have to test if the candidate is dominated while in the green one we have to verify that partial solutions that are already inserted are not dominated. The order 3 is defined for each pair (s 1 , s 2 ) as follows:

s 1 3 s 2 if and only if: C(s 1 ) + T (s 1 ) ≤ C(s 2 ) + T (s 2 )
The last order has the advantage to allow less dominance tests to execute without promoting any of the two criteria.

Computational experiments and results

Experimental design

To get the description of space object orbits, we use the most complete source of orbital element information available to the public today which is the two-line elements (TLEs) set given by NORth American aerospace Defense command (NORAD). NORAD distributes orbital information on a very large number of space objects in the form of two-line elements sets (TLEs). Each TLE consists of two lines of 69 characters, containing an object identifier, the epoch and a set of orbital elements that have been estimated by NORAD. The NORAD element sets are "mean" values obtained by removing periodic variations in a particular way [START_REF] Hoots | Models for propagation of NORAD element sets[END_REF]. In order to obtain good predictions, these periodic variations must be reconstructed in exactly the same way they were removed by NORAD . Hence, The NORAD element Table 1: The percentage of trajectories verifying the triangle inequality sets must be used with the propagation algorithm: SGP4. In the present article, we use this algorithm to predict debris positions.

Results

The goals of the experiments are:

• to determine the best order to sort partial solutions.

• to determine the best order for numbering debris in our approach.

• to analyse the impact of using fathoming rules.

• to compare Pareto fronts obtained with different time steps.

• to view some Pareto fronts and compare them with inital solutions pool.

Preprocessing

Verification of the validity of the triangle inequality

We consider five triplets of debris orbits. For each of these triplets (i, j, k), we compute costs and durations of efficient trajectories allowing the shuttle to go from i to k through j given a date t i at which the shuttle reaches orbit i. We compute the same set of trajectories when the shuttle reaches i at t i + 0.25 × P i, t i + 0.5 × P i, t i + 0.75 × P i and t i + P i. The idea is to show that these sets of trajectories are dominated by at least one transfer of A il (t i ). For this purpose, we compute the percentage of each of these sets dominated by at least one transfer among A il (t i ). Results are reported in Table 1. As shown in the table, most of these transfers are dominated which legitimizes the triangle inequality hypothesis.

Determination of transfer parameters

The computation time is directly related to the cardinality of the set of possible transfers. Reducing the computation time requires decreasing this set cardinality. The cardinality of A ij (t i ) depends on:

• The maximum waiting duration allowed on each orbit.

Debris pairs

The % of ND transfers The % of ND transfers The % of ND transfers obtained after waiting obtained after waiting obtained after waiting one period two periods three periods Table 2: The percentage of efficient transfers progression

• The number of possible transfer durations.

• Choice of the maximum transfer duration

One may allow the shuttle to take a long time in order to achieve a transfer, but the loss on time must be compensated by a significant gain on cost. We consider in Figure 6, five pairs of LEO debris and study the evolution of the transfer cost depending on the transfer duration. As shown in the figures, most transfers that take more than 10 days to be performed are dominated. In the following, we assume that an interorbital transfer takes at most 10 days to be achieved.

• Choice of the maximum waiting duration

The departure possibilities number must be limited. For this purpose, we test all the possible waiting durations using an increment of 30 minutes and less than ten times the period of the departure orbit. We then determine the costs of the transfers that depart from each of these positions and take complete days to be achieved without exceeding one month. We retain only non dominated transfers. The percentage of the non dominated transfers already obtained after a number of periods is then computed. We show in Table 2 the evolution of the efficient transfers duration obtained with eight different pairs of orbits.

As shown in the table, the majority of non dominated transfers is already obtained after testing departure positions allowed by waiting one period. Therefore, in the following, we assume that the waiting time in any orbit cannot exceed the period of this orbit.

The impact of partial solutions sorting order and of debris numbering

We compare here the performance the algorithm when partial solutions are inserted in the three orders and when we number debris in their orbit in- As shown in Table 3, the best order to sort partial solutions is 2 which is the lexicographic order of (|V |, cost, duration). Furthermore the best way of numbering the orbits is the inclination order. The lexicographic order promotes partial solutions visiting more debris which accelerates the construction of complete solutions which improves the solutions pool and strengthen the first fathoming rule. Furthermore numbering debris in their inclination order promotes least cost solutions to appear first.

Impact of fathoming rules

We use here the 2 to sort partial solutions, and we number debris in their inclination order. To view the impact of fathoming rules, we compare the computation time with and without each fathoming rule. As shown in Table 4 the use of fathoming rules decrease notably the computation time in the two cases. The best computation time is obtained when both of the fathoming rules are used.

Impact of time step

In order to view the impact of the time step, we insert partial solutions in the order 2 , use both of dominance relations and number debris in their orbit inclination order. For several instance sizes, we generate the set of non dominated solutions obtained with each instance and compute the number of non dominated vectors, and the corresponding computation time. We increase the time step and compute the set of non dominated vectors obtained using greater time steps. The idea is to compare the quality of these sets of non dominated debris. We use the -approximation relation introduced in [START_REF] Zitzler | Performance assessment of multiobjective optimizers: An analysis and review[END_REF]: for a given accuracy > 0, an (1 + )-approximation is a subset of solutions which contains, for each efficient solution, a solution that is at most at a factor (1 + ) on both objective values. Thus, to evaluate the quality of the Pareto front obtained with greater time steps, we compute the least value of such that the second Pareto front is an -approximation of the first one. As shown in Table 5, the use of greater time steps decreases the computation time but deteriorates notably the quality of non dominated set.

Studying some Pareto fronts

In this paragraph, we propose to view some Pareto fronts obtained with several instance sizes. In each case we compare the quality of the non dominated vectors set and the solutions pool used in the beginning of the algorithm. To compare these two sets quality a worst case criterion is used. Hence we compute the least value of such that the solutions pool is an -approximation of the final set of non dominated vectors.

In Figure 7, we present the set of non dominated vectors obtained with a six debris instance (solutions marked with '+') and performances of the pool of solutions initially generated for the same instance (solutions marked with points). As can be seen in the figures, solutions that visit debris according to their inclination order present a good approximation of the set of non dominated vectors. Hence, the least value of the for which the solutions pool is an -approximation of the final set of non dominated vectors is 0.49 for the 6 debris instance, 0.33 for the seven's 0.22 for the eight's and 0.13 for the nine debris instance. In addition, the order of the visit in the set of non dominated vectors obtained in each case is not very different from the inclination order. Thus, most of the non dominated vectors differ from the set of solutions pool by less than three swaps.

Conclusions

The purpose of this work was to propose a new modeling of the problem of removing potentially dangerous space debris. The problem was presented as a time dependent biobjective traveling salesman problem. The BiTDTSP has not been studied in the literature before.

To efficiently solve the BiTDTSP, we proposed a branch and bound approach. The use of dominance relations avoids exploring some partial solutions that would not lead to non dominated vectors, and the limitation of the 'on orbit waiting durations' gives a good approximation of the set of non dominated vectors in less time. As a perspective of this work, we will study more deeply the interesting area in each non dominated solution set by using smaller time steps. As removing all space debris is too expensive another perspective of this work is to model and efficiently solve the problem of determining the most troublesome debris.

Figure 2 :

 2 Figure 2: Progression of the graph of the studied objects

Figure 3 :

 3 Figure 3: Transfer possibilities when the shuttle reaches orbit i at t i

Algorithm 1 :j ← j + 1 if 1 elseAlgorithm 2 :

 1112 Branch and bound algorithm Output: Set of non dominated vectors N D begin L ← {(∅, 0, 0, 0)} Initialize N D a set of candidate non dominated vectors while L = ∅ do /*assume N D = {p 1 , p 2 , ..., p r } where solutions are ordered by non increasing costs*/ Select s ∈ L L ← L\{s} if s is complete then Dominated ← f alse j ← 1 while j ≤ r and C(p j ) ≤ C(s)and Dominated = f alse do if T (p j ) ≤ T (s) then Dominated ← true if C(s) = C(p j ) and T (s) ≤ T (p j ) then N D ← N D\{p j } Dominated = f alse then N D ← N D ∪ {s} while j ≤ r and T (s) ≤ T (p j ) do N D ← N D\{p j } j ← j + Generate Γ s Update(L, Γ s , N D) return N D Procedure Update (L, L , N D) Input: Two lists of partial solutions L and L and list of pool solutions N D Output: Working list L updated begin /*assume L = {s 1 , s 2 , ..., s q } and N D = {p 1 , p 2 , ..., p s }*/ for i ← 1 to q do Dominated ← f alse j ← 1 while Dominated = false and j ≤ s do if p j D 1 s i then Dominated ← true j ← j + 1

Figure 4 :

 4 Figure 4: The insertion of a new solution in the solution pool

Figure 5 :

 5 Figure 5: The insertion of a new partial solution with respect to 2 order

Figure 6 :

 6 Figure 6: Changes in transfer cost depending on travel duration

Figure 7 :

 7 Figure 7: Non dominated vectors and pool solutions for different instances of debris

  (c ij (t i , w, p), d ij (t i , w, p)).

	The transfer possibilities can be seen in Figure 3.
	t	i	+	t δ +	t δ	trans	i t δ t +	trans				t t i δ +
													t	i
													t	i	+	t δ	+	3	t δ .	trans
								t	i	+	trans i t t δ . 2 δ t +	+	t δ . 2	trans	t	i	+	t δ . 3	trans

Table 3 :

 3 The obtained computation time using different orders

	Instance	With D 1	With D 2	With D = D 1 ∪ D 2
	size	B & B	global	B & B	global	B & B	global
		time	time	time	time	time	time
	3	3	32	3	35	2	26
	4	3	145	3	159	2	132
	5	5	756	6	680	5	630
	6	23	1373	22	1327	18	1250
	7	67	5643	62	5472	52	4199
	8	156	8340	134	7983	105	7600
	9	463	10376	423	10562	371	9341
	10	714	11540	702	11582	658	10321

Table 4 :

 4 The impact of fathoming rules use on computation time clination order or in their orbit altitude order. The time step is set to 15 minutes.

Table 5 :

 5 The impact of time step on computation time

* This research has been partially supported by the project ANR-09-BLAN-0361 ' GUaranted Efficiency for PAReto optimal solutions Determination'(GUEPARD).


Agnese Maria

Luca Levrini

Laura Rampazzi

Carlo Dossi

Lanfredo Castelletti

Marta Licata

Cristina Corti
email: cristina.corti@uninsubria.it

Barbara Lazzati A

Laura Rampazzi B

Carlo Dossi C

Marta Licata E

The diet of three medieval individuals from Caravate (Varese, Italy). Combined results of ICP-MS analysis of trace elements and phytolith analysis conducted on their dental calculus

Keywords: Medieval diet, dental calculus, phytolith, optical microscopy, SEM-EDX, trace elements, ICP-MS

de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

The Diet of Three Medieval Individuals from Caravate (Varese, Italy). Combined Results of ICP-MS Analysis of Trace Elements and Phytolith Analysis Conducted on Their Dental Calculus

INTRODUCTION

Paleodiet is defined as the field of study of anthropology that, through the examination of a variety of indicators (skeleton and teeth, dental calculus, hair, nails, gall stones, calcified tumours, soft tissue residuals, etc.), collectable from human remains found in archaeological sites, attempts to reconstruct the diet of a specific past population.

Teeth are in close contact with foodstuff that is introduced with the diet. Food and oral health are thus two aspects that are closely related. Food has a direct role on the balance of oral flora, which, if altered, may lead to phenomena such as caries, formation of dental calculus, and periodontal disease. The dental calculus, or tartar, is a hard concretion, which is formed by calcification of plaque. It is located on the surfaces of the teeth, on both (buccal and lingual) sides of the teeth, in a supra-or sub-gingival position.

The presence of pathogenic bacteria and tartar in the oral cavity can lead to cavitary lesions, of varying depth, or loss of the element itself intra vitam, due to resorption of the supporting tissues such as bone and periodontal ligament. Bacteria can easily remain trapped in the crystalline matrix, so dental calculus can be considered as a precious source of preserved human microbiota [START_REF] Adler | Sequencing ancient calcified dental plaque shows changes in oral microbiota with dietary shifts of the Neolithic and Industrial revolutions[END_REF]. Genetic analyses on ancient bacterial DNA can be used to trace dietary changes, health and movements of past populations. All these elements are important in understanding the type of diet, the processing of food and its cooking, as well as the general state of health of the individual.

Food can also affect the degree of wear of the tooth crowns. Therefore, in addition to a macroscopic study, a micro-configuration analysis of wear allows us to understand the type of diet (hard or soft) and whether based mainly on vegetables, grains or meat [START_REF] Lalueza | Dietary information through the examination of plant phytoliths on the enamel surface of human dentition[END_REF][START_REF] Ungar | Dental microwear and diets of African early Homo[END_REF].

Dental calculus is constituted by both organic and inorganic components. The organic fraction is made of proteins, carbohydrates and lipids and amounts to about 12% of the total weight; the remaining 88% has an inorganic nature, and it is made of calcium phosphate, traces of magnesium phosphate and other elements. The composition is very similar to that of bones and teeth [START_REF] Henry | Microfossils in calculus demonstrate consumption of plants and cooked foods in Neanderthal diets (Shanidar III, Iraq; Spy I and II, Belgium)[END_REF][START_REF] Lieverse | Diet and the aetiology of dental calculus[END_REF][START_REF] Lindhe | Clinical periodontology and implant dentistry[END_REF][START_REF] Perez | Microscopic X-ray fluorescence analysis of human dental calculus using synchrotron radiation[END_REF]. Its mechanism of formation is still not fully understood, but it surely involves the deposition of calcium and phosphate from the saliva and gingival fluids into the matrix of biofilm that constitutes oral plaque [START_REF] Hayashizaki | Site specific mineral composition and microstructure of human supra-gingival dental calculus[END_REF]. Many authors linked the formation of calculus with the amount of carbohydrates consumed, whilst others showed a connection with the ingestion of proteins. The controversy is still open [START_REF] Arnay-De-La-Rosa | The Ba/Sr ratio, carious lesions, and dental calculus among the population buried in the church La Concepcion (Tenerife, Canary Islands)[END_REF]. For what concerns the time frame of calculus accumulation, the rate of deposition depends on a number of aspects, such as genetic predisposition, diet and oral healthcare [START_REF] Mickleburgh | New insights into the consumption of maize and other food plants in the pre-Columbian Caribbean from starch grains trapped in human dental calculus[END_REF].

Dental calculus is mineralized, so it commonly survives in archaeological human remains, in the form of a hard deposit on teeth [START_REF] Lieverse | Diet and the aetiology of dental calculus[END_REF]. Different calcium phosphates can be found, partly depending on the degree of 'maturation' of the pathological concretion. In particular, hydroxyapatite, dicalcium phosphate dihydrate, octacalcium phosphate and whitlockite are generally observed [START_REF] Hayashizaki | Site specific mineral composition and microstructure of human supra-gingival dental calculus[END_REF][START_REF] Perez | Microscopic X-ray fluorescence analysis of human dental calculus using synchrotron radiation[END_REF]. The crystals of phosphates can also exchange ions with the oral fluids and adsorb organic compounds [START_REF] Abraham | A case study of elemental and structural composition of dental calculus during several stages of maturation using SRXRF[END_REF][START_REF] Perez | Microscopic X-ray fluorescence analysis of human dental calculus using synchrotron radiation[END_REF]. Divalent cations can enter the lattice and replace calcium and many different anions, such as silicate, can substitute phosphates.

Whilst the major elements in dental calculus are Ca, P and Mg, its minor components are very variable and not still extensively studied [START_REF] Hayashizaki | Site specific mineral composition and microstructure of human supra-gingival dental calculus[END_REF][START_REF] Molokhia | Studies on the composition of human dental calculus. Determination of some major and trace elements by Instrumental Neutron Activation Analysis[END_REF][START_REF] Perez | Microscopic X-ray fluorescence analysis of human dental calculus using synchrotron radiation[END_REF].

Some articles about chemical analysis of calculus can be found in the literature, generally focused on mineralogical phases of modern individuals, but sometimes also on elemental composition and on ancient samples. It is reported, for example, the use of X-ray diffraction (XRD) [START_REF] Hayashizaki | Site specific mineral composition and microstructure of human supra-gingival dental calculus[END_REF], Fourier transform infrared spectroscopy (FTIR) [START_REF] Hayashizaki | Site specific mineral composition and microstructure of human supra-gingival dental calculus[END_REF], fluorescent X-ray spectroscopy (XRF) [START_REF] Hayashizaki | Site specific mineral composition and microstructure of human supra-gingival dental calculus[END_REF] and synchrotron radiation analysis by X-ray fluorescence (SRXRF) [START_REF] Abraham | A case study of elemental and structural composition of dental calculus during several stages of maturation using SRXRF[END_REF][START_REF] Perez | Microscopic X-ray fluorescence analysis of human dental calculus using synchrotron radiation[END_REF][START_REF] Sanchez | SRXRF analysis with spatial resolution of dental calculus[END_REF]. Data were also collected using scanning electron microscopy associated with X-ray spectrometry (SEM-EDX) [START_REF] Charlier | The microscopic (optical and SEM) examination of dental calculus deposits (DCD). Potential interest in forensic anthropology of a bio-archaeological method[END_REF][START_REF] Hardy | Neanderthal medics? Evidence for food, cooking, and medicinal plants entrapped in dental calculus[END_REF], gas chromatographymass spectrometry [START_REF] Hardy | Neanderthal medics? Evidence for food, cooking, and medicinal plants entrapped in dental calculus[END_REF] and instrumental neutron activation analysis (INAA) [START_REF] Molokhia | Studies on the composition of human dental calculus. Determination of some major and trace elements by Instrumental Neutron Activation Analysis[END_REF].

As previously ascertained, dental calculus can be a source of information about diet, through the analysis of chemical trace elements. In recent years, considerable progress has been made in dietary studies, thanks to the use of instruments such as Inductively Coupled Plasma Mass Spectrometry (ICP-MS), by which it is possible to perform a fast and precise analysis of trace elements contained in various types of samples, including teeth and bones. The evaluation of concentration and the isotopic abundances of certain elements to obtain information about past populations is now possible [START_REF] Budd | The distribution of lead within ancient and modern human teeth: Implications for long-term and historical exposure monitoring[END_REF][START_REF] Busetto | Dietary investigation by trace element content in bones of ancient inhabitants of Northern Italy[END_REF][START_REF] Corti | On the use of trace elements in ancient necropolis studies: Overview and ICP-MS application to the case study of Valdaro site, Italy[END_REF][START_REF] Djingova | On the possibilities of inductively coupled plasma mass spectrometry for analysis of archaeological bones for reconstruction of paleodiet[END_REF][START_REF] Klepinger | Nutritional assessment from bone[END_REF][START_REF] Lambert | Chemical analysis of excavated human bone from Middle and Late Woodland sites[END_REF][START_REF] Schutkowski | Diet, status and decomposition at Weingarten: trace element and isotope analyses on Early Mediaeval skeletal material[END_REF][START_REF] Simonetti | In situ elemental and Sr isotope investigation of human tooth enamel by Laser Ablation-(MC)-ICP-MS: successes and pitfalls[END_REF][START_REF] Szostek | Trace element analysis of Bronze Age skeletal and crematory graves from Southern Poland for diet reconstruction[END_REF][START_REF] Szostek | Trace elements concentrations in human teeth from a Neolithic common grave at Nakonowo (Central Poland)[END_REF][START_REF] Tomczyk | Dental caries and chemical analyses in reconstruction of diet, health and hygienic behaviour in the Middle Euphrates valley (Syria)[END_REF][START_REF] Webb | Inductively coupled plasma-mass (ICP-MS) and atomic emission spectrometry (ICP-AES): Versatile analytical techniques to identify the archived elemental information in human teeth[END_REF].

The most important elements used in palaeonutritional studies are strontium and zinc. Strontium is contained in foods such as cereals and vegetables and can provide an indication of the consumption of carbohydrates, while zinc comes from foods such as red meat, fish, egg yolk, dairy products, vegetables and seafood, and indicates the consumption of proteins [START_REF] Beck | Bivariate analysis of trace elements in bone[END_REF][START_REF] Busetto | Dietary investigation by trace element content in bones of ancient inhabitants of Northern Italy[END_REF][START_REF] Janos | Evaluation of elemental status of ancient human bone samples from Northeastern Hungary dated to the 10th century AD by XRF[END_REF][START_REF] Lambert | Ancient human diet from inorganic analysis of bone[END_REF][START_REF] Rheingold | Strontium and zinc content in bones as an indication of diet[END_REF]. Lambert et al. reported values of strontium in bones in the range 400-500ppm for herbivores, 150-400ppm for omnivores and 100-300ppm for carnivores; they also reported low concentrations of zinc in herbivorous animals (90-150ppm), intermediate for omnivores (120-220ppm) and higher for carnivores (175-250ppm) [START_REF] Lambert | Ancient human diet from inorganic analysis of bone[END_REF]. The reason why the amount of these elements in bones can change is linked to the element's capability to replace small percentages of calcium in normal processes of resorption and redeposition, which is typical of bone tissue and depends on the type of diet. Lead is another important element in this analysis, due to the fact that its value increases in the presence of post-mortem contamination, so it can be used as a marker of the 'good preservation' of skeletal samples. Finally, isotopic analysis of carbon and nitrogen has been applied to tooth and bone collagen to indicate the source of certain vegetable foods and to determine the prevalence of a vegetarian or carnivorous diet.

Similar to bones, dental calculus (if not removed) is in a continuous state of formation and remodelling and can potentially record elemental information during life, whilst enamel tends to be more resistant to the incorporation of elements, due to its peculiar stability after formation [START_REF] Budd | The distribution of lead within ancient and modern human teeth: Implications for long-term and historical exposure monitoring[END_REF][START_REF] Webb | Inductively coupled plasma-mass (ICP-MS) and atomic emission spectrometry (ICP-AES): Versatile analytical techniques to identify the archived elemental information in human teeth[END_REF]. In agreement with the approach of [START_REF] Scott | Stable carbon and nitrogen isotopes of human dental calculus: a potentially new non-destructive proxy for paleodietary analysis[END_REF] dental calculus can be a suitable material to investigate paleodiet in such cases when museum curators do not authorize the use of archaeological teeth and bones. Dental calculus is not an inherent original part of the skeleton, but just an 'add-on' so its sampling and analysis can be considered as a 'non-destructive' technique. [START_REF] Scott | Stable carbon and nitrogen isotopes of human dental calculus: a potentially new non-destructive proxy for paleodietary analysis[END_REF] also found encouraging clues to support the link between diet and dental calculus stable isotope composition, which were a lead for our study.

Moreover, further features of the paleodiet can be revealed by the analysis of phytoliths [START_REF] Blatt | Dirty teeth and ancient trade: Evidence of cotton fibres in human dental calculus from Late Woodland, Ohio[END_REF][START_REF] Buchet | Comparison between the distribution of dental microstriations and plant material included in the calculus of human teeth[END_REF][START_REF] Dudgeon | Diet, geography and drinking water in Polynesia: microfossil research from archaeological human dental calculus, Rapa Nui (Easter Island)[END_REF][START_REF] Henry | Microfossils in calculus demonstrate consumption of plants and cooked foods in Neanderthal diets (Shanidar III, Iraq; Spy I and II, Belgium)[END_REF][START_REF] Henry | Using plant microfossils from dental calculus to recover human diet: a case study from Tell al-Raqa'i, Syria[END_REF][START_REF] Lalueza | Phytolith analysis on dental calculus, enamel surface, and burial soil: information about diet and paleoenvironment[END_REF][START_REF] Menendez | A first phytolith approximation to diet study in the archaeological site Gascon 1, Pampean Region, Republica Argentina[END_REF][START_REF] Piperno | Phytoliths: a comprehensive guide for archaeologists and paleoecologists[END_REF][START_REF] Power | Assessing use and suitability of scanning electron microscopy in the analysis of microremains in dental calculus[END_REF][START_REF] Wesolowski | Evaluating microfossil content of dental calculus from Brazilian sambaquis[END_REF]. Phytoliths are granules of amorphous silica (95%), water (4%) and trace elements (1%), of variable size between 5 and 200µm, produced by plants. They are deposited in the oral cavity while chewing fruit and vegetables and then incorporated into dental calculus during its formation. Their morphologies are specific to particular plants, therefore, from the recovery of an isolated phytolith, it is possible to recognize taxa of higher rank, such as families, tribes or subfamilies of the plant that formed it, thus obtaining information about the diet. In particular, the phytoliths research and identification point to the food eaten in the last times before death, are nondestructive, cheap and give information on specific plants (Boyadjian et al., 2007). Hardy et al. [START_REF] Hardy | Neanderthal medics? Evidence for food, cooking, and medicinal plants entrapped in dental calculus[END_REF] considered both phytoliths and organic chemical compounds found in dental calculus of a Neanderthal population. Carbohydrate markers, fatty acid esters and combustion markers were directly linked to the ingested food and to the method of cooking. A similar approach was used by Buckley et al. [START_REF] Buckley | Dental calculus reveals unique insights into food items, cooking and plant processing in prehistoric Central Sudan[END_REF].

Therefore, the aim of our work was to set up an analytical protocol, for the investigation of phytoliths and the elemental composition of dental calculus. To the best of our knowledge, this is the first use of ICP-MS for trace element analysis on dental calculus. The research was conducted on samples coming from medieval skeletons that were excavated in Caravate (Varese, Italy).

ARCHAEOLOGICAL CONTEXT This study was conducted by examining skeletons that were buried adjacent to the church of St. Augustine in Caravate (Figure 1), an Italian town of 2632 inhabitants near Varese, Lombardy.

The origin of the site name 'Caravate' is uncertain: it is thought to derive either from the Celtic name of person 'Carevus', or from pre-Latin carabus or caravos, which derives from karra, i.e. stone (or from caravum, heap of stones). It is also possible that this name derives from the ancient Caravè, a place where hemp was macerated [START_REF] Armocida | Il nostro paese Caravate. Viganò Grafiche[END_REF]. The latter hypothesis is supported by some findings of fossils, which suggest that the area surrounding this town was originally a large swamp in which hemp was placed to soak. The history of this town extends from the Celtic domination up to the Roman. The reason for so much attention on the site by so many people probably lies on its strategic position, which made it easy to defend. Although Christianity was practiced in Caravate since the 4th century AD, it was not until the subsequent conquest of the village by the Lombards that the building of churches was promoted, like the one dedicated to St. Augustine.

This church was built between the 11th and 12th centuries AD and was later renovated several times. During some recent restoration work a subsidence of the land suggested the presence of a cemetery on the north side of the building. Since then an archaeological survey was carried out in four phases, from August-October 2002 and again in January-March 2003 under the direction of Dr. M. A. Binaghi. After a simple clean up of the site and the levelling of the land, the first two burials came to light in a tomb structured and covered with a stone slab.

These burials were aligned with each other and oriented along the east-west axis in the Romanesque tradition, about 2m away from the wall of the church. After removing the cover stone, two skeletons of adult individuals in good condition were found. A third burial was revealed, leaning against the north wall of the church. At the time of the discovery, the tomb, whose small size suggested it was built for a child, was already emptied and its structure overturned. A fourth burial was later found but was full of debris and, unlike the others, the tomb was built of squared limestone blocks but did not have any coverage.

The clean up of the entire area north of the church of St. Augustine was then completed and other graves and niches were discovered, aligned along the east-west axis, but built in different times.

Among all the tombs discovered, those against the wall of the church were the oldest, whilst those farther away from the church were more recent. The different periods of construction were also evident based on the fact that the oldest tombs were well structured, whilst the subsequent ones showed the re-use of previous tombs. This was confirmed by the discovery of a skeleton of an adult in a tomb originally belonging to a child: to be able to push the body within the available space, the body was placed on its side and several bones in the neck, arms and legs were broken. To expand the small tomb as much as possible the closure to the east was also torn down.

Besides the re-use of previously existing graves, materials from the tombs were later also recycled to satisfy the demand of a growing population. This was easily visible in burial number 5, which was partially superimposed on the oldest loculus and some of the original tomb components were used.

At the end of the survey in the areas north and east of the site, ten tombs had been found, seven still sealed and three emptied. The seven skeletons were in good condition and belonged both to adults and to adolescents. In addition to these, a square stone bearing an inscription was found and is currently under study.

MATERIALS AND METHODS

Skeletons and samples of dental calculus

A total of seven samples were collected from three skeletons, named respectively Tb6 US 115, Tb7 US 118 and Tb9 US 129.

They were the only ones available because the rest had already been cleaned, catalogued and studied by archaeologists. All the dental calculus found on the teeth of these three skeletons was collected, subdividing the various samples depending on the tooth and the surface where calculus was collected. Five samples were taken from the lingual surface of the teeth, which faces the inside of the oral cavity and the tongue, and two from the buccal surface, which faces the lips.

Two of the samples were taken from the first skeleton Tb6 US 115: one from the lingual surface of the element 3.4 (permanent mandibular left first premolar) and one from the lingual surface of the element 3.5 (permanent mandibular left second premolar).

From the second skeleton, one sample was collected from the lingual surface of the element 4.6 (permanent mandibular right first molar).

Eventually four samples were collected from the third skeleton: one from the buccal surface of the element 3.3

(permanent mandibular left canine), one from the buccal surface of the element 4.4 (permanent mandibular right first premolar), one from the lingual surface of the element 4.6 (permanent first mandibular right molar) and one from the lingual surface of the element 4.7 (permanent mandibular right second molar).

Collected samples are summarised in

Reagents

Ultrapure water (18MΩ resistivity, <5ppb TOC) from a Milli-Q Gradient A10 system (Millipore) was used throughout the analytical procedure. Ultrapure HCl and HNO 3 were produced by sub-boiling distillation in a quartz apparatus (DuoPUR, Milestone). Glacial CH 3 COOH (Suprapur grade from Merck) was used to clean the samples. A 10mg/l multi-element standard solution (Merck Certipur ICP multi-element standard solution XXI)

was used to prepare standard solutions at lower concentrations. Another solution (Multi-element Tune A solution, Analytika) was used to optimise the ICP-MS settings.

All the manipulations were performed in a Class 100 laminar-flow hood to prevent external contamination.

Low-density polyethylene (LDPE) sample and reagent bottles and vials were washed and stored in 2% HNO 3 and rinsed with ultrapure water before use.

Sampling and sample preparation

A modified version of the protocol of C. Lalueza was followed [START_REF] Lalueza | Phytolith analysis on dental calculus, enamel surface, and burial soil: information about diet and paleoenvironment[END_REF]. The teeth were cleaned with a soft toothbrush and abundant deionised water in order to eliminate possible soil residues. Small fragments of dental calculus were then removed with a sterile curette and kept in clean LDPE vials.

Optical microscopy and SEM-EDX for phytolith analysis

The samples were soaked in a solution of 20% HCl for 6-8 hours and then disintegrated with the aid of ultrasonic waves. The residuals later underwent three series of washings with ultrapure MilliQ water, centrifugation and removal of the liquid part of the samples. This would eliminate most of the chlorides, which otherwise would have had the tendency to crystallise.

The samples were then transferred into plastic coverslips and left to dry in a stove overnight at 70°C. Once dried, they were observed through optical and electronic microscopes.

An optical Nikon Eclipse E600 microscope was used, equipped with a digital acquisition system (Nikon Digital Sight D5-U1) and controller software (Nikon NIS-Elements F).

A field emission gun scanning electron microscope (FEG-SEM) FEI XL-30 was used along with an Everhart-Thornley detector for secondary electrons (SE), annular detector for backscattered electrons (BSE) and X-ray spectrometer (EDAX) detector EDAX Sirion 200/400 with controller software EDAX Genesis 5.21. Prior to the analysis, samples were thinly coated with gold or carbon.

Elemental analysis of the samples was performed with EDAX and each chemical element of interest was matched to a different colour in order to obtain colour maps useful for displaying the phytoliths. When a silica particle was found, it was enlarged, photographed and its morphology was observed. Phytoliths were then identified with the help of a collection of recent specimens existing at the 'Laboratorio di Archeobiologia' in Como.

ICP-MS for trace elements analysis

Samples were taken with PTFE-coated stainless steel tweezers in order to prevent contamination, and then rinsed with a 5% CH 3 COOH solution to eliminate the superficial layers that could have been contaminated by soil during burial or by metallic tools during excavation or sampling [START_REF] Bentley | Human mobility at the Early Neolithic settlement of Vaihingen, Germany: evidence from strontium isotope analysis[END_REF][START_REF] Price | Immigration and the ancient city of Teotihuacan in Mexico: a study using strontium isotope ratios in human bone and teeth[END_REF].

Samples were then rinsed with ultrapure water and desiccated in a drying oven at 40°C for 3 hours.

Samples were weighed (Tb6 US 115 8.17mg; Tb7 US 118 9.13mg) and placed in clean LDPE bottles: then 1ml of ultrapure HNO 3 67.4% was added. Bottles were closed and placed in an ultrasonic bath for 30s to enhance the dissolution. After this step, the resulting samples were completely disaggregated.

The solutions obtained were then gravimetrically diluted to 3g with ultrapure water and stored at 4°C until analysis.

The solutions were then diluted at 1:200 and at 1:15 with a solution of 2% HNO 3 in ultrapure water just before the analysis to assure that the analytes were within the calibration range. Two final solutions were obtained for every sample.

A procedural blank solution was prepared with ultrapure 67.4% HNO 3 and ultrapure water.

An ICP-MS (Thermo Elemental X-Series II), equipped with a Cetac ASX-260 autosampler and a concentric nebuliser was used. Optimisation of the instrument settings was performed daily with a 10µg/l multi-standard solution, as recommended by the manufacturer.

The optimised parameters used for analysis are: argon plasma gas flow 13l/min; auxiliary gas flow 0.8l/min; nebuliser gas flow 0.90l/min; ICP-RF power 1400W; dwell time 10,000ms; acquisition mode peak jumping;

sample uptake rate 0.4ml/min; sweeps per reading 80; sampling depth 150 arbitrary units; instrumental repetitions 5 runs.

Multistandard solutions at 0.05, 0.1, 0.5, 1, 5, 10, 20 and 50µg/l in 2% HNO 3 were prepared. The following isotopes were chosen: 7 Li, 9 Be, 23 Na, 25 Mg, 27 Al, 39 K, 51 V, 52 Cr, 55 Mn, 59 Co, 60 Ni, 65 Cu, 66 Zn, 69 Ga, 75 As, 85 Rb, 88 Sr, 111 Cd, 137 Ba, 208 Pb, 209 Bi, 238 U.

A procedural blank was also analysed.

RESULTS AND DISCUSSION

Phytolith analysis by Optical Microscopy and SEM-EDX

Two samples were taken: one from the lingual or inner surface of the mandibular second right premolar (element 3.5) of the skeleton Tb6 US 115 and another from the buccal or outer surface of the mandibular right canine (element 3.3) of the skeleton Tb9 US 129. Both samples were analysed with EDAX.

The results of this analysis showed that the elemental spectra of these two samples were similar: all spectra were dominated by a peak of chlorine, caused by the use of hydrochloric acid to dissolve the calculus, and in all of them there were peaks of other elements such as potassium, sulphur, phosphorus, carbon, nitrogen and oxygen, whose origin could be either organic or inorganic. Also, aluminum and silicon were detected, whose origin is definitely exogenous. Due to the high content of silicon in the samples, it was possible to hypothesize the presence of phytoliths, which consist predominantly of this element.

For what concerns the research of phytoliths, it must be noted that their probability to remain trapped inside the oral cavity during mastication of fruit and vegetables is rather poor and that the same calculus formation takes place in relatively short times (days or weeks).

Thanks to EDAX it was possible to create maps associating specific chemical elements to different colours: red has been assigned to carbon, green to aluminium and blue to silicon.

The maps created with EDAX made it possible to identify the particles containing silicon, because they were blue: once photographed, their morphology was studied by means of optical and electronic microscopy. An example of a map can be observed in Figure 2.

The observation of the samples with optical and electronic microscopy led us to identify some phytoliths from the Poaceae family (also named Gramineae) (Figure 3), with elongated, polyhedral or thorn-like shapes.

Phytoliths are present in the cells of the chaff, of the rachis of the ear, and in the stem and leaves. When using naked grain cereals or those that have been hulled but thoroughly cleansed from their glumes (husked), there is little chance that the phytoliths would be introduced with food in the mouth. The Gramineae phytoliths appear to be scarce and difficult to identify due to the lack of characteristic features for identification. Most of these phytoliths were fragmented and occurred very sporadically. What is preserved in the dental calculus is often not part of the actual food but is instead sporadically deposited from an external source, or perhaps through food preparation and manipulation before, during and after cooking.

For example, it is in theory possible to find fragments of wood from wooden containers or cutlery, chips of charcoal from the fireplace (micro-charcoals), inorganic substances such as mineral silica, or tectosilicates from a totally exogenous source, among other things [START_REF] Buckley | Dental calculus reveals unique insights into food items, cooking and plant processing in prehistoric Central Sudan[END_REF][START_REF] Mickleburgh | New insights into the consumption of maize and other food plants in the pre-Columbian Caribbean from starch grains trapped in human dental calculus[END_REF].

Some phytoliths from wood of the class of Dicotyledons and perhaps Conifers (e.g. for the symbolism of conifers Rabanus Maurus [START_REF] Petrocchi | Il simbolismo delle piante in Rabano Mauro e altri studi di storia medievale[END_REF]) were also found: some smooth rods, some polyhedral fragments (Figure 4) and one small probable fibre. Obviously, these are part of unintentionally ingested elements, received with food or directly into the mouth by other means.

ICP-MS trace element analysis

The procedure set up for ICP-MS analysis required a sample to be dissolved, amounting to at least 5mg. Only the Tb6 US 115 and Tb7 US 118 samples satisfied this requirement as enough sample was saved after the treatment for the research of phytoliths.

The analytical results are summarised in Table 2 where the concentrations of all the analysed elements are reported. In the following discussion only the elements which are relevant from the dietary point of view will be considered.

Due to the lack of reports in the literature about using ICP-MS on fragments of dental calculus, it was decided to follow a sample pre-treatment similar to the one used on tooth enamel and bone, because they are all constituted predominantly of calcium compounds, independently from their mechanism of formation [START_REF] Lindhe | Clinical periodontology and implant dentistry[END_REF].

Both samples (Tb6 US 115 and Tb7 US 118) gave similar results in terms of elemental composition. Their concentrations are comparable to those that would be expected from an analysis of dental enamel or dental calculus. The two biological materials have different mechanisms of incorporation of trace elements, the first one is biological, whilst the second is a physical entrapment and can be seen as a stochastic process. An approach similar to [START_REF] Hardy | Neanderthal medics? Evidence for food, cooking, and medicinal plants entrapped in dental calculus[END_REF], who analysed organic chemical compounds found in dental calculi and directly linked them to the ingested food and to the method of cooking, was followed. Moreover, in the lack of an extensive literature about elemental composition of dental calculus, the same approach of Scott and Poulson was followed and the measured values were mostly compared with those related to teeth [START_REF] Scott | Stable carbon and nitrogen isotopes of human dental calculus: a potentially new non-destructive proxy for paleodietary analysis[END_REF].

Tb6 Just a few examples for the most relevant elements, taken from literature, are reported in Table 3.

Zinc seems to be an important element in our investigation: its values are higher than those given by Lambert and colleagues for bones belonging to carnivores [START_REF] Lambert | Ancient human diet from inorganic analysis of bone[END_REF]. Zinc is particularly abundant in foods such as red meat, dairy products, nuts, legumes (for example lentils), crustaceans (crabs and shrimps), fish, egg yolk, oats and barley [START_REF] Allmäe Raili | The content of chemical elements in archaeological human bones as a source of nutrition research[END_REF][START_REF] Beck | Bivariate analysis of trace elements in bone[END_REF][START_REF] Busetto | Dietary investigation by trace element content in bones of ancient inhabitants of Northern Italy[END_REF][START_REF] Lambert | Ancient human diet from inorganic analysis of bone[END_REF]. This element has the ability to replace calcium in bones, which is, obviously, related to its amount in the diet. Dental calculus is predominantly formed of calcium compounds and its formation depends strongly on diet and oral hygiene.

Trace metals are incorporated in bones due to a biological process, whilst in dental calculus a physical inclusion of materials in a growing matrix of calcium phosphates is the predominant process. Such a high concentration of zinc in our samples could then suggest a relevant intake of proteins, as for what concerns teeth [START_REF] Szostek | Trace elements concentrations in human teeth from a Neolithic common grave at Nakonowo (Central Poland)[END_REF][START_REF] Webb | Inductively coupled plasma-mass (ICP-MS) and atomic emission spectrometry (ICP-AES): Versatile analytical techniques to identify the archived elemental information in human teeth[END_REF].

Copper can be found in foodstuffs such as offals, red meat, fish, crustaceans, molluscs, legumes and nuts. The concentrations measured in these samples are compatible with a large consumption of animal proteins [START_REF] Dobrovolskaya | Upper Palaeolithic and Late Stone Age human diet[END_REF].

It can also be seen that strontium, a marker of a vegetarian diet, shows concentrations comparable to those tabulated by Lambert and colleagues for carnivores. A reduced consumption of cereals and other vegetables could thus be hypothesised [START_REF] Lambert | Ancient human diet from inorganic analysis of bone[END_REF]. Strontium can also derive from sea fish, shellfish and molluscs, and this evidence becomes particularly relevant when associated with a high level of zinc [START_REF] Giorgi | Oligoelements and isotopic geochemistry: a multidisciplinary approach to the reconstruction of the paleodiet[END_REF].

Manganese should be another index of the consumption of vegetables, but the very high values reported in these samples (Tb6 US 115 334.26±2.52mg/kg; Tb7 US 118 188.18±2.95 mg/kg) cannot be interpreted straightforwardly. It is reported that manganese is particularly abundant in legumes, green vegetables, cereals and nuts [START_REF] Allmäe Raili | The content of chemical elements in archaeological human bones as a source of nutrition research[END_REF][START_REF] Aschner | Nutritional aspects of manganese homeostasis[END_REF]. In that epoch, in that region, walnuts and hazelnuts were largely consumed in the diet, and even used to extract oil, so they could be a source of this metal. Another study on the chemical analysis of ancient dental calculus found a high level of manganese in a sample from an Etruscan-Celtic Necropolis [START_REF] Charlier | The microscopic (optical and SEM) examination of dental calculus deposits (DCD). Potential interest in forensic anthropology of a bio-archaeological method[END_REF]. In that case, it was attributed to the use of painted ceramic vessels, since the presence of manganese in that kind of pigments is very common. Manganese could also derive from clays used to produce kitchenware or from drinking water.

Lead is used to evaluate post-mortem contamination of the samples: because of its very low concentration values in the analysed calculus, it can be assumed that contamination did not occur. During the Middle-Ages, the use of domestic wares containing this metal, such as glazed vessels, or of drinking water coming from pipelines with lead was common (Allmäe Raili, 2012). However, the low lead levels found indicated that glazed vessels and metal piping was not used for drinking water.

CONCLUSIONS

This study has shown that the interdisciplinary study of dental calculus from ancient skulls may give important information on diet.

The dissolution of dental calculus and observation of the solid residue showed that phytoliths were mainly derived from Poaceae, Conifers and Dicotyledons (Urticals). In addition to these, it was possible to observe several fragments of burned phytoliths, micro-charcoals, a diatom fragment, fibres, starch granules, and many fragments of quartz, flint and mica.

In this study, few phytoliths, in particular of Poaceae, were observed, which is often the case [START_REF] Buchet | Comparison between the distribution of dental microstriations and plant material included in the calculus of human teeth[END_REF][START_REF] Henry | Microfossils in calculus demonstrate consumption of plants and cooked foods in Neanderthal diets (Shanidar III, Iraq; Spy I and II, Belgium)[END_REF][START_REF] Henry | Using plant microfossils from dental calculus to recover human diet: a case study from Tell al-Raqa'i, Syria[END_REF][START_REF] Lalueza | Phytolith analysis on dental calculus, enamel surface, and burial soil: information about diet and paleoenvironment[END_REF][START_REF] Lalueza | Dietary information through the examination of plant phytoliths on the enamel surface of human dentition[END_REF].

This may occur both from the low probability that, after ingestion of fruits and vegetables, such particles remained long enough in the oral cavity to stay trapped inside the calculus during its relatively fast formation, and from the kind and the amount of plant consumption.

Another study [START_REF] Licata | Gli scheletri di Caravate (VA): analisi antropologica e paleonutrizionale di una popolazione alto-medievale[END_REF] carried out by some of the authors on the same skulls studied here, examined the wear on the dentition and showed that these teeth have streaks and abrasion typical of people whose diet is based mainly on fish consumption.

Despite the small number of samples that were available for analysis, the ICP-MS showed that trace elements contained in dental calculus exhibited concentrations similar to that in samples taken from skeletons from other sites. Furthermore, the concentrations of the various elements matched well with those found in bones and teeth, with the exception of zinc. Its very high concentration could be due to a high consumption of proteins along with plant consumption that is difficult to measure, but that could be partly identified, within the limits of phytolith analysis.

This hypothesis was supported by the location of the town of Caravate, just next to Lake Maggiore, and by the high content of zinc, found by analysis with ICP-MS, which is an index of a protein-based diet.

Fragments of wood and fibres were relatively abundant and larger, and are more typical than fragmentary phytoliths frequently found in many charred food remains within the wall of the vessel of the prehistoric era and later.

In conclusion, it can be hypothesised that the diet of the Medieval inhabitants of Caravate was based on the consumption of both carbohydrates (from plants belonging to the families of Dicotyledons, Monocotyledons such as the Poaceae, also known as Gramineae, and Coniferae) and proteins, mostly coming from fish.

Figure 1 -

 1 Figure 1 -(A) -The Church of St. Augustine in Caravate, Varese, Italy; (B) -The tomb of individual Tb7 US 118

Figure 2 -

 2 Figure 2 -EDAX mapping of the sample 3.5 L Tb6 US 115

Figure 3 -

 3 Figure 3 -(A) -Trapezoidal winding phytolith (left) of the family of Poaceae, sample 4.6L Tb7 US 118; (B) -Fragmented long-shaped phytolith, width 17μm, sample 4.6L Tb7 US 118; (C) -Trichome or thorn, elongated and hollow, measuring 15μm in diameter just above the base and probably forming part of a remain of cereals, sample 4.6L Tb9 US 129

Figure 5 -

 5 Figure 5 -(A) -Fragment with conchoidal fracture that appears to derive from non-crystalline silicon dioxide (flintstone), sample 3.3V Tb9 US 129 -(B) -Fragmented rectangular perforated plate, 1.3µm with a double row of 4 holes (diameter 0.3µm) arranged orthogonally, sample 3.4L Tb6 US 115 -(C) -Objects of ellipsoidal-round shape, measuring between 40 and 60μm, sample 3.5L Tb6 US 115 -(D) -Not well definable rounded body with traces of small spines on the surface and diameter of 45μm, sample 4.6L Tb9 US 129

Table 1 .

 1 Also included in Table1is a summary of the analytical techniques used depending on the amount of available calculus.

		dental element and surface of
	individual		analysis
		provenance
		3.4 (lingual)	optical microscopy
	Tb6 US 115
		3.5 (lingual)	optical microscopy, SEM-EDX, ICP-MS
	Tb7 US 118	4.6 (lingual)	optical microscopy, ICP-MS
		3.3 (buccal)	optical microscopy, SEM-EDX
		4.4 (buccal)	optical microscopy, SEM-EDX
	Tb9 US 129
		4.6 (lingual)	optical microscopy, SEM-EDX
		4.7 (lingual)	optical microscopy, SEM-EDX

Table 1 -

 1 List of the seven samples collected and analysed, with the techniques applied

Table 2 -

 2 Measured concentrations of elements in the two dental calculus samples analysed

	US 115 3.5 lingual Tb7 US 118 4.6 lingual
	mean	std dev	mean	std dev
	(mg/kg)	(mg/kg)	(mg/kg)	(mg/kg)
	Li 0.67	0.01	0.33	0.01
	Be 0.14	0.01	0.06	0.01
	Na 3877.11	48.96	3700.95	55.89
	Mg 1851.86	14.21	1428.73	25.87
	Al 663.73	4.50	220.06	3.19
	K 364.04	19.09	539.56	25.02
	V 5.62	0.10	4.39	0.11
	Cr 3.54	0.04	7.76	0.29
	Mn 334.26	2.52	188.18	2.95
	Co 0.88	0.01	0.55	0.02
	Ni 12.29	0.21	12.08	0.63
	Cu 34.60	0.36	19.28	0.50
	Zn 867.86	8.60	456.15	2.83
	Ga 4.06	0.04	3.44	0.14
	As 0.59	0.09	0.36	0.08
	Rb 1.53	0.03	0.29	0.01
	Sr 176.30	1.64	146.04	1.49
	Cd 0.61	0.02	0.33	0.02
	Ba 79.15	0.81	68.55	1.04
	Pb 1.69	0.01	1.23	0.05
	Bi 0.25	0.01	0.11	0.01
	U 2.24	0.04	1.52	0.06

Table 3 -

 3 The two mechanisms are deeply different, but the final results both depend on food or other ways of exposition and can be considered as similar. Trace elements in teeth and dental calculus: concentration of the most relevant elements found in literature

					Pb
	Reference	Description of samples	Zn (mg/kg)	Sr (mg/kg)
					(mg/kg)
		Kalama, Egypt, contemporary population,
			472±382	143±25.3	25.5±30.3
		whole deciduous teeth
		Solís, Mexico, contemporary population,
			112±29.9	39.4±26.3	3.7±3.0
	(Webb et al., 2005)	whole deciduous teeth
		Tell Abraq, United Arab Emirates, Bronze
	Average concentration ± TCI		168±23	2197±309	0.32±0.28
		Age population, whole permanent teeth
	(confidence interval, 95%)
		New York African Burial Ground, USA,
		18th century population, whole permanent	150±51	242±88	32±42
		teeth
	(Szostek and Głąb, 2001)
		Nakonowo, Poland, Neolithic population,
			147.25±22.73 23.94±4.38	0.24±0.17
	Average concentration ±	whole permanent teeth
	standard deviation
		Lapa da Furada, Sesimbra, Portugal,	273	430	2
	(Carvalho et al., 2000)	Chalcolithic population (dentine
			180	450	2
	Median concentrations	root, dentine
		crown, enamel)	90	330	2
	(Abraham et al., 2002)
			55±58
		contemporary population, dental calculi
	Average concentration ±
		(supragingival, subgingival)
	standard deviation		280±230
		Middle Euphrates valley (Syria), Late
				90.78±27.44
		Roman Age, whole permanent teeth
	(Tomczyk et al., 2013)	Middle Euphrates valley (Syria), Islamic
	average concentration±	Age, whole permanent teeth		90.15±32.95
	standard deviation
		Middle Euphrates valley (Syria), Modern
				54.29±20.71
		Islamic Age, whole permanent teeth
	(Molokhia and Nixon, 1984)	Manchester area, contemporary patients of
	average concentration±	the University Dental Hospital, dental	221.16±63.39 37.09±3.91
	standard deviation	calculi

The examination of the samples with optical microscopy showed several dark charcoal-like fragments, which were identified as burned [START_REF] Parr | Effect of fire on phytolith coloration[END_REF] or carbon coated phytoliths. Other fragments, like micro-charcoals derived from sparks raised during the combustion in the hearth, were included in the cooking vessel and consequently found among the samples. Some fibres were also found in the samples: some ribbon-like plates, other cylindrical elements with or without chippings as well as ribbon-like (helicoidal) twisted fibres. Some other particles were found in the dental calculus samples, as can be seen in Figure 5.

A large number of inorganic components, including micas, quartzes and flint, which are likely to derive from the soil, were also found. These components could be trapped in the calculus because of the cleaning of vegetables or from containers. These components could also be the result of a residue from the milling of cereals with millstones made of sandstone, conglomerates, breccias and crystalline rocks such as granite (e.g. Figure 5a). The nature of these inorganic components was consistent with the geology of the area of Caravate.

A fragmented rectangular perforated plate (1.7μm x 1.3μm) with a double row of 4 holes (diameter 0.3μm) arranged orthogonally was found in sample 3.4L Tb6 US 115 (Figure 5b). A similar object can be found both in Dicotyledons (for example Ericaceae) and in Conifers. However, given the small size of the frustule, similar structures could also be present in the shrine of diatoms, either Centrales (marine) or Pinnales (freshwater).

The observation by optical microscope in transmitted light of the sample 3.5L Tb6 US 115 highlighted objects of ellipsoidal-round shape, measuring between 40 and 60μm (Figure 5c). Some other smaller round-oval bodies (about 10μm) with a sort of hilus characteristic of starch granules were also found but not further investigated.

This study was conducted on pictures previously taken, so a deeper investigation was not possible.

A not well-defined rounded body was found in the sample 4.6L Tb9 US 129, with traces of small spines on the surface and a diameter of 45μm (Figure 5d). It is not a pollen grain but it may be a cystolith (i.e. a crystal of calcium carbonate produced by particular plant cells called lithocites), probably an urticals tribe cystolith. Urtica sp. is an edible plant and was used as 'spinach' and for soups [START_REF] Maurizio | Histoire de l'alimentation vegetale depuis la prehistoire jusqu'a nos jours[END_REF].