Damien F G Minenna
email: damien.minenna@univ-amu.fr

Yves Elskens
email: yves.elskens@univ-amu.fr

Frédéric André
email: frederic.andre@thalesgroup.fr

Electron-wave momentum exchange and time domain simulations applied to traveling wave tubes

Keywords: wave-particle interaction, momentum, power, momentum, traveling wave tubes, TWT, Gel'fand β-transform, time domain simulations, hamiltonian dynamics

We presented at the 2013 and 2015 IVEC conferences a new discrete model for the traveling wave tube (TWT) as an alternative to PIC and frequency algorithms. Since then, we upgraded the model to consider hamiltonian aspects of the dynamics (better ensuring physical conservation laws). We also provide a new formulation for the electromagnetic momentum and power in tubes, offering further opportunities to characterize wave-electrons interactions.

I. INTRODUCTION

The wave output power (along with its ratio to the electric power used by the device) is a major industrial criterion for traveling wave tube (TWT) design, especially considering satellite payload for space systems. A recent hamiltonian approach [START_REF] André | Hamiltonian description of self-consistent wave-particle dynamics in a periodic structure[END_REF], [2] provides a new insight to the electromagnetic power.

Besides, rigorous simulations of the electrons distribution are needed to design depressed collectors in amplification devices (such as TWTs, klystrons, gyrotrons, free electron lasers).

Optimizing TWTs requires high accuracy codes reaching regimes where the electron dynamics suffers significant perturbations (trapping, chaos). In order to model those perturbations, and to observe wave oscillations and electron distribution functions, we settle for a time domain model. Numerical simulation is difficult to implement because the number of parameters exceeds the available computing power. Our project aims at building a fast and accurate algorithm to complement particle-in-cell (PIC) schemes and specialized frequency models.

II. CURRENT MODELS

Time simulations often rely on PIC codes (like e.g. CSTparticle, NEPTUNE and MAFIA), which are too slow (hours on supercomputers) for design activities, forcing industrials to rely on specialized frequency models (like MVTRAD, CHRIS-TINE, . . . ). These models are industrial standard but, based on our experience, sometimes fail to accurately predict the TWT performances (predicting oscillation phenomena, especially for electron velocity modulations, is difficult with frequency models). One reason for these discrepancies may be their not being fully consistent with Maxwell equations [3]- [START_REF] André | First comparison of new TWT discrete model with existing models[END_REF].

III. A DISCRETE HAMILTONIAN MODEL, POWER AND

MOMENTUM

To increase rapidity (to the order of seconds for a small PC configuration) and accuracy of numerical simulations (contrasting with PIC and frequency codes), we choose an approach deriving from the hamiltonian formalism and from a field decomposition [START_REF] André | Hamiltonian description of self-consistent wave-particle dynamics in a periodic structure[END_REF], [START_REF] André | First comparison of new TWT discrete model with existing models[END_REF], [START_REF] André | Hamiltonian description of the electron and its radiated field in a periodic structure[END_REF]. This approach ensures physical conservation laws [START_REF] Hairer | Geometric Numerical Integration[END_REF], enables a drastic reduction of the number of evolving variables and leads to a compact time incrementation algorithm (for example, our simulations use hundred degrees of freedom to represent the electromagnetic waves, where PIC codes employ millions).

The radiofrequency (RF) field momentum density along the TWT is expressed from the flux of the Poynting vector, and provides a new way to formulate the electromagnetic power [2], [START_REF] Minenna | [END_REF], differing from the traditional envelope approximation used in frequency models, as well as way to formulate the conserved electron-wave momentum of the system.

IV. TWT SIMULATIONS

In our current algorithm DIMO, cold fields are implemented using one-dimension eigenfunctions for the electric field depending only on experimentally known data (dispersion diagram, tube geometry, impedances). We solved an eigenfunctions normalization issue [START_REF] André | First comparison of new TWT discrete model with existing models[END_REF]. The electron evolution takes space charge into account. With our algorithm DIMO, we can also plot the distribution function of electrons in the z direction of the tube. Figures 1 and2 are plotted from the same initial conditions.

Figure 1 displays the RF electric field and the electron velocity relative change v(z, t)/v 0 -1 (arbitrary scales) in one dimension z at a typical time t during operation for an industrial TWT. Note the significant deceleration of the electrons observed near the tube collector to the right. Figure 2 shows the electromagnetic power along this TWT, computed both with our momentum density (Poynting flux) and from a Fourier analysis (black) as done in frequency models.

Oscillations, reflections and multitone operation are currently investigated. We are also working on two-and threedimensional versions of the model to provide a complete electron velocity and position distribution.

Fig. 1 .

 1 Fig. 1. Electric field (red) and electron velocity relative variation (black, with arbitrary scales), along the TWT z axis.

Fig. 2 .

 2 Fig.2. RF power computed from the envelope approximation (black) and from the Poynting flux on the hamiltonian-based model (blue and red). The blue curve is plotted at t and the red curve is plotted at t + T period /3 (viz. with a third-period phase shift) along the TWT axis. The depletion at the center is caused by an attenuator.


Josef Wilczek
email: josef.wilczek@uhk.cz

Fabrice Monna

Nicolas Navarro

Carmela Chateau-Smith

A computer tool to identify best matches for pottery fragments

Keywords: archaeology, best match, automation, pottery, ICP (iterative closest point), DCT (discrete cosine transform), RDP (Ramer-Douglas-Peucker), RTC (radius, tangent, and curvature)

Archaeologists working with pottery spend a considerable amount of time on a fundamental task -providing precise descriptions of pottery fragments. This study presents a survey of existing computational solutions to identify the best matches for a given fragment, based on its shape. Four methods (ICP, DCT, RDP, and RTC) are compared, using a pottery dataset from Graufesenque (southern France), dated to the Roman Period. The first three methods produced successful and very similar results for rim fragments (within the five best candidates for 95% of the dataset). The ICP algorithm produced the best overall results for rim fragments, and can also be used for non-rim fragments. A practical computer application, including all the above methods, was developed in R programming language, with an easy-to-use graphical interface, and is now made freely available to the archaeological community for future studies, and further development.

Introduction

Pottery is one of the most abundant materials present in archaeological excavations. It provides information about chronology, and the evolution of technique and style. It may also provide evidence of social organisation. Unlike precious artefacts belonging only to the elite, ceramics have been used by all social strata. The socio-economic dynamics of ancient populations can be reconstructed, based on pottery features: clay, fabrication technique, shape, decoration, spatial distribution, discovery context, etc. (e.g. [START_REF] Buko | Ceramika wczesnopolska. Wprowadzenie do badań[END_REF][START_REF] Orton | Mathematics in Archaeology[END_REF][START_REF] Orton | Pottery in Archaeology[END_REF][START_REF] Rice | Pottery analysis. A Sourcebook[END_REF]. However, pottery is fragile, and can be damaged or even destroyed by post-depositional processes, thus reducing the information available for archaeological inquiries.

Among archaeological investigations, one strategy is to use a set of descriptors (e.g. clay, colour, decoration, etc.) with pre-defined classes for the precise characterization of pottery fragments. Here, the main focus is on shape descriptors, which are generally related to period, origin, function, and/or aesthetics (e.g. [START_REF] Bahn | Archaeology. Theories, Methods and Practice[END_REF].

The specialist often tries to find the best match for a given fragment from within a wellestablished classification system, which can be composed of morphological classes (e.g. [START_REF] Macháček | Studie k velkomoravské keramice. Metody, analýzy a syntézy, modely. Ústav archeologie a muzeologie[END_REF][START_REF] Ness | Classification Systems with a Plot: Vessel Forms and Ceramic Typologies in the Spanish Atlantic[END_REF][START_REF] Vaginay | L'habitat gaulois de Feurs (Loire). Fouilles récentes[END_REF][START_REF] Venclová | Výroba a sídla v době laténské[END_REF]. This attribution process can be time-consuming, particularly when processing thousands of fragments. To overcome these problems, many quantitative methods have been proposed for the automatic retrieval of pottery fragments.

Almost all archaeological pottery vessels can be considered to be rotationally symmetrical. These 3D objects are thus usually represented by a 2D profile, corresponding to any crosssection between the vessel and the plane passing through its rotational axis. Quantitative methods are therefore often based on the calculation of similarities between 2D profiles, usually represented by polar or cartesian coordinates [START_REF] Liu | An XML-based information model for archaeological pottery[END_REF][START_REF] Maiza | Automatic Classification of Archaeological Potsherds[END_REF], or expressed as a function of radius (e.g. [START_REF] Jičín | K možnostem srovnávání tvarů na základě podobnosti[END_REF][START_REF] Karasik | Computerized morphological classification of ceramics[END_REF][START_REF] Mom | SECANTO-The Section Analysis Tool[END_REF][START_REF] Smith | The Pottery Informatics Query Database: A New Method for Mathematic and Quantitative Analyses of Large Regional Ceramic Datasets[END_REF], or its derivatives, such as tangent (e.g. [START_REF] Gilboa | Towards computerized typology and classification of ceramics[END_REF][START_REF] Karasik | Computerized morphological classification of ceramics[END_REF][START_REF] Leese | An Approach to the Assessment of Artifact Dimension as Descriptors of Shape[END_REF][START_REF] Main | Accessing Outline Shape Information Efficiently within a Large Database II: Database Compaction Techniques[END_REF][START_REF] Main | Accessing Outline Shape Information Efficiently within a Large Database[END_REF][START_REF] Saragusti | Quantitative analysis of shape attributes based on contours and section profiles in artifact analysis[END_REF][START_REF] Smith | The Pottery Informatics Query Database: A New Method for Mathematic and Quantitative Analyses of Large Regional Ceramic Datasets[END_REF], or curvature [START_REF] Gilboa | Towards computerized typology and classification of ceramics[END_REF][START_REF] Hristov | A Software System for Classification of Archaeological Artefacts Represented by 2D Plans[END_REF][START_REF] Karasik | New Typological Analyses of Early Bronze Age Holemouth Jars from Tel Arad and Southern Sinai[END_REF][START_REF] Saragusti | Quantitative analysis of shape attributes based on contours and section profiles in artifact analysis[END_REF]. Profiles may sometimes be represented by b-spline coefficients associated with segments of the profile curve (e.g. [START_REF] Adler | Computer Aided Classification of Ceramics: Achievements and Problems[END_REF][START_REF] Hlaváčková-Schindler | Fitting of a Closed Planar Curve Representing a Profile of an Archaeological Fragment[END_REF][START_REF] Kampel | Rule based system for archaeological pottery classification[END_REF][START_REF] Fitzgibbon | Robust registration of 2D and 3D point sets[END_REF], 2002[START_REF] Kampel | On Estimating the Position of Fragments on Rotational Symmetric Pottery[END_REF][START_REF] Laflin | Use of a Sinclair Spectrum for Shape Analysis[END_REF][START_REF] Schurmans | Advances in Geometric Modeling and Feature Extraction on Pots, Rocks and Bones for Representation and Query via the Internet[END_REF], by polylines [START_REF] Lucena | Efficient classification of Iberian ceramics using simplified curves[END_REF], or by dominant feature points that are extracted using the medialness measurement [START_REF] Piccolli | Towards the automatic classification of pottery sherds: two complementary approaches[END_REF].

Similarities between two pottery fragments can be calculated from Euclidean distances between their profiles [START_REF] Maiza | Automatic Classification of Archaeological Potsherds[END_REF], or from the similarity of the coefficients expressing their shape [START_REF] Gilboa | Towards computerized typology and classification of ceramics[END_REF]. The best match for a given fragment (and hence its morphological class) is obtained by minimising differences with potential candidates in a referential database.

The above-mentioned approaches have already been tested on real-world artefacts, with a high success rate for classification [START_REF] Karasik | Computerized morphological classification of ceramics[END_REF], and for best-match retrieval [START_REF] Lucena | Efficient classification of Iberian ceramics using simplified curves[END_REF]. It should, however, be noted that several major issues have not yet been fully addressed. For example, if the position of the outline in relation to the rotational axis is not appropriately constrained, the shape of the pottery fragment or vessel may be drastically deformed, leading to misidentification. Methods requiring successive derivatives generally suffer from numerical instability [START_REF] Karasik | 3D scanning technology as a standard archaeological tool for pottery analysis: practice and theory[END_REF].

Although the idea underlying morphological correspondence between profiles appears simple, achieving perfect correspondence between individual points on profiles is challenging, even more so when working on fragments. Many approaches have been evaluated solely on the outer surface of the profile, or only on rim fragments. Although the mathematical basis of all these approaches has already been explored (e.g. [START_REF] Hristov | A Software System for Classification of Archaeological Artefacts Represented by 2D Plans[END_REF][START_REF] Liu | An XML-based information model for archaeological pottery[END_REF][START_REF] Mom | SECANTO-The Section Analysis Tool[END_REF][START_REF] Mom | SECANTO: A Retrieval System and Classification Tool for Simple Artifacts[END_REF][START_REF] Smith | The Pottery Informatics Query Database: A New Method for Mathematic and Quantitative Analyses of Large Regional Ceramic Datasets[END_REF], no practical solution is currently available to the broader archaeological community for routine use. It is important for such proposals to be tested on a wider variety of artefacts, as no method is likely to be universally applicable.

This paper implements supervised classification, i.e. the attribution of a given fragment to one of the predefined classes. The aim of this study is to compare three existing approaches for the mathematical matching of pottery fragments based on morphology, together with a new method based on Discrete Cosine Transform (DCT). The goal is not to reconstruct complete vessels, but to identify the best matches to the fragment within the referential dataset, by shape similarity, thus indicating which shape label or labels would best suit the fragment. Attribution accuracy is evaluated on an already labelled real-world dataset. All the approaches tested are made available as a set of functions encoded in R programming language (R Core Team, 2019). An easy-to-use graphical interface was also developed, using Shiny GUI [START_REF] Chang | shiny: Web Application Framework for R[END_REF], and is made freely available to the archaeological community, to simplify data retrieval.

Material and methods

Corpus

The collection of ceramic vessels found at Graufesenque (southern France), dating from the Roman period (first to mid-second century AD), is published as an illustrated paper catalogue [START_REF] Genin | La Graufesenque (Millau, Aveyron) -Volume II -Sigillées lisses et autres productions[END_REF]. The fact that vessels in this catalogue are already labelled by a widely adopted morphological classification scheme (e.g. Brulet et al., 2012;Passelac and Vernhet, 1993;Py et al., 1993;Schucany et al., 1999;Tyers, 1996), makes it suitable for method comparison. The test corpus identified within this vast collection contains complete vessels only, and includes all morphological classes (and sub-groups) represented by ten or more vessels in the catalogue. The test corpus is thus composed of 319 vessels, including plates, bowls, cups, and goblets (Fig. 1; Supplementary Materials SM1), already divided into 14 subgroups, belonging to 10 morphological classes (Table 1).

Profile acquisition and data preparation

The drawings of the 319 vessels thus selected were scanned at 600 dpi from the paper catalogue. Profile outlines were extracted using a modified set of functions described in Claude (2008), and positioned according to their rotational axis. The outlines were automatically divided into outer and inner segments, using the rim as a reference point (Fig. 2:A).

To evaluate the potential of the four methods under comparison, two test subsets of synthetic data (rim fragments and non-rim fragments) were then produced from these outlines, by virtually fragmenting the profiles of the original complete vessels. Rim fragment analysis was based only on shape, and not on size.

For rim fragments, all outlines were size-normalised using baseline registration [START_REF] Bookstein | Morphometric tools for landmark data: Geometry and Biology[END_REF], sending the rim point and the point of intersection between rim plane and the rotational axis to the (- 1,0) and(0,0 Non-rim fragments, representing 40-45% of the profile length of the original vessel, were extracted approximately from the middle of the profile (Fig. 2:C). Each non-rim fragment was expressed by 250 points (125 each for outer and inner segments).

Best match searching algorithms

Iterative Closest Point (ICP)

The Iterative Closest Point algorithm (ICP; [START_REF] Liu | An XML-based information model for archaeological pottery[END_REF][START_REF] Maiza | Automatic Classification of Archaeological Potsherds[END_REF] is widely applied in graphics and computer vision for 3D model alignment (e.g. [START_REF] Besl | A method for registration of 3-D shapes[END_REF][START_REF] Fitzgibbon | Robust registration of 2D and 3D point sets[END_REF][START_REF] Turk | Zipped Polygon Meshes from Range Images[END_REF]. To identify the best match, a complete profile (target), is positioned in the coordinate system, and the fragment to be aligned (source), is iteratively translated, scaled, and rotated (Fig. 3:A), to minimise (using a Simulated Annealing algorithm; [START_REF] Bélisle | Convergence Theorems for a Class of Simulated Annealing Algorithms on Rd[END_REF][START_REF] Cortez | Modern Optimization with R[END_REF][START_REF] Hajek | Cooling schedules for optimal annealing[END_REF] the sum of root-mean-square deviations (RMSDs) of the distances between the points of the source and those of the target (Fig. 3:B):

min !,#,$,% ⎝ ⎛ & 1 𝑀 ) 𝑑+𝑃 & , 𝑇 ' / ( ) &*+ + & 1 𝑁 ) 𝑑(𝑃 , , 𝑇 -) ( ) ,*+ ⎠ ⎞
For the source profile, 𝑎 and 𝑏 correspond to translation along the r and z axes, 𝜑 to rotation, and 𝑆 to scaling; 𝑀 and 𝑁 are the total number of points on the inner and outer segments of the source profile. The expression 𝑑+𝑃 & , 𝑇 ' / is the distance between the 𝑗-th point, 𝑃, on the source segment and its closest point on the inner target segment, 𝑇 ' . The expression 𝑑(𝑃 , , 𝑇 -) is the distance between the 𝑘-th point, 𝑃, on the source segment and its closest point on the outer target segment, 𝑇 -.

Note that scaling (𝑆), and translation along the z axis (𝑏) are the only rigid transformations that do not alter the original shape of the vessel, unlike 𝑎 and 𝜑, which are strongly constrained by the position of the rotational axis. Ancient pottery vessels are often considerably fragmented, and were not always perfectly regular, thus making it difficult to determine the precise rotational axis to be used for profile extraction. In real-life situations, the researcher may therefore decide to relax the constraints of the rotational axis position a little, by allowing 𝑎 and 𝜑 to vary to some extent, although some shape modification will occur.

Discrete Cosine Transform (DCT)

Outlines can also be treated by Discrete Cosine Transform (DCT), used in signal treatment for compression and de-noising. This Fourier-based method transforms open outline coordinates (Fig. 4:A) into a set of harmonics (Fig. 4:B), each represented by a pair of coefficients that may be used as shape variables (e.g. [START_REF] Dommergues | The Discrete Cosine Transform, a Fourier-related Method for Morphometric Analysis of Open Contours[END_REF][START_REF] Forel | Morphometry of Middle Bronze Age palstaves by Discrete Cosine Transform[END_REF]. The higher the number of coefficients, the more precise the reconstruction of the outline (Fig. 4:B). Shape information can generally be preserved with a low number of harmonics [START_REF] Hurth | Examination of the tooth morphospace of three Mimomys lineages (Arvicolinae, Rodentia) by elliptical Fourier methods[END_REF][START_REF] Schmittbuhl | Elliptical Descriptors: Some Simplified Morphometric Parameters for the Quantification of Complex Outlines[END_REF]. The similarity between two profiles is then expressed as the RMSD between their corresponding coefficients.

The calculation of DCT is strongly constrained by the position of the rotational axis. Thus the researcher does not control variation of position, size, and orientation of the source fragment, as applied in ICP. Note that these facts also concern the two remaining algorithms: radius, tangent, and curvature, and Ramer-Douglas-Peucker polyline.

Radius, tangent, and curvature (RTC)

The profile curve (Fig. 5 The similarity between two profiles can then be obtained by calculating the sum of the RMSDs of these three functions, weighted to stress the relative importance of different parts of the profiles if required (e.g. [START_REF] Chateau-Smith | Differentiation of ceramic chemical element composition and vessel morphology at a pottery production center in Roman Galilee[END_REF][START_REF] Hristov | A Software System for Classification of Archaeological Artefacts Represented by 2D Plans[END_REF][START_REF] Karasik | Computerized morphological classification of ceramics[END_REF][START_REF] Smith | The Pottery Informatics Query Database: A New Method for Mathematic and Quantitative Analyses of Large Regional Ceramic Datasets[END_REF].

Ramer-Douglas-Peucker (RDP) polyline

Pottery profile coordinates can also be approximated by a polyline with a fixed number of segments, obtained with the Ramer-Douglas-Peucker (RDP) algorithm [START_REF] Douglas | Algorithms for the reduction of the number of points required to represent a digitized line or its caricature[END_REF][START_REF] Lucena | Efficient classification of Iberian ceramics using simplified curves[END_REF]. Given a profile curve with endpoints a and b, this algorithm seeks the most distant point on the curve from the segment ab. Once this point, c, is located, the segment ab is then replaced by two new segments, ac and cb (Fig. 6:A). The entire procedure is then repeated until the desired number of segments (𝐿) is obtained on the profile outline (Fig. 6:B-D). The points generated in this way then serve as shape variables.

Similarity between two profiles is expressed as the RMSD between their corresponding points.

Evaluating fragment retrieval

All four algorithms are a priori suitable for rim fragment retrieval, but only the ICP algorithm can be used to identify the best matches for non-rim fragments. The leave-one-out procedure was used to evaluate the algorithms: one vessel was selected from the test corpus of 319 complete vessels. A fragment extracted from this vessel was matched to the remaining 318 complete vessels in order to identify the best matches. This procedure was performed for all 319 vessels. Rim fragment retrieval was evaluated with the strategy used by [START_REF] Lucena | Efficient classification of Iberian ceramics using simplified curves[END_REF] and [START_REF] Martínez-Carrillo | Morphometric Analysis Applied to the Archaeological Pottery of the Valley of Guadalquivir[END_REF], where fragment attribution is considered correct when the original class label is present among the 𝑘 best matches. The traditional values for classification performance evaluation (with 𝑘 equal to 1, 3, and 5) were used here for rim fragments. Stricter conditions were applied for non -rim fragments, where fragment attribution was considered correct only if all 𝑘 best matches corresponded to the original class label (with 𝑘 values ranging from 1 to 10).

The acceptable ranges for parameter values used to evaluate algorithms are given in Table 2.

The optimal values required for very good approximation of the original fragment profile were 20 DCT harmonics and 20 RDP polyline segments (see Supplementary Materials SM2).

The time required for ICP calculation depends on the maximum number of iterations used in the Simulated Annealing optimisation algorithm. Rim fragments were optimised with 1000 iterations. To speed up the calculation time for non-rim fragments, the initial raw position of the source fragments on target complete vessels was estimated with the maximum number of iterations set to 500. The ten best candidates were then selected, and the optimisation procedure was repeated with the maximum number of iterations increased to 5000.

Programming

The code for archaeological pottery identification was written in R language, version 4.0.3 (R Core Team, 2021), with the aid of libraries 'sp' [START_REF] Pebesma | Classes and Methods for Spatial Data[END_REF], 'MASS' [START_REF] Ripley | Package "MASS[END_REF][START_REF] Venables | Modern Applied Statistics with S[END_REF], and 'kmlShape' [START_REF] Genolini | K-Means for Longitudinal Data using Shape-Respecting Distance[END_REF]. The graphical user interface was programmed with the 'shiny' package [START_REF] Chang | shiny: Web Application Framework for R[END_REF], combined with RStudio (RStudio Team, 2019). All software and packages used here are freely available. The application, with manual and sample data, is provided as Supplementary Materials SM3, and is also accessible via the public Git repository (https://github.com/jwilczek-dotcom/RACORD).

Results and discussion

Rim fragment retrieval

Rim retrieval results can be seen in Table 3 (see also Supplementary materials SM4 for confusion matrices obtained with 𝑘 = 1). The RTC representations, although used with success to create a classification of rim fragments dated to the Iron Age from the region of Levante [START_REF] Karasik | Computerized morphological classification of ceramics[END_REF], did not correctly classify ca. 20 % of the test dataset, while the other three methods all produced better results, reaching 95.9% for ICP with 𝑘 = 5. Interestingly, these successful retrieval rates, based only on rim fragments, were very close to those obtained from analyses of whole ceramic profiles [START_REF] Lucena | Efficient classification of Iberian ceramics using simplified curves[END_REF].

Examples of visual outputs of rim fragment retrieval can be seen in Figure 7. This figure

shows that, when a direct correspondence for the fragment is present in the corpus, the correct morphological class is identified, and the superimposition of the source fragment and target complete vessel is almost perfect (Fig. 7:1-4). Retrieval was not considered successful for rim fragments (i) with no direct correspondence in the corpus (Fig. 7:5-6), or (ii) which were attributed to a sub-group of the same morphological class (Fig. 7:7-8), or (iii)

which could have been attributed to several different classes (Fig. 7:9-10). However, the first problem would easily be identified by the archaeologist, who could therefore decide not to attribute that fragment to a specific class. The remaining two problems clearly illustrate the underlying limits of morphological attribution in archaeology.

Non-rim fragment retrieval

Non-rim fragments were judged to be attributed correctly only if all 𝑘 best candidates belonged to the same sub-group. This choice was made to show the practical use of ICP for matching, and to stress that a fair proportion of non-rim fragments, but not all, can be attributed with very high accuracy. Examples of several correctly and incorrectly labelled fragments, for 𝑘 = 5, are shown in Figure 8. Very good results can be seen in Figure 8:1-5: the superimposition of source fragments with target complete vessels is almost perfect, and all candidates belong to the same sub-group as the fragment to be labelled. Fragments without a direct correspondence in the dataset, i.e. which are not well aligned on the best candidates, can easily be identified and filtered out by the archaeologist (Fig. 8:6-7), thus minimising errors in attribution.

However, in some cases, the source fragment may seem to be well labelled when it is not -it is perfectly aligned with vessels from a different sub-group than the original (Fig. 8:8). The reason is that such fragments do not possess features characteristic of the original class (here the specific shape of the rim) and/or that some parts of vessel profiles are the same for several sub-groups, and so they cannot be discriminated solely on geometric information. As in the case of rim retrieval, such objects show clearly the limits of automatic attribution.

Concluding remarks

The four algorithms presented here propose a set of potential best matches for a given fragment, ordered by shape similarity. For three of these algorithms, the results obtained were very good, both qualitatively and quantitatively. The ICP algorithm may be recommended for all types of fragments, because it gives the best labelling outputs, and can easily be applied to all parts of the profile (rims, bodies, and bases), and to both the outer and inner parts of the profile. The basic geometric transformations on which this algorithm is based are easy to grasp, allowing the user to make pertinent choices (e.g. relaxing constraints on the rotational axis). Shape attribution with ICP is rapid, as it takes less than one second to align a given fragment. This algorithm could also be adapted to 3D data to handle fragments that are not perfectly symmetrical around the rotational axis (e.g. deformed fragments, cubic shapes, fragments with plastic features, etc.).

All the algorithms presented here were implemented in a simple computer application, with the open-source R software (Fig. 9). This program provides an extension to (semi-)automatic systems dedicated to pottery fragment orientation and profile extraction, based on 3D models (e.g. [START_REF] Karasik | 3D scanning technology as a standard archaeological tool for pottery analysis: practice and theory[END_REF][START_REF] Mara | Orientation of Fragments of Rotationally Symmetrical 3D-Shapes for Archaeological Documentation[END_REF][START_REF] Wilczek | Computer-assisted Orientation and Drawing of Archaeological Pottery[END_REF].

The user has full control over the criteria used for best-match retrieval (e.g., the range of searching values in the optimisation step). Criteria should be set in relation to the presupposed quality of the rotational axis position and orientation, the size of the objects stored in the dataset, and the number of points sampled along profile outlines. By arbitrarily fixing several criteria (e.g. the percentage of unique classes within 𝑘-best candidates), problematic fragments (e.g. typical of many classes) can be identified, and the entire process can be fully automated. The proposed class labels can then be verified, both visually by inspection of superimposed profiles, and quantitatively from the RMSD values. The quality of the reference database and existing typological schemes will obviously affect the quality of the output, whatever the algorithm used. The archaeologist always retains full control over the final and most important options with regard to fragment attribution (i.e. whether attribution is possible, and if so, into which class).

This freely available tool can be maintained, developed, and tested by the archaeological community. It can be adapted to suit the requirements of various classification strategies.

The best-match retrieval procedures can also be integrated into other existing computeraided systems for the documentation, retrieval, and classification of archaeological pottery fragments (e.g. PIQD [START_REF] Smith | The Pottery Informatics Query Database: A New Method for Mathematic and Quantitative Analyses of Large Regional Ceramic Datasets[END_REF]; ArchAIDE [START_REF] Gualandi | ArchAIDE Archaeological Automatic Interpretation and Documentation of cEramics[END_REF]; GRAVITATE [START_REF] Phillips | GRAVITATE: Geometric and semantic matching for cultural heritage artefacts[END_REF]; and more recent works [START_REF] Angelo | A 3D information framework for automated archaeological pottery archival[END_REF]). The tool can also calculate the minimum distance (i.e. similarity) between two aligned profiles. Calculated over a set of homologous profiles, these distances can be used for unsupervised classification, or to calculate shape variability, in order to explore artefact variability over space and time.

  ) coordinates (Fig.2:Bi). Rim fragments were then obtained by virtually cutting the outer and inner profile parts at a distance of 50 percent of the rim radius from the rim point, expressed by 100 equally spaced points (50 for the inner and 50 for the outer segment; Fig.2:Bii). Virtual rim fragments therefore measured 11-35% of the profile length of the original complete vessel.

  :A) can be represented by three mathematical functions: radius, tangent, and curvature. The radius function of the profile outline represents the distance between the rotational axis and each point on the profile (Fig. 5:B). The tangent function represents the angle between the tangent of each point on the profile and the rotational axis (Fig. 5:C). The curvature function represents the rate of change of this angle (Fig. 5:D).

FIGURESFigure 1 .

 1 FIGURES

Figure 2 .

 2 Figure 2. Profile acquisition and data preparation. A) Scaling the outline to the original size of the vessel, positioning it according to the axis of symmetry (i.e. z axis), and r axis, and segmenting it to the outer (red) and inner (blue) segments using the rim (red arrow) as a reference point. Bi) Size normalisation of the profile by baseline registration (here the rim point and the point of intersection between rim plane and the axis of the symmetry are set to (-1,0) and (0,0) respectively; note that the coordinate system is now expressed in rim radius units). Bii) Virtual extraction of the rim fragment: the outer and inner part of the profile is virtually broken at the distance of 0.5 rim radius along the curvilinear abscissa from the rim point, and resampled by 100 equally spaced points, 50 for the outer and 50 for the inner part. C) Virtual extraction of the non-rim fragment: the outer and inner part of the profile are extracted approximately in the middle of the vessel, and resampled by 250 points, with 125 points equally spaced along the outer and inner parts.

Figure 3 .

 3 Figure 3. Alignment process between source fragment (red polygon) and complete target vessel (grey polygon). A) The fragment is translated along the r and z axes (𝑎, 𝑏), rotated (𝜑), and scaled (𝑆) to minimise the distance with the complete vessel. B) Final alignment of the source on the target.

Figure 4 .

 4 Figure 4. Discrete cosine transform. A) The original cartesian coordinates of the outline. B)The original cartesian coordinates of the outline are decomposed into a set of harmonics. A given number of harmonics can be used to reconstruct the approximation of the original contour (here 4, 6, 8, and a full set of harmonics were used for reconstruction).

Figure 5 .

 5 Figure 5. Radius, tangent, and curvature. The first point on the outline, 𝑝, and its corresponding representations are highlighted in red. A) Position of the fragment according to rotational axis and tangent angle calculation. The tangent angle 𝜃 of the point (here 𝑝) is calculated as an angle between the tangent line passing through the point, 𝑡, and the line parallel to the rim of vessel, 𝑥. B) The radius function represents the distance between each point on the profile and the rotational axis. B) Tangent function represents the angle between the tangent of each point on the profile and the line parallel to the rim (see A for details). C) Curvature function represents the rate of change of the tangent function.

Figure 6 .

 6 Figure 6. Ramer-Douglas-Peucker polyline algorithm. A) Given a profile curve (black dots)with endpoints a and b (red dots), this algorithm seeks the most distant point on the curve, from the segment ab (dashed red line). Once this point, c, is located (red point), the segment ab is replaced by two new segments, ac and bc (red lines). 2-4) The procedure is repeated until the desired number of segments is obtained (here 5 segments).

Figure 7 .

 7 Figure 7. Evaluation of rim fragment retrieval. 1-4) Example of well-labelled fragments. 5-10) Example of incorrectly labelled fragments. Source fragments correspond to black outlines. Original complete vessels of source fragments are represented by red polygons, with correct sub-groups in red. Complete target vessels are represented by grey polygons, with correct sub-groups in grey. Rims of all fragments and vessels are scaled to the size of the radius.

Figure 8 .

 8 Figure 8. Evaluation of non-rim fragment retrieval, using ICP, for 𝑘 = 5. 1-5) Example of welllabelled fragments. 6-8) Example of incorrectly labelled fragments. A) Source fragments (darker red polygon bounded by black outlines) and their position on the original complete vessels (lighter red polygons). Source fragment sub-group code is shown in red. B) Source fragments (darker red polygons bounded by black outlines) are superimposed on the five most similar target complete vessels (grey polygons). Sub-group code for the five best complete vessels is shown in grey. Scale 1/3.

Figure 9 .

 9 Figure 9. Screenshot of the application. The source rim fragment to be labelled (G07_042_11; red polygon) is shown superimposed onto the most similar target complete vessel (G07_084_06; grey polygon) found in the test corpus, containing 319 individuals. The target vessel sub-group (here D18B) can be used to label the source fragment.

  Table 4 shows the percentage of classifiable fragments, and the accuracy of attribution, as 𝑘 increases. With five best candidates (𝑘 = 5),

	approximately 25% of fragments can be judged classifiable, with very high accuracy (93.7%
	here).

Table 1 .

 1 Test corpus used in the study (from[START_REF] Genin | La Graufesenque (Millau, Aveyron) -Volume II -Sigillées lisses et autres productions[END_REF]. For morphological class and

	TABLES
	Class	C23	D15/17		D18		D24/25		D27		D33	D35	G	R5	R8
	Sub-group	C23	D15/17B D17A D18B	D18C	D24/25	D27A	D27B	D27C	D33A	D35	G	R5	R8A
	Total	10	21	11	52	11	61	10	49	11	15	22	15	17	14
	sub-group codes: C -Curle, G -Goblet, D -Dragendorff, R -Ritterling. See Figure 1 for
	corresponding images.
		Algorithm				Rim fragments				Non-rim fragments
											𝑎 ∈ (-1 𝑐𝑚, +1 𝑐𝑚)
					𝑎 ∈ (-0.1, +0.1)					𝜑 ∈ (-3 ! , +3 ! )
	Iterative Closest Point		𝜑 ∈ (-3 ! , +3 ! ) 𝑆 ∈ (-2%, +2%)				𝑆 ∈ (-2%, +2%) 𝑏 ∈ (0 𝑐𝑚, 𝑇 " 𝑐𝑚)
					𝑎 values are expressed in rim radius units		𝑎 and 𝑏 values are expressed in real units
										𝑇 " corresponds to the height of the target profile
	Discrete Cosine Transform			20 harmonics					-
	Radius, tangent, and curvature	Equal weights for radius, tangent, and curvature			-
	Ramer-Douglas-Peucker			20 segments					-

Table 2 .

 2 Parameters and ranges of acceptable values used to evaluate algorithms. See Supplementary materials SM2 for visualisation of the reconstruction quality for 20 DCT harmonics and 20 RDP polyline segments.

	𝒌

Table 3 .

 3 Percentage of well-labelled rim fragments, when the correct class is among the 𝑘 best candidates for the four algorithms tested. The best result for each 𝑘 is highlighted in bold. See Supplementary materials SM4 for visualisation of confusion matrices obtained with

	𝑘 = 1.

Table 4 .

 4 Attribution for non-rim fragments, based on the 𝑘 best candidates. Fragments were judged classifiable only if all 𝑘 best candidates belonged to the same sub-group. Accuracy indicates the percentage of fragments judged classifiable that were attributed to the correct sub-group.

Data

Processing, Visualization, and Transmission, pp.

SUPPLEMENTARY MATERIALS

SM1. Information about the test corpus. A description of rim diameters, sub-groups, groups, and bibliographic references for a set of 319 vessels from the site of Graufesenque (southern France) dated to the Roman Period.

SM2.

The reconstruction quality of pottery fragment outlines, based on an increasing number of DCT harmonics (A), and RDP segments (B). The original outline corresponds to the black points; the shapes reconstructed by a given number of harmonics and segments are expressed as a blue curve and a red curve. Evaluations were performed using 20 harmonics and 20 segments (see black box).

SM3.

The suite of functions for best-match identification implemented in the RACORD application, with instruction manual, and examples. SM4. Four confusion matrices of rim fragment labelling when the 1 st best candidate is the correct sub-group (for 𝑘 = 1). Each table represents the confusion matrix between the original sub-groups and sub-groups obtained with corresponding method. Correct correspondences are highlighted in bold. Dashed rectangles delimit sub-groups of vessels belonging to the same group. See also Table 3.


Mao Dang

Thi Nguyen

Hanh My

Yuri Diep

Thi Nhung Ferreira Da Silva

Dongquy Vu

Chi Hoang

Ha Nhan

Thuc

Bao Quoc
email: buiquocbao@tdtu.edu.vn

Patrick Bui

Perré

Dang Mao Nguyen

Thi My

Hanh Diep

Yuri Ferreira Da Silva

Thi Nhung Vu

Dongquy Hoang

Chi Nhan

Ha Thuc

Bao Bui

Patrick Perré

Three-dimensional pore characterization of poly(lactic)acid/bamboo biodegradable panels

Keywords: biocomposites, 3D pore properties, X-ray tomography, bamboo fibers/powders, image processing

Introduction

Buildings are responsible for about 36% of global energy consumption and generate about 40-50% of global greenhouse gas emissions during construction [START_REF] Jang | Hybrid LCA model for assessing the embodied environmental impacts of buildings in South Korea[END_REF][START_REF] Hong | Framework for establishing the optimal implementation strategy of a fuel-cell-based combined heat and power system: Focused on multi-family housing complex[END_REF]. In this context, a succession of thermal regulations (RT2000, RT2005, RT2012) has been introduced with increasing thermal objectives to reduce the energy consumption of buildings and promote green construction. Similarly, important evolution, the last norm, the RE2020 regulation, introduced the concept of eco-design, by introducing new environmental standards for future construction to reduce the global emissions of the building sector. Europe is no exception, with carbon neutrality targeted by it in 2050.

Moreover, France takes part in the Energy Transition for Green Growth law, therefore requiring major changes within the sector. The environmental regulations for the new construction of buildings in France (RE2020), in force from January 2022, are fully in line with this trend. Its objectives are to continue to reduce energy consumption, strengthen performance in summer comfort and continue to promote the production of decentralized renewable energy. For this reason, one of the most effective ways to reduce building energy consumption and emissions is to carry out energy performance monitoring and diagnostics for all existing buildings as well as the development of new buildings using low-emitting component materials. Another way to reduce energy consumption in buildings would be to develop biocomposite materials [START_REF] Liu | The development history and prospects of biomass-based insulation materials for buildings[END_REF][START_REF] Palumbo | Bio-based insulation materials and their hygrothermal performance in a building envelope system (ETICS)[END_REF][START_REF] Maskell | Properties of bio-based insulation materials and their potential impact on indoor air quality[END_REF] to renew and improve indoor air quality, properly orient the building, use renewable energy and select suitable building materials. In this context, the energy recovery of waste and renewable resources plays an important role. Composites derived from agricultural by-products and bio-binders not only solve the problem of environmental pollution through the reuse of agricultural wastes but also contribute to protecting the environment by reducing CO2 emissions during production and use as well as reducing the need for petroleum-based composites. In particular, bamboo fibers are considered potential materials for making effective biocomposites because of their large availability, sustainability, durability, and hygrothermal properties [START_REF] Nguyen | Building bioinsulation materials based on bamboo powder and bio-binders[END_REF][START_REF] Nguyen | Hygrothermal properties of bio-insulation building materials based on bamboo fibers and bio-glues[END_REF]. The biocomposites fabricated from bamboo fibers were reported as environmentally friendly materials [START_REF] Rao | Outdoor bamboo-fiber-reinforced composite: Influence of resin content on water resistance and mechanical properties[END_REF][START_REF] Jawaid | Improvements in the thermal behaviour of date palm/bamboo fibres reinforced epoxy hybrid composites[END_REF][START_REF] Hu | Facile one-step fabrication of all cellulose composites with unique optical performance from wood and bamboo pulp[END_REF][START_REF] Sánchez | Physical-mechanical properties of bamboo fibers-reinforced biocomposites: Influence of surface treatment of fibers[END_REF][START_REF] Bai | Wetting mechanism and interfacial bonding performance of bamboo fiber reinforced epoxy resin composites[END_REF]. Most of the studies evaluated the physicochemical properties, morphology, mechanical properties, and water absorption capacity of these materials [START_REF] Rao | Outdoor bamboo-fiber-reinforced composite: Influence of resin content on water resistance and mechanical properties[END_REF][START_REF] Sánchez | Physical-mechanical properties of bamboo fibers-reinforced biocomposites: Influence of surface treatment of fibers[END_REF][START_REF] Bai | Wetting mechanism and interfacial bonding performance of bamboo fiber reinforced epoxy resin composites[END_REF][START_REF] Hu | A rapid, green method for the preparation of cellulosic self-reinforcing composites from wood and bamboo pulp[END_REF][START_REF] Qiu | Low velocity flexural impact behaviors of bamboo fiber reinforced composite beams[END_REF][START_REF] Supian | Mechanical and physical performance of date palm/bamboo fibre reinforced epoxy hybrid composites[END_REF][START_REF] Yusoff | Tensile and flexural properties of polylactic acid-based hybrid green composites reinforced by kenaf, bamboo and coir fibers[END_REF]. Some of them were considered insulation materials with a low thermal conductivity due to their low-density and highly porous structure.

Although playing a fundamental role in thermal insulation performance, physicochemical, and hygrothermal properties, the evaluation of porosity and pore characteristics in bio-insulation materials has not received much attention. Yet, it is wellknown that porosity and connectivity between pores have a direct effect on mass transfer properties [START_REF] Louërat | Heat and moisture diffusion in spruce and wood panels computed from 3-D morphologies using the Lattice Boltzmann method[END_REF]. In particular, the tridimensional observation and quantification of pores of these materials are still very limited, and most works focus on concrete-based materials. For example, the porosity could be obtained by identifying size parameters empirically from tomography images of concrete. The segmentation was performed by thresholding, and the porosity has been evaluated on a distinction between spherical and non-spherical pores [START_REF] Plessis | A review of X-ray computed tomography of concrete and asphalt construction materials[END_REF]. Moreover, the porosity calculated from 2D and 3D images had no difference from experimental values obtained by Archimedes' principle-based method, performed by immersing the concrete in water [START_REF] Wen | Pore characteristics analysis and numerical seepage simulation of antifreeze permeable concrete[END_REF]. These results evidence the validity of pore quantification by image processing techniques when the pore size distribution is consistent with the scans resolution. Besides, pore volumes have also been separated by pore throats. In this method, the shape of the pores was assigned to values from 1 to 6, considering pore throats as ellipses. Hence, the major axis divided by the minor axis was calculated to express how close the pores are [START_REF] Yu | Study on the pores characteristics and permeability simulation of pervious concrete based on 2D/3D CT images[END_REF]. Moreover, the porosity has been determined by the average pores area of about 20 slices of 2D tomography images divided by the slice area, and no significant differences were found compared to the porosity from 3D images [START_REF] Yu | Study on the pores characteristics and permeability simulation of pervious concrete based on 2D/3D CT images[END_REF][START_REF] Zhang | Numerical study on seepage flow in pervious concrete based on 3D CT imaging[END_REF].

To our knowledge, the pore characterization has not been extensively studied to provide an overview of biocomposites-based bamboo fibers. Therefore, in this work, we studied the porous properties of biocomposites from bamboo fibers and bamboo powders with PLA/PEG matrix using the X-ray tomography technique. The images obtained from this high-end device are processed to observe and quantify the 3D pores inside the materials.

The distribution and orientation of bamboo fibers in the whole materials are also studied and reported.

Materials and methods

Materials

Poly(lactic) acid 4043D with 94% L-lactic acid content, a number-average molecular weight (Mn) of 67 kDa, the density of 1.24 g/cm 3 , melt flow index of 0.6 g/min and the melting point of 145-160°C was purchased from NatureWorks LLC Co. Ltd (Minnesota, USA). Polyethylene glycol (PEG, Mw = 4000 g.mol -1 ) was used as a plasticizer of pure PLA and purchased from Merck Co. (Germany). The bamboo fibers (Bambusa stenostachya) were collected at the Scientific Research Center for the Conservation of Natural Resources (PhuAn Bamboo Park, Vietnam). The bamboo fibers, from 1.0 to 2.0 cm in length and from 0.16 to 0.52 mm in diameter, were prepared by the mechanic method, as reported in a previous study [START_REF] Nguyen | Hygrothermal properties of bio-insulation building materials based on bamboo fibers and bio-glues[END_REF]. The bamboo powders, about 0.1-0.2 mm in diameter), were crushed from bamboo fibers as reported in the previous study [START_REF] Nguyen | Building bioinsulation materials based on bamboo powder and bio-binders[END_REF].

Bamboo fibers treatment

To remove impurities, unexpected substances and improve the mechanical properties of biocomposites, bamboo fibers were pretreated with a NaOH solution. Firstly, the bamboo fibers were soaked in a 2% NaOH solution for 24 h at room temperature [START_REF] Bui | Mechanical and thermal properties of a biocomposite based on polyvinylchloride/epoxidized natural rubber blend reinforced with rice husk microfiller[END_REF]. Then, they were taken out and washed thoroughly with distilled water until a neutral solution was obtained after washing (pH  7). Then, the fibers were dried at 60°C for 24 h before using for further experiments.

Fabrication of biocomposite materials

Biocomposites containing pure PLA, and bamboo fibers or bamboo powders were manufactured by the melting process method using a twin-screw Haake Polydrive device (Germany) with different components, as listed in Table 1. 10 wt.% PEG was used as the effective plasticizer for PLA as reported in the previous study [START_REF] Nguyen | Valorization of starch nanoparticles on microstructural and physical properties of PLAstarch nanocomposites[END_REF]. Firstly, the components were preliminarily mixed and filled into a chamber mixer at a rotation speed of 60 rpm, at 170 o C for 10 min. Then, the mixture was shaped by an aluminum mold of 12×12×0.3 cm 3 using a hydraulic press at 185 o C, 50 kg/cm 2 for 15 min. The composite plates were then cooled to room temperature and kept in the mold for 24 h before removing them to avoid warping. Finally, the plates were then prepared for further characterizations. 500 to 4000 cm -1 . The bamboo fibers were dried to remove the water vapor before being analysed.

Morphology assessment

The interfaces between matrix and bamboo fibers and the morphology of biocomposites were assessed using a scanning electron microscope (SEM, JSM 6600, JEOL Co., Japan).

The cryo-fractured specimens were first coated with platinum and observed after freezing in liquid nitrogen.

3D microstructure of biocomposites

The tridimensional microstructure of the composites was analyzed by combining X-ray tomography (EasyTom XL 150/160, RX-solutions, France) and image processing in the software Avizo 2019.1 (Thermo Fisher Scientific, USA). To obtain X-ray projections, samples were scanned using a micro-source, a CCD detector of 2004 ×1336 pixels and a voxel size of 5 μm. For these acquisition parameters, each scan lasted approximately 30 min and about 1100 projections were recorded. These images were processed in X-Act software 21.04.3 (RX-Solutions, France) and reconstructed to obtain image stacks, which were further analyzed in Avizo [START_REF] Perré | A macroscopic Washburn approach of liquid imbibition in wood derived from X-ray tomography observations[END_REF]. 3D volume visualization of the stacks was performed by applying the Volume Rendering tool (Colormap = "Grayscale"). Further, subvolumes (3 x 3 x 6 mm) were extracted from the original datasets and processed using a 3D Gaussian filter (Standard deviation = 1 px, kernel size = 2). After filtering, segmentation was performed in 3D using the Auto Thresholding tool (Auto Threshold Low, Mode = min-max, Criterion = Entropy) to identify sample pores. The obtained binary images were analyzed using the Material Statistics function to calculate the pore volume fraction in the composites. In addition to that, the Labeling tool was applied to assign pores as individual particles, allowing the quantification of their volumes with the Label Analysis function. As fibers presented a similar intensity to the composite matrix, their segmentation could not be performed by automatic thresholding. Hence, in order to visualize their distribution in the composites, the Volume rendering tool was applied using a different colormap ("Physics") and a reduced opacity (0.02).

Thermal stability

The thermal stability of biocomposites was evaluated by Thermogravimetric analysis (TGA) using a TGA Q500 (TA Instruments, USA) under a nitrogen atmosphere. The heating rate was 10°C/min with a temperature range from 25 to 600°C.

Flexural properties

The flexural properties of biocomposites were characterized using a mechanical testing machine (Comtech, Taiwan). The specimen's dimension and testing protocol were prepared and performed according to ASTM D790. The obtained flexural results were the mean values of at least 5 measurements for each composite.

Water absorption

The water absorption of biocomposites was measured according to the ASTM D1348 standard. The specimens were immersed in water and mass change was monitored over time. Every 24h, the specimens were collected and their mass was weighed after the excess of water was wiped off from the surface. The measurement was stopped when the specimen mass was constant (mass variation between two measurements less than 5%).

The amount of water absorbed into specimens was an average value of at least 5 samples and calculated by the formula:

WA (%) = 100 × (m -m 0 ) m 0
where WA (%) is water absorbed into specimen; m (g) is mass of specimen after immersing in water; and m0 (g) is initial mass of specimen.

Results and discussions

Bamboo fibers treatment

FTIR technique is used to assess the chemical composition and structural changes in bamboo fibers before and after treatment with NaOH solution, as exhibited in Figure 1. The FTIR spectrum of bamboo fibers before being treated with NaOH solution shows low-intensity overlapping peaks, which do not allow us to distinguish the typical vibrations of cellulose, hemicellulose, and lignin. However, the FTIR spectrum of the bamboo fibers treated with NaOH solution presents discrete peaks with high intensity, such as the band around 3000-3700 cm -1 , characteristic of the -OH groups present in the main constituents of bamboo fibers.

A low-intensity peak at 1732 cm -1 could not be observed in the FTIR spectrum of treated bamboo fibers. This peak could be associated with the oscillation of -C = O groups of acetyl and esters of hemicelluloses or with constituents from bamboo lignin, such as esters, ferulic and p-coumaric acids [START_REF] Li | Cellulose whiskers extracted from mulberry: A novel biomass production[END_REF][START_REF] Elanthikkal | Cellulose microfibres produced from banana plant wastes: Isolation and characterization[END_REF]. A peak appears at a wavelength of 1508 cm - A change of a small peak is identified at 896 cm -1 , corresponding to -C-H deformation vibration in cellulose [START_REF] Elanthikkal | Cellulose microfibres produced from banana plant wastes: Isolation and characterization[END_REF]. Thus, the FTIR spectrum of the bamboo fibers after treating with 2 wt.% NaOH solution shows that the impurities and a part of lignin and hemicelluloses were removed from bamboo fibers, which was intended to improve the bamboo fiber/matrix interface.

Furthermore, the FTIR spectra of biocomposites are exhibited in Figure 2. Accordingly, the oscillations observed from the FTIR spectra of biocomposites reinforced with treated bamboo fibers and powders do not differ because of the same interaction between bamboo and matrix PLA/PEG. This is also observed for the FTIR spectra of biocomposites with different bamboo fiber contents. However, the peak at 1735 cm -1 representing carbonyl (C=O) stretching vibration peak does not appear in treated bamboo fibers but it is observed in all biocomposites. This result proves the creation of the chemical bonding between functional groups from PLA/PEG matrix and bamboo fibers. On the other hand, the flexural properties of the PLA-composite reinforced with bamboo powders show a significant increase compared to with the reinforcement by bamboo fibers. This implies that bamboo powders exhibit a higher reinforcing role than bamboo fibers in PLA/PEG composite. Accordingly, the flexural modulus and flexural strength of PLA/PEG-composite reinforced of 50 wt.% bamboo powders are raised up 6219 MPa and 61.5 MPa, respectively, which are higher than that of PLA/PEG composite reinforced with bamboo fibers at 40 wt.%. This result may be surprising as the shape factor is important in composites, and we would expect the fiber composites to perform better than the powder composites. This observation is therefore probably explained by a better interface between matrix and particles in the case of powder: the smaller sizes of bamboo powders allow a more efficient mixing during the melting processing in the mixing chamber [START_REF] Nguyen | Building bioinsulation materials based on bamboo powder and bio-binders[END_REF]. Moreover, with a smaller size than bamboo fibers, the bamboo powders are able to fill in the large voids in the composites, reducing the pore volume in the composite bulk and at interfaces between PLA/PEG matrix and bamboo powders. This result is consistent with the quantitative results of pore volume by tomography technique and image processing displayed in Figure 6.

However, when the bamboo powders content in composite increases to 60 wt.%, its flexural modulus still increases up to 6444 MPa. However, its flexural strength tends to decrease to 59.1 MPa. This result shows a decrease in flexural strength when the reinforcement content in the composite is too high that the PLA/PEG matrix no longer exhibits as good adhesion as the composite with lower bamboo powder content (PP50). It means that the PLA/PEG matrix is not sufficient to wet all bamboo powders in this composite [START_REF] Nguyen | Hygrothermal properties of bio-insulation building materials based on bamboo fibers and bio-glues[END_REF][START_REF] Yusoff | Tensile and flexural properties of polylactic acid-based hybrid green composites reinforced by kenaf, bamboo and coir fibers[END_REF][START_REF] Chin | Thermal and mechanical properties of bamboo fiber reinforced composites[END_REF]. Therefore, in this study, the biocomposites reinforced with the high bamboo powders contents (> 50 wt.%) are not necessary and useful for further studies.

3D visualization of bamboo fibers orientation in composites

3D volume morphology

The 3D volume of the selected composites is reconstructed from a series of 2D images of each composite. The components of the composites including bamboo fibers, pores, and PLA/PEG matrix are distinguished by their density in the X-ray projections. Figure 4 shows 3D volume reconstructions, corresponding to three PLA composites reinforced with bamboo fibers (PF30 and PF40) and bamboo powders (PP50). The 3D volume morphology of the composite containing 30 wt.% (PF30) bamboo fibers exhibits very few pores, which are scattered throughout the sample. This is due to the low content of bamboo fibers, so the matrix has completely covered the whole fibers and filled the gaps inside this composite. When the fiber content is increased to 40 wt.% (PF40), the pores appear in larger size due to the rigidity of bamboo fibers. For the composite containing 50 wt.% bamboo powders (PP50), the pores are more abundant but evenly distributed throughout the sample. However, the pores seemed to be smaller in size compared to composites containing 40 wt.% bamboo fibers. This is explained by the small size of bamboo powders allowing the matrix to better penetrate during fabrication.

3D fibers orientation

X-ray tomography and image processing techniques also allow us to visualize the distribution and orientation of bamboo fibers within the composites, as presented in The bamboo fibers are observed with different sizes in the form of bundles of interwoven or stacked fibers. Besides, spaces between the adjacent fibers could also be observed.

When the bamboo fibers content was increased to 40 wt.% (Figure 5b), the fibers remained randomly dispersed and more fiber agglomeration could be observed due to an increase in the number of fibers in this composite. The fibers sizes could still be visualized, however, compared to the PF30, the spaces between the fibers could not be identified. Indeed, the bamboo fibers are mainly presented as folded lines and curves in the composite. Bamboo fibers distribution could be affected by mixing time and fiber content, while the fiber orientation depends on several factors, such as fabrication method and coarse aggregate content [START_REF] Qin | Microscopic Texture of Polypropylene Fiber-Reinforced Concrete with X-Ray Computed Tomography[END_REF]. In the composite containing 50 wt.% bamboo powders (Figure 5c), the bamboo powders density is high. Particles tend to agglomerate and the distance between the grain boundaries is not observed. Moreover, the particles size is also indistinguishable. This is due to the excessive bamboo powders loading and the applied tomography resolution, which does not allow the observation of bamboo powder size.

3D visualization of pore distribution

The pore size distribution and pore volume fraction are two fundamental parameters to express the pore characteristics. The pore size distribution has been studied in 3D and the pixel (or voxel) number count is the most convenient size measurement method. The 3D visualization of segmented pores in the composites is presented in For the PF30 sample, the pore quantity is low, the pores are unevenly distributed inside the composites, and small pores are clustered together in longitudinal arrays. This is explained by the porosity formed at the interface between the bamboo fibers and the matrix. This result is confirmed through the inserted 2D SEM image and depicted through the red arrows in Figure 6. The pore distribution in three composites are also investigated as reported in Figure 7. In order to avoid the huge number of small pores, the categories were quantified by their cumulated volumes. For a better appreciation of these distributions, we proposed both a linear and log scales. This figure confirms that the total porosity is very small for sample PF30, with a slightly decreasing distribution of volume per category over the total range, 0 to 10 mm 3 . Sample PF40, depicts a huge increase of small pores (5000 mm 3 instead of 170 mm 3 for the first category: 0-0.2 mm 3 ), but also large volumes of medium-size pores and a significant presence of large pores, as can be seen on the images. At this level of solid fraction, this is explained by the rigidity of fibers that might preserve cavities, either during forming in the internal mixer or by subsequent relaxation. In addition, the pore features could also be affected by the characteristics of fibers. For instance, the shape of pores showed a tendency from spherical to oval with fibers due to the distortion of the matrix [START_REF] Ríos | Analysis of the tensile fracture properties of ultra-high-strength fiber-reinforced concrete with different types of steel fibers by X-ray tomography[END_REF]. The bar graph of PP50 composite confirms the large porosity, but with regular spatial distribution of pore size. Compared to the previous distributions, the huge increase of small-and medium-size pores (up to 4 mm 3 ) is obvious. Above this value, contrary to the previous samples, a rapid drop in category volumes is observed for large pores, especially visible on the linear scale.

This partition and amount or porosity points out the role of the particle size: large and long particles are likely to produce a small number of large pores due to the stiffness of biomass. They build up weakness points at high mass fractions and explain the poor mechanical behavior, which seems contradictory with the positive effect of the shape ratio of particles. On the contrary, with powder, the porosity might be large, again as a function of mass fraction, but the pores remain of average size and remain evenly distributed in space. In conclusion, the particle size matters more than the shape ratio for high particle/PLA ratios. This explains the apparently paradoxical behavior of the mechanical properties.

The total volume pore per category of three samples was also plotted in Figure 7 to correct the bias in number of pores, generally much larger for small pores. The partition of volume per category confirms a small porosity for PF30, which is consistent with total value reported in Table 2. PF40 has an increased porosity but keeps a similar balance between small and large pores. The bamboo powders composite, PP50, presents the largest total porosity, but, contrary to the fiber sample, the proportion of small pores is larger, even when quantified in terms of volume rather than number. This depicts a much regular pore distribution, likely to explain the best mechanical performances. This is even more meaningful if, in addition to this observation, we keep in mind that the larger the pore volume, the higher the presence of interconnected pores, which are counted as one large pore by the software. This is also confirmed by the increased kinetic in water uptake result.

The pore volume fraction of the 3 biocomposites is also calculated based on the cumulated pore volume from the composites and presented in Table 2. The PF30 sample has a very low porosity, about 0.87%, the PF40 sample has a porosity of about 2.73% and the PP50 sample has a porosity of about 12.12%. When the bamboo fibers content in composites is increased, the porosity increases because of increasing in the rigidity of fibers. This result is also consistent with the concrete materials in the literature, indeed, the porosity of the concrete containing polypropylene fibers increased alongside fiber's content [START_REF] Simões | Influence of fibres on the mechanical behaviour of fibre reinforced concrete matrixes[END_REF]. Furthermore, these values are consistent with the number and volume of pore observed in Figure 7. It could be observed that the porosity of biocomposites is increased with biomass content, whatever it is fiber or powder. However, the smaller pore size and better distribution obtained with powder composite is likely to explain the best mechanical performances.

Thermal stability of biocomposites

The thermal stability of the PLA/PEG matrix and its composites is shown through the TGA curves, as illustrated in The powders appear to be agglomerated in the tomography images. Maybe the small pores in the image could be caused by the emission of volatile compounds during the composite fabrication, generating bubbles in the matrix that are trapped after cooling and solidification of the polymer. It could be explained due to mechanical relaxation, indeed, when the level of fiber increases and pore formation occurs due to the relaxation of bend fiber, that leads to large pore because bending deformation is large. However, in the PP50 sample, the relaxation occurs more to release cell compression, which generates smallerdistributed pores of similar volume.

Moreover, maybe the NaOH treatment only modifies fibers on the surface. The interior, which is revealed after the powders are crushed, may still present volatile substances, which we see in the TGA and can be emitted during extrusion, generating bubbles that could result in the porous structure observed by tomography.

Water absorption

The water absorption capacity of the biocomposites is shown in Figure 9. The PLA/PEG matrix exhibits a good water resistance, corresponding to the lowest and seemingly negligible water absorption, less than 0.6% for 11 days. The water absorption capacity of biocomposites increases with the increasing loading contents and experimental time.

Indeed, the composite containing 30 wt.% bamboo fibers absorbs about 9% and reaches the water saturation after 9 days of test. The same water saturation period for the sample which contains 40 wt.% bamboo fibers with the water absorption maximum of about 11%. The composite containing 50 wt.% bamboo powders (PP50) has the highest hygroscopicity with a water absorption amount of about 13% and reached water saturation after 5 days. These results imply that bamboo fibers and bamboo powders have increased the hygroscopicity of the biocomposites. Interestingly, the quantity of absorbed water directly depends, with values of 30%, 27.5% and 26% respectively for PF30, PF40 and PP50.

Conclusion

Biocomposites from bamboo fibers, bamboo powders and PLA/PEG matrix have been fabricated by melting method. The flexural properties of the materials are optimal at 40 wt.% bamboo fibers, but is higher, 50 wt.%, for bamboo powders. This result was explained by 3D observations, using an x-ray tomography. 3D volume reconstruction, fiber orientation and pore characteristics of composite materials such as volume fraction and size distribution have been determined by image processing using the Avizo software. Accordingly, the porosity of the materials increases with the content of the bamboo reinforcement phase. More specifically, the composite material containing 50 wt.% bamboo powders has numerous pores of average size, well distributed spatially, while that of bamboo fibers have a few numbers of large pores, likely to explain the difference of mechanical behavior. Consistently with the water-uptake capacity of bamboo, the water absorption of the composite material increases linearly with the content of the bamboo reinforcement phase. In addition, the average-size pores in the powder composite are sufficiently numerous to percolate and significantly accelerate the dynamic of water uptake.

1 .

 1 Fourier transform infrared spectroscopy (FTIR)The structure change of bamboo fibers after being treated with 2% NaOH solution was investigated by Fourier transform infrared (FT-IR) spectrometer Equinox 55 (Bruker Co., Germany) with an absorption mode at a resolution of 2 cm -1 in wavenumber ranges from

Figure 1 .

 1 Figure 1. FTIR spectra of bamboo fibers before and after treating with NaOH solution

1 ,

 1 representing the aromatic -C = C-stretching of the aromatic rings of lignin. The strong cellulose bands around 1113 cm -1 correspond to antisymmetric -C-O-Cbridge stretching.

Furthermore, the intensity

  of this 1735 cm -1 peak does not change with bamboo fiber contents in biocomposites. It implies that the more addition of bamboo fibers does not change the chemical bonding of the PLA/PEG matrix and bamboo fibers during internal melt processing.

Figure 2 .Figure 3 .

 23 Figure 2. FTIR spectra of biocomposites

Figure 4 :

 4 Figure 4: 3D volume reconstruction of PLA composites reinforced with: (a) 30 wt.%

Figure 5 .

 5 Figure 5. Accordingly, the images of the PLA/PEG composite containing 30 wt.%

Figure 5 :

 5 Figure 5: 3D bamboo fibers orientation in composites: (a) PF30, (b) PF40, and (c) PP50

Figure 6 .Figure 6 :

 66 Figure 6: 3D visualization of pore in the composites: (a) PF30, (b) PF40, and (c) PP50

Figure 7 :

 7 Figure 7: Pore distribution and total volume pore per category in the composites with

Figure 8 .Figure 8 :

 88 Figure 8: Thermal decomposition curves of biocomposites

Figure 9 :

 9 Figure 9: Water absorption ability of biocomposites

Table 2 :

 2 The pore volume fraction of composites

	Material	Porosity (%)
	PF30	0.87
	PF40	2.73
	PP50	12.12

  The thermal stability of PLA/PEG reinforced with 50 wt.% bamboo powders (PP50) is reduced compared to the PLA/PEG matrix. Accordingly, the first decomposition is observed before 150 o C, which is the decomposition of low molecular weight and volatile substances. In fact, the bamboo exists in powders form with small sizes (< 20 micrometers), containing compounds which are easily decomposed by heat. Therefore, the decomposition weight of this step is about 5%. The thermal decomposition starts at about 260 o C and ends at 350 o C, this is thought to be the thermal decomposition of PLA/PEG mixed with bamboo powders with a lower temperature than PLA/PEG only.

This result is due to the lower thermal stability of bamboo powders compared to PLA/PEG matrix. Moreover, the first parts of the curves are completely similar for net PLA/PEG and PF40, as if PLA is likely to protect fibers. The large porosity of PP50 could explain that this protective effect doesn't exist for powder. The second decomposition step starts at about 350 o C and ended at 500 o C, corresponding to the thermal decomposition of the heat-stable components in bamboo powders. The residual of this sample after thermal decomposition is about 20%.

Acknowledgments

This research is funded by Vietnam National University HoChiMinh City (VNU-HCM) under grant number QGB11.ĐA 01.

Conflicts of Interest:

The authors declare no conflict of interest.


Damien F G Minenna
email: damien.minenna@univ-amu.fr

Yves Elskens

Member, IEEE Jérôme Puech

Alexandre Poyé

Fabrice Doveil

Telma Pereira

DIMOHA : Traveling-wave tube simulations including band edge and multiple-carriers operations

Keywords: traveling-wave tubes TWTs, DIMOHA, simulation, band edge, multiple carriers, intermodulation, AM/AM, AM/PM, nonlinear signals, time domain, wave-particle interaction, Hamiltonian dynamics, N-body dynamics, communication systems, large-signal amplifiers

A new kind of time domain simulation for travelingwave tube (TWT) is presented. The DIMOHA algorithm should be considered as an alternative to current particle-in-cell (PIC) and frequency approaches. To assess its validity, a comparison with a commercial Ku-band tapered helix TWT is performed. It shows an excellent agreement for the output power and the nonlinear phase shift at various frequencies. The algorithm is also tested near the band edge and for multitone simulations, both with success.

I. INTRODUCTION

Traveling-wave tubes (TWTs) (see our recent historical review [START_REF] Minenna | The Traveling-Wave Tube in the History of Telecommunication[END_REF]) are high-power amplifiers mostly used for space communications. The physical process at the heart of the device is the wave-particle interaction, where by electrons transfer their momentum to amplify telecommunication signals, as, e.g., in Landau damping. A synchronization process [2] occurs when the phase velocity of the wave (determined by the tube geometry) is close to the particle speed.

Nowadays, it is still difficult to simulate the output wave since TWTs modulate signals in amplitude (this is the role of an amplifier), in phase (knowing that the information is coded on this phase), and in frequency (harmonic and intermodulation generations). In addition, the large power of high frequency TWTs and their broad band spectrum generate critical instabilities. For example, spurious frequencies are generated very far from the drive frequency in the saturated nonlinear regime (where electron motion perturbations become significant) because of drive-induced oscillations.

Moreover, to improve the design of depressed collectors, we need accurate simulations of the electron distributions, allowing to lower the electric power used by the device, especially considering satellite payload for space systems. II. THE ALGORITHM Currently, there are two kinds of numerical approach to simulate the wave propagation in TWTs, each with typical benefits. The first kind relies on kinetic models in time domain such as particle-in-cell (PIC) codes based on kinetic equations (like, e.g., CST-PARTICLE, NEPTUNE or MAFIA). The second kind is based on specialized frequency models (like, e.g., MVTRAD, CHRISTINE or TESLA). In contrast, we use a specialized time domain model relying on an Nbody (a.k.a. many body or molecular dynamics) description [3], [START_REF] Minenna | Proc. IEEE Int. Vacuum Electron. Conf. (IVEC)[END_REF]. By definition, this description is deemed slower and computationally more expensive than PIC since all particles should be considered. We overcome this issue thanks to a drastic degrees-of-freedom reduction [START_REF] André | [END_REF].

Our algorithm DIMOHA (for HAmiltonian DIscrete MOdel) is based on a rigorous theoretical model that respects Maxwell equations and Hamiltonian conservation properties, even with the number of degrees-of-freedom reduction. The heart of this model is the momentum exchange between the wave and particles, within the Hamiltonian formalism using a field decomposition. A particularity of DIMOHA is that it allows one-dimensional simulations, where generally time domain codes require to simulate the full tube in three dimensions. But its main advantage is that it considers all the fields in time (not just the wave envelope) as well as the electron dynamics (positions and velocities). This enables DIMOHA to work in nonlinear regime (trapping, chaos) and to inject any wanted type of wave, like multi-tones and telecommunication signals. A comparison with MVTRAD has already shown [3] a perfect agreement at power saturation in the middle of the transmission band of a basic tube without taper nor attenuators.

III. COMPARISON WITH MEASURES

The current version, that supports parallel computing (multicore and multi-processor) on desktop computers or clusters, is able to perform a full simulation of industrial tubes. Figure 1 displays the RF output power (AM/AM) and the nonlinear phase shift (AM/PM) at F = 11.7 GHz for a commercial Ku-band tapered helix TWT, comparing DIMOHA with direct measurements. The computing time is less than 10 minutes to reach steady state on a small desktop computer (less than 800 US$). It is a good progress compared to the 12 h needed with CST-PARTICLE [6].

A recent comparison [START_REF] Minenna | Degrees-of-freedom reduction for wave-particle systems : Comparison of Pierce equivalent circuit and recent discrete model for traveling-wave tubes[END_REF] with the Pierce four-waves theory has shown that our theoretical model should be more accurate near band-edge [START_REF] Ryskin | [END_REF]. Indeed, Figure 2 As shown, the cathode potential V h of MVTRAD is shifted by 0 to -100 V to adjust either the power output or the phase shift but not both at the same time. The phase shift of DIMOHA is identical to measurements (Fig 1).

for the same tube. A comparison with MVTRAD (that suffers the same shortcoming as Pierce's theory) is performed. The small variations near both edges between DIMOHA and the measurement result from our incomplete management of tube features, and we are working to reduce this discrepancy. Electromagnetic waves for telecommunication signals are not simple sinusoidal waves and DIMOHA is able to simulate them without limitations. We propose here a single example where two-tones are injected in the same Ku-band TWT. of the bench and this TWT is not adapted in this band. The computation time is about 5 h for 100 ns.

IV. PERSPECTIVES We are currently working to improve the multitone simulations (including more complex nonlinear signal analysis) and to improve the modeling of tube defects (taper, attenuators, hot voltage standing-wave ratio (VSWR), etc.) for the characterisation of oscillations and reflections. In addition, a possible three-dimensional version of the algorithm is investigated. We also work on the electron distribution function with the aim to improve depressed collectors. Finally, we consider using the algorithm to study nonlinear and chaotic regimes in various aspects, including simulations of the 4 meters long TWT at Aix-Marseille Université, CNRS.

Fig. 1 .Fig. 2 .

 12 Fig. 1. Output power and nonlinear phase shift for a Ku-band tapered helix TWT and for DIMOHA at F = 11.7 GHz.

Fig. 3 .

 3 Fig. 3. Comparison of the output spectrum in two-carriers regime of DIMOHA (after t = 100 ns) and a Ku-band tapered helix TWT. Input: F 1 = 11.65 GHz and F 2 = 11.75 GHz at P in,1,2 = -10 dBm.

ACKNOWLEDGMENT

The authors thank Lénaïc Couëdel, Guillaume Fuhr and Geoffroy Soubercaze-Pun for fruitful discussions.