Jérémy Frey
email: jeremy.frey@inria.fr

Renaud Gervais
email: renaud.gervais@inria.fr

Stéphanie Fleck
email: stephanie.fleck@univ-lorraine.fr

Fabien Lotte
email: fabien.lotte@inria.fr

Martin Hachet
email: martin.hachet@inria.fr

Teegi: Tangible EEG Interface

Keywords: Tangible Interaction, EEG, Spatial Augmented Reality, Learning ACM Classification H.5.1 Multimedia Information Systems: Artificial, augmented, and virtual realities, H.5.2 User Interfaces: Interaction styles, H.1.2 User/Machine Systems: Human information processing, I.2.6 Learning: Knowledge acquisition

We introduce Teegi, a Tangible ElectroEncephaloGraphy (EEG) Interface that enables novice users to get to know more about something as complex as brain signals, in an easy, engaging and informative way. To this end, we have designed a new system based on a unique combination of spatial augmented reality, tangible interaction and real-time neurotechnologies. With Teegi, a user can visualize and analyze his or her own brain activity in real-time, on a tangible character that can be easily manipulated, and with which it is possible to interact. An explorative study has shown that interacting with Teegi seems to be easy, motivating, reliable and informative. Overall, this suggests that Teegi is a promising and relevant training and mediation tool for the general public.

INTRODUCTION

Electroencephalography (EEG) measures the brain activity of participants under the form of electrical currents, through use of a set of electrodes connected to amplifiers and placed on the scalp [START_REF] Niedermeyer | Electroencephalography: basic principles, clinical applications, and related fields[END_REF]. This technology is widely used in medicine for diagnostic purposes. It is also increasingly explored in the field of Brain-Computer Interfaces (BCI), the goal of which is to enable a user to send input commands to interactive systems without any physical motor activities, by using brain activity alone [START_REF] Wolpaw | Principles and Practice[END_REF]. BCI is an emerging research area in Human-Computer Interaction (HCI) that offers new opportunities for interaction, beyond standard input devices [START_REF] Tan | Brain-Computer Interaction: Applying our Minds to Human-Computer Interaction[END_REF]. These emerging technologies are becoming increasingly more popular. It feeds into fears and dreams in the general public where many fantasies are linked to a misunderstanding of the strengths and weaknesses of such new technologies. No, it is not possible to read thoughts! But what can be done exactly? Our motivation is to provide a tool that allows one to better learn how EEG works, and to better understand the kinds of brain activity that can be detected in EEG signals. Beyond the knowledge of the brain that a user can acquire, we believe that a dedicated tool may help demystify BCI, and consequently, it may favor the development of such a promising field.

We followed a multidisciplinary approach, combining Human-Computer Interaction (Spatial Augmented Reality, Tangible User Interfaces), Neurotechnologies (EEG, brain signal processing) and Psychology/Human sciences (Human Learning and Representations, Scientific Mediation) to design an interactive multimedia system that enables novice users to get to know more about something as complex as EEG signals and the brain, in an easy, engaging and informative way. Our final goal is to enhance learning efficiency and knowledge acquisition by letting users actively and individually manipulate and investigate the concept to be learned [START_REF] Vosniadou | Designing learning environments to promote conceptual change in science[END_REF], i.e. EEG signals.

This gave birth to Teegi (Tangible EEG Interface), a physical character that users can manipulate in a natural way to observe and analyze their own brain activity projected in realtime on the character's head (see Figure 1). Beyond the technical description of Teegi, this paper depicts an explorative study we conducted, which provides an experimental basis for discussions and future works. Our major contribution for this paper is the design of the first system to make EEG signals and brain activity easily accessible, interactive and understandable. This work is based on theoretical foundations, technical developments, and preliminary investigations.

NEUROIMAGING AND EEG

EEG signals are small electrical currents (in the µV range) that can be measured on the surface of the scalp [START_REF] Niedermeyer | Electroencephalography: basic principles, clinical applications, and related fields[END_REF]. They reflect the synchronous activity of millions of neurons from the brain cortex (i.e., the outer layer of the brain). Compared to alternative neuroimaging techniques, such as MagnetoEn-cephaloGraphy or functional Magnetic Resonance Imaging, EEG is simultaneously cheap, portable and provides good time resolution. Because of these advantages, EEG has been used for many years in medicine, e.g., for the diagnosis of sleep disorders or epilepsy [START_REF] Niedermeyer | Electroencephalography: basic principles, clinical applications, and related fields[END_REF]. More recently, with the advance of computer processing performance, it became possible to measure and analyze in real-time the content of these EEG signals. This paved the way for the rise of BCI which uses real-time analysis and decoding of EEG signals in order to identify the mental state of the user and translate it into a command for an application [START_REF] Wolpaw | Principles and Practice[END_REF].

The currently available tools used to visualize and analyze such signals are tailored for experts with a deep understanding of the brain, EEG principles and EEG signal processing [START_REF] Niedermeyer | Electroencephalography: basic principles, clinical applications, and related fields[END_REF]. Although EEG visualizations are intended for experts only, the general public is often compelled by how the brain works and how its activity is measured. Anyone wondering about brain injuries, epilepsy, sleep or learning disorders, aging, etc. may want to seek further knowledge about how the brain works. Currently, the public is increasingly exposed to neurotechnologies due to the availability of consumer grade EEG devices, such as the Emotiv EPOC or the Neurosky Mind-Wave. Consequently, it has become necessary to design tools and user interfaces which will allow the general public to visualize, understand and interact with EEG signals. For instance, Mullen et al. proposed a software solution to process EEG signals collected using wearable EEG devices and visualize them in 3D [START_REF] Mullen | Real-time modeling and 3D visualization of source dynamics and connectivity using wearable EEG[END_REF]. This software enables the user to estimate brain activity sources and connectivity but is still mainly designed for brain signal and neuroscience experts, and not so public-friendly. Another recent work, more suited to lay persons, is the "Portable Brain Scanner" [START_REF] Stopczynski | The Smartphone Brain Scanner: A Portable Real-Time Neuroimaging System[END_REF]. This system makes use of a consumer-grade EEG device (the Emotiv EPOC) and a smartphone to provide a cheap and portable solution enabling anyone to visualize the sources of their brain activity on their smartphone in 3D. Another more attractive work, which is the most closely related to our Teegi system, is the "Mind-Mirror" system [START_REF] Mercier-Ganady | The Mind-Mirror: See Your Brain in Action in Your Head Using EEG and Augmented Reality[END_REF]. This system combines Augmented Reality (AR), 3D Visualization, and EEG to enable users to visualize their own brain activity in real-time superimposed to their own head, thanks to a semi-transparent mirror-based AR setup. This short review of the existing literature about making EEG accessible to the general public revealed that this is still a vastly unexplored area. Moreover, these solutions do not take into account any representation that the general public may have regarding the brain and EEG signals -many lay people do not even know what EEG signals are -in order to provide suitable visualizations and interaction devices to better understand these concepts. Some rare studies have indicated that misconceptions about brain functions prevail in general public [START_REF] Dekker | Neuromyths in education: Prevalence and predictors of misconceptions among teachers[END_REF][START_REF] Herculano-Houzel | Do you know your brain? A survey on public neuroscience literacy at the closing of the decade of the brain[END_REF][START_REF] Simons | What people believe about how memory works: A representative survey of the US population[END_REF]. These works stress the importance of popular scientific communication and indicate that communication efforts should be focused on increasing public awareness. It is important to note that the existing works mentioned above are mostly centered on visualization, with little or no interaction possibilities to manipulate and understand the EEG signals in real-time and in a friendly way. This further deters the general public from understanding brain activity [START_REF] Vosniadou | Knowledge Acquisition and Conceptual Change[END_REF]. Therefore, with the aim to enhance general public awareness, our work associates technical innovation and user-centered design.

INTRODUCING TEEGI

Founding principles

Design choices were made according to pedagogical principles. It has long been recognized that learner-centered education is much more effective than transmission-based education, even in informal situations [START_REF] Wellington | Formal and informal learning in science: The role of the interactive science centres[END_REF]. According to the constructivist paradigm, people create unique personal meanings by reflecting on interactive learning experiences. Therefore, people/learners should investigate and manipulate in order to become conscious of complex phenomena, change their misconceptions and construct scientific knowledge [START_REF] Vosniadou | Designing learning environments to promote conceptual change in science[END_REF]. In association, meaningful models play an important role in this type of learning processes [START_REF] Fleck | An Augmented Reality Environment for Astronomy Learning in Elementary Grades: An Exploratory Study[END_REF]. This motivated the design of an anthropomorphic interface that can be freely manipulated.

Our user-centered interactive media uses Spatial Augmented Reality (SAR) and tangible interaction. SAR, introduced by Raskar et al. [START_REF] Raskar | Shader Lamps: Animating Real Objects With Image-Based Illumination[END_REF], adds dynamic graphics to real-world surfaces by the means of projected light. Many systems were designed using projectors to add "painted" surface [START_REF] Raskar | Shader Lamps: Animating Real Objects With Image-Based Illumination[END_REF] or to give the illusion of virtual elements actually being there [START_REF] Wilson | Steerable augmented reality with the beamatron[END_REF][START_REF] Benko | Miragetable: freehand interaction on a projected augmented reality tabletop[END_REF]. A related approach is Tangible User Interface (TUI). TUI is concerned with providing tangible (i.e. physical) representations to digital information and controls [START_REF] Shaer | Tangible User Interfaces: Past, Present, and Future Directions[END_REF]. One of the strengths of tangibles is their situatedness: the interaction takes place in a real-world environment that often hides most of the technological aspects to expose physical interaction components only. They are particularly well suited for mediation purposes as they tend to be more inviting compared to mouse-screen based interfaces [START_REF] Horn | Comparing the use of tangible and graphical programming languages for informal science education[END_REF]. SAR and TUIs are often found together [START_REF] Underkoffler | Urp: a luminous-tangible workbench for urban planning and design[END_REF][START_REF] Piper | Illuminating clay: a 3-D tangible interface for landscape analysis[END_REF]. They are very complementary in that they both take place in the real world, in a common canvas. The tangible interaction allows for a hands-on approach by offering different input affordances (as well as physical constraints) to the user while the SAR technology allows for a flexible and situated way to give feedback. SAR can also be used as an affordable way to embed dynamic graphics on a physical surface that would otherwise require curved displays [START_REF] Brockmeyer | Designing Curved Display Surfaces with Printed Optics[END_REF] or rear projection [START_REF] Benko | Sphere: multi-touch interactions on a spherical display[END_REF].

There are examples of systems that use either tangible or AR principles to interact or review physiological data. Hinckley et al. [START_REF] Hinckley | Passive real-world interface props for neurosurgical visualization[END_REF] designed a system which used tangible props to do neurosurgical planning. A small tangible head was used in conjunction with a plastic plane to select the cutting planes to be visualized on a screen. Also mentioned above, the "Mind-Mirror" [START_REF] Mercier-Ganady | The Mind-Mirror: See Your Brain in Action in Your Head Using EEG and Augmented Reality[END_REF] is the work closest to Teegi. However, with Teegi, the data is not co-localized with the data source. It provides flexibility and easier visualization as the users can change viewpoints by tangible interactions instead of rotating their head while keeping their eyes on the mirror. This "out-of-body" visualization also enables collaboration where multiple users can explore the data.

General description

Teegi is a tangible interface that enables users to visualize and analyze a representation of their own brain activity recorded via an EEG system in real-time and displayed on a physical character. After some processing of the raw signals, a dedicated visualization is projected directly on top of the character. This character is tracked, which allows us to co-locate the projection with the character's head, at any time. Hence, the user can easily visualize a realistic modeling of the EEG signals in any part of the scalp by manipulating the character, while maintaining a good spatial topology of the observed data. Teegi was purposely given a child-like appearance, as well as animated eyes (also projected) that blink at the same time as the users do, in order to breathe life to the character and enhance attractiveness. Indeed blinking can be easily detected in electrodes neighboring the eyes.

Three different filters can be applied to the raw data (see the technical section for details) enabling users to investigate influences of motor motions, visual activities or meditation, on their brain activity in real-time. To remain consistent with the tangible philosophy of this project, we decided to control the filters by way of small tangible characters (mini-Teegis) that can be moved on a "filter area", which is highlighted on the table by a projected halo (see Figure 3). For example, if a user wants to apply a filter that will allow her to better see what happens when moving her hand, she just needs to take the dedicated mini-Teegi, i.e. the one with the colored hands, and to move it to the filter area. Then, by moving her right hand, she should see changes in EEG amplitude on the left hemisphere of Teegi's head, as illustrated in Figure 4. The manipulation of Teegi requires a motor activity. Therefore, when the motor filter is on, manipulating Teegi will obviously lead to observable changes in brain activity. At the end of this paper, we present an explorative study we conducted to obtain feedback about the main features of Teegi. However, Teegi is not limited to these first features. In the next section, we describe additional interaction metaphors we have explored, and that may benefit more advanced users. These advanced features were not evaluated during the study.

Advanced features

Visualizing the raw signal recorded on each electrode of the EEG is not very informative for the general public. However, this can be instructive for students who are learning EEG signal processing and analysis. In our approach, we can display on the table these raw data, as shown in Figure 5 (left). This creates a visual link between what is recorded with the EEG system, and the visualization that is provided on Teegi's head. This is possible because we know the rough position of the user, and the exact position of Teegi. When applying a filter, as described in the previous section, the user can see the effect of his or her action on the signal (see Figure 5, right). Compared to a standard approach where everything takes place on a screen, we believe that such a spatial and tangible approach might ease the understanding of the filters' effect. Another dimension we explored is the use of tangible actions to control some parameters of the EEG signal processing. As an example, we have implemented a technique where the user can control the amplitude of the visualization color map by moving a tangible object on the table (figure 6). This could be useful to reveal tiny fluctuations of EEG signals. With such interaction techniques, the whole analysis could be conducted without the use of a screen or a mouse, which remains consistent with the tangible philosophy of the project. Finally, we developed a solution that highlights the relationship between EEG signals and localized cortical sources, that is where the signals come from inside the brain. Using sLORETA inverse modeling [START_REF] Pascual-Marqui | Standardized low resolution brain electromagnetic tomography (sLORETA): technical details[END_REF] and Brainstorm to compute the kernel matrix [START_REF] Tadel | Brainstorm: a user-friendly application for MEG/EEG analysis[END_REF], we obtained a model of the cortex containing 2002 voxels linked to the 32 EEG electrodes we used. We can then project in real time the activity which arises from the outer regions of the cortex on an object representing the brain, alongside with Teegi (Figure 7). Since both Teegi and the brain proxy are tracked, it becomes possible to manipulate two synchronized representations of the same brain activity (the source at the surface of the brain and the measures on the scalp). This opens way to mediation activities that are more advanced all the while keeping the simplicity and ease of use brought forth by using SAR and tangible interaction.

TECHNICAL DESCRIPTION

EEG

We designed different EEG signal processing pipelines that each create a specific visualization tailored to identify specific elements in the signal. The details of these pipelines are transparent to the user. Each pipeline corresponds to a mini-Teegi filter. In particular, we set up the following EEG signal processing pipelines:

1. Wide-band EEG activity: EEG signals were band-pass filtered between 3 Hz and 26 Hz, in order to filter DC drift and part of the artifacts (e.g., facial muscle activity [START_REF] Fatourechi | EMG and EOG artifacts in brain computer interface systems: A survey[END_REF]) that may pollute them. Their power is then computed before being displayed. This corresponds to unspecific brain signals, hence they were labeled as "raw" signals. 2. Sensorimotor activity: EEG signals were first band-pass filtered in the β band (16-24Hz), a brain rhythm highly involved in sensorimotor tasks [START_REF] Pfurtscheller | Event-related EEG/MEG synchronization and desynchronization: basic principles[END_REF]. Then, they were spatially filtered, i.e., the signals from several neighboring EEG sensors were combined in order to enhance the signal of interest. In particular, we used and displayed Laplacian spatial filters around electrodes C3, C4 and Cz. This enabled the users to visualize EEG activity changes due to movements of the left hand, right hand and feet. Indeed, it is known that the power of EEG signals in the β rhythm decreases in electrodes C3/Cz/C4 during right hand/feet/left hand movements respectively, and increases just after the end of this movement [START_REF] Pfurtscheller | Event-related EEG/MEG synchronization and desynchronization: basic principles[END_REF]. 3. Visual activity: EEG signals were band-pass filtered in the α band (8-12 Hz), then only electrodes P3, Pz, P4, PO3, POz, PO4, O1, Oz and O2 (located on the back of the head, above the neck) were selected and displayed. These electrodes are indeed located over the visual cortex of the brain, i.e., the brain area in charge of visual information processing. The amplitude of the α rhythm is actually known to increase while the user is closing his/her eyes and is thus not processing any visual information [START_REF] Niedermeyer | Electroencephalography: basic principles, clinical applications, and related fields[END_REF]. To ensure that the user could perceive this increase after he/she reopened his/her eyes, the visualization was delayed by 0.5s. 4. Meditation: on a more exploratory note, we used the synchronization between the signals from the anterior and posterior cortex (AFz/Pz), which was measured in a 7-28 Hz band with instantaneous phase locking value [START_REF] Lachaux | Studying single-trials of phase synchronous activity in the brain[END_REF]. There are different outcomes (increase/decrease in synchronization) depending on meditation type. Mindfulness and body focus practices decrease the synchronization while transcendental practice increases it [START_REF] Lehmann | Reduced functional connectivity between cortical sources in five meditation traditions detected with lagged coherence using EEG tomography[END_REF].

EEG signals were acquired with a 32-channels EEG device (made of two g.tec g.USBAmp EEG amplifiers). This professional-grade system ensured that our prototype had a good signal-to-noise ratio and accurate electrode location, avoiding unneeded uncertainties. Signals were processed in real-time using OpenViBE [START_REF] Renard | An Open-Source Software Platform to Design, Test and Use Brain-Computer Interfaces in Real and Virtual Environments[END_REF]. For pipelines 1 to 3, the displayed colors correspond to signal power strength; for pipeline 4 they correspond to the degree of synchronization.

Spatial Augmented Reality

In order to create an augmented character, we have designed a tabletop augmentation setup (see Figure 8). Teegi itself is a 25cm high Trexi DIY toy. The mini-Teegis are also 10cm high Trexis. The main program handling the whole installation was created with vvvv [1]. The primary projected content (Teegi augmentation and GUI display) is handled with a single wide lens projector ProjectionDesign F20SX of resolution 1024x768 located over the table in a top-down orientation. The tracking of Teegi is achieved with an OptiTrack V120:Trio. It runs at 120 FPS with an overall latency of 8.3ms and a precision of 0.8mm. The OptiTrack is located in the same configuration as the main projector and both devices are calibrated together manually. The tracking data is sent to vvvv using OptiTrack's NatNet protocol. Teegi's eyes are projected using a second projector (Vivitek Qumi Q2) that is located on the side of the table.

The filter selection is done using a Sony PSEye web camera pointed at the position of the program selection GUI. Each mini-teegi representing a filter has a fiducial marker attached to it. The library ARToolkitPlus [START_REF] Wagner | Artoolkitplus for pose tracking on mobile devices[END_REF] is used to detect which marker is currently selected.

The EEG signals are processed by the OpenViBE software [START_REF] Renard | An Open-Source Software Platform to Design, Test and Use Brain-Computer Interfaces in Real and Virtual Environments[END_REF] that also generates a grayscale texture of the scalp signals. This texture is then exported to a local shared samba folder which is then fetched and remapped to an appropriate color scale in vvvv before being mapped to Teegi's head. In addition, the raw EEG signals are sent to vvvv over VRPN for display purposes (see Figure 5).

EXPLORATIVE STUDY AND DISCUSSIONS

Protocol

We conducted an explorative study where participants had to manipulate Teegi following a given scenario. The objectives of this study were to i) evaluate the general usability of the interface and ii) obtain initial feedback about the relevance of the approach to help users understand EEG signals and the brain. Ten participants (6 males, 4 females, mean age 28.6 (SD=9.7)) took part in this study. Pre-tests confirmed they were rather naive on the subject. They manipulated the version of Teegi described in the General Description section (no advanced features). The general procedure was as follows:

1. Pre-tests: The participant answered a first questionnaire assessing his or her representation of the brain. The participant then filled in different forms to measure his or her previous knowledge; one form per studied brain process (motor, vision and meditation). 2. Setting-up: The experimenter positioned the EEG cap on the participant's head. In parallel, the participant, guided by the experimenter, was made aware of the four didactic "cards" explaining the different filters i.e., Motor, Vision, Meditation and Raw. Each card was comprised of an image of the mini-Teegi associated with the filter along with basic instructions to follow (e.g. the Motor card indicated to the participants to move their hands or feet while staying relaxed). There were also two cards describing the two types of visualization participants could face, signal strength and synchronization. Once the participant was equipped, a quick calibration phase occurred. While Teegi was still inactive, participants were asked to close their eyes for a few seconds, and to move their hands and feet in order to identify the baseline activity for visualization. 3. Personal Investigation: The participant was asked to freely manipulate Teegi as well as the filters to be able to answer the following questions:

• What happens when you move your hands or feet?

• What happens when you close your eyes?

• What happens when you meditate?

During the whole study the participant sat comfortably in a chair. To avoid the occurrence of muscle artifacts that may pollute the signals, the user was instructed to stay relaxed and to refrain from making strong head movements. 4. Post-tests: The participant answered the questions above on dedicated forms, the same that were given at the beginning of step 1. Finally, he or she filled in a user survey questionnaire based on a 7-point Likert scale.

The whole session lasted approximately 1.5 hours per participant, with 15 to 20 minutes of hands-on time with Teegi. Each session was video-recorded. Video segments were separately visualized and labeled with the corresponding behavior (i.e. tangible and visual interactions, emotional expressions, and investigation strategies) using The Observer XT R 11.5 (Noldus, Info Tech, Wageninen, The Netherlands). After the session, the experimenter had an informal talk with the participant. He corrected the answers, making sure the participant was not leaving with false knowledge, and explained in more detail some aspects of the system (e.g. relationship between visual filter and attentional states, the various effects of meditation, ...). This phase lasted from 30 min to 1 hour depending on the participant's curiosity.

Results and discussions

To better understand the inherent strengths of Teegi towards learning, we assessed three main aspects of Teegi: its technical reliability, its relevance to ease understanding for nonexperts, and the User eXperience (UX) it provides. This evaluation is based on 1) the results of the questionnaire that are summarized in Figure 9, 2) the analysis of the video recordings, and 3) the analysis of the forms the participants filled in to assess their pre and post-knowledge of the brain and EEG.

Technical reliability

Participants unanimously reported that the whole system worked properly. The quality of the SAR display is valued by the participants. In particular, they reported that the resolution was appropriate, and they did not report problems of offset between the display and the physical character. Participants declared that they were not disturbed by occlusion problems. The mild temporal delay between their action and their consequences seems not to be an issue.

Manipulations of Teegi were numerous and frequent. Teegi was touched or moved on average 25% of the session's duration, twice per minute. These manipulations consisted mostly of rotations, and to a lesser extent of lifting Teegi to enhance visual perception. Two participants reported difficulties in grasping Teegi while the remaining 8 were comfortable with the form of the character. Video analyses did not show difficulties for the manipulation of Teegi. Similarly, applying filters by manipulating the mini-Teegis seemed easy for the participants.

Relevance of the interface to ease understanding

The participants reported that they understood the visualization associated with the filters. Video analyses indicated that they systematically used all filters several times (3 times per session on average) for a similar duration (Raw filter : 30.4% (SD 13.3) of session duration; Motor filter: 26.0% (8.3); visual filter: 16.9% (5.5); meditation filter: 26.6% (8.7)). Interestingly, the visual activity filter seemed slightly easier to understand than the other filters. Moreover, video analyses indicated that the participants did not have difficulty observing the signals on Teegi's head, as soon as they found the right location to observe. Overall, participants reported that they were able to use Teegi without any difficulties.

All participants completed the required tasks. They used instruction cards 5 times per session on average. They reported that they could focus on the tasks rather than on the mechanisms used to achieve them. This suggests that Teegi is a rather transparent interface. Regarding learning of brain processes and EEG, participants reported that they believed they had learned while doing the study. This was confirmed by the results of the pre-and post-test assessments (see Figure 10). These assessments focused on the recognition and the understanding of brain activation during Motor activities, Visual activities and Meditation. Understanding was marked as acquired if 1) the activated areas were correctly localized and 2) the explanations of the brain process were correct. It was marked as under way if only 1) or 2) was satisfied but not both; and as not acquired if neither 1) nor 2) were satisfied. The marks obtained by the participants improved after using Teegi. Overall, this suggests that Teegi offers many interesting features to ease learning and mediation.

All our results indicate that Teegi clearly promotes real-time tangible interactions, which contributes to enhancing awareness. Constructivism and inquiry-based science education principles indicate that, to become conscious of complex phe- nomena and construct scientific knowledge, people/learners have to experiment by interacting with and physically manipulating the content [START_REF] Vosniadou | Knowledge Acquisition and Conceptual Change[END_REF]. This is particularly true for brain activity that is difficult to understand because it cannot be sensed [START_REF] Damasio | Descartes' error: emotion, reason, and the human brain[END_REF], contrary to other activities (e.g. respiratory) that are perceived through sensory-motor mechanisms. Hence, brain activities need to be conceptualized, and the success of learning processes strongly depends on the interface. Teegi, which has been largely promoted by the participants, seems to fulfill this function.

User eXperience

The general experience with Teegi was rated as pleasant, attractive and stimulating, and participants did not feel stressed or oppressed. Overall, participants reported that they liked interacting with Teegi. The emotion expression analyses confirmed those statements. They showed that on average participants expressed curiosity and questioning about Teegi feedback during almost 20% (20.1% SD=9.1) of the manipulation duration. Other emotion expressions observed for all participants were joy and pleasure (e.g. smile, laugh, joyful verbal expression. . . ). They occurred during almost 10% (9.8% SD= 6.7) of the interaction duration with Teegi. Surprise emotions were observed but less frequently. Interestingly, boredom, weariness expressions occurred rarely (only for 2 users) and only at the end of the manipulation time. We did not observe any occurrence of exasperation or irritation. These results suggest a high level of acceptance for Teegi. This is a fundamental requirement for a tool aiming at improving access to knowledge.

Behavior observations indicated that the majority of participants spoke with Teegi and used morphological zones specific to human interactions while manipulating it. For example, they held its hands and held it up by the waist as one would do with a child. Some users spoke in the first person when they observed changes on the character's scalp for example "so, when I move my hands, I light up on the sides"; many said aloud that Teegi was their own image, for example "so, Teegi is me!". This identification suggests that an activation of associations between the perceived character's personality and self-perception may have occurred [START_REF] Paiva | Learning by Feeling: Evoking Empathy with Synthetic Characters[END_REF]. It is known that identification can be associated with increasing loss of self-awareness, and its temporary replacement with elements of the perceived character's personality [START_REF] Cohen | Defining Identification: A Theoretical Look at the Identification of Audiences With Media Characters[END_REF]. Therefore a human shaped, child-like character, made lifelike by animated projected eyes, could enhance both empathy and implicit selfperception of one's own brain activity, as provided by our interactive media. The anthropomorphic appearance of Teegi could explain the motivation and positive UX reported by the users. All these hypotheses would be the aim of a more extensive UX study.

Regarding visual attention, the participants were apparently paying attention to Teegi most of the time (83.3%, SD 7.6). This supports the fact that Teegi mobilized user attention. It also indicates a cognitive user engagement. Personal investigations were permanent (only 1.9% of inactivity was measured during the session duration; SD=1.7). Behavior analyses indicate that participants made predictions, hypotheses and tested them by conducting experiments. Numerous trial and error strategies were frequently used. This clearly indicates personal active control of the task and inquiry processes. Overall, Teegi stimulates investigations and encourages persistence in task completion.

CONCLUSIONS

In this paper, we presented Teegi, a tangible interface that makes EEG understandable to non-expert users. Our main contribution is the interface itself, which is built from both theoretical foundations, notably from human learning and scientific mediation and technical developments, including spatial augmented reality, tangible interaction and real-time neurotechnologies. We also demonstrated that this interface was well accepted by a first pool of users. In the future, we plan to make a more in-depth investigation into how well users are able to learn about EEG and brain activity with Teegi. To this end, we will conduct dedicated experiments with students and/or visitors in scientific museums. We would also like to precisely evaluate how Teegi benefits learning compared to standard approaches. For more advanced users, ad-hoc tangible filter creation could prove to be of great interest, adding flexibility to the overall system. Finally, it is known that BCI requires the user to learn to control his/her own brain activity to input computer commands [START_REF] Wolpaw | Principles and Practice[END_REF], which is a long and tedious task. We expect Teegi to be a motivating and informative way to support this training in the future.

Figure 1 .

 1 Figure 1. Teegi (Tangible EEG Interface) is a friendly interactive character that users can manipulate to observe and analyze their own brain activity in real-time.

Figure 2 (

 2 left and center) shows some typical visualizations of EEG signals used by experts, i.e., EEG signal traces and a 2D topographic map. More complex visualizations have been proposed, such as 3D topographic maps (Figure2, right), but they require many mouse inputs to be observed from all angles, which make them inconvenient to use in practice.

Figure 2 .

 2 Figure 2. (Left) A trace of EEG signals collected from multiple sensors. 2D (center) and 3D (right) topographic maps. (Screenshots from Open-ViBE [23]). The first two views are traditionally used by experts.

Figure 3 .

 3 Figure3. Three mini-Teegis can be used to apply high-level EEG filters to highlight brain processes associated to Motor, Vision and Meditation activities. To do so, the user simply needs to move the desired mini-Teegi into a specific zone projected on the table (green circle).

Figure 4 .

 4 Figure 4. Examples of the displayed visualizations on Teegi for each of the provided filters. Once a filter is active, the brain area corresponding to the selected and processed activity is highlighted in colors while the remaining EEG signals are displayed in grayscale.

Figure 5 .

 5 Figure 5. (Left) The raw EEG readings are displayed going from the user to the filter area and then rerouted towards Teegi. (Right) When a mini-Teegi (i.e. a filter) is active, the corresponding filtered signals are displayed between the filter area and Teegi instead.

Figure 6 .

 6 Figure 6. A moving tangible cursor is controlling the amplitude of the visualization color map.

Figure 7 .

 7 Figure 7. Using an inverse model, the cortical activity and EEG measures are presented together to users.

Figure 8 .

 8 Figure 8. Diagram of the installation. (A) ProjectionDesign F20SX projector (B) Sony PSEye web camera (C) OptiTrack V120:Trio (D) Vivitek Qumi Q2 projector (E) Teegi (F) Program selection zone and mini-Teegis.

Figure 9 .

 9 Figure 9. Results of the questionnaire (selected questions). Note that purple (resp. orange) bars indicate questions measuring Teegi's qualities (resp. limitations).

Figure 10 .

 10 Figure 10. Marks obtained by the participants during the pre-and posttest assessments. See text for details.


Jean-Baptiste Barreau
email: jean-baptiste.barreau@univ-rennes1.fr

Théophane Nicolas
email: theophane.nicolas@inrap.fr

Emilien Guillaume Bruniaux
email: guillaumebruniaux1989@gmail.com

Emilien Petit

Quentin Petit
email: quentin.petit@irisa.fr

Yann Bernard
email: yann.bernard@univ-rennes1.fr

Ronan Gaugne
email: ronan.gaugne@irisa.fr

Valérie Gouranton
email: valerie.gouranton@irisa.fr

Jean-Baptiste Barreau

J-B Barreau

Ceramics Fragments Digitization by Photogrammetry, Reconstructions and Applications

INTRODUCTION

Today, photogrammetry and 3D modelling are an integral part of the methods used in archeology and heritage management. They provide answers to scientific needs in the fields of conservation, preservation, restoration and mediation of architectural, archaeological and cultural heritage [START_REF] Barreau | The West Digital Conservatory of Archaeological Heritage project[END_REF] [START_REF] Arles | La photogrammétrie appliquée à l'archéologie minière : premier essai[END_REF] [7] [START_REF] Varea | Photogrammétrie et muséologie -Les sèvres du musée de Rouen[END_REF]. Photogrammetry on ceramic fragments was one of the first applications contemporary of the development of this technique applied in the archaeological community [START_REF] Kalantari | Implementation of a Low-Cost Photogrammetric Methodology for 3D Modelling of Ceramic Fragments[END_REF]. More recently and due to its democratization, it was applied more generally to artifacts [START_REF] Nicolae | Photogrammetry applied to problematic artefacts[END_REF]. Finally joined today by the rise of 3D printing [START_REF] Rutland | Non-contact 3D laser scanning as a tool to aid identification and interpretation of archaeological artefacts; the case of a Middle Bronze Age Hittite dice[END_REF] [10], it can restore fragmented artifacts [START_REF] Arbace | Innovative uses of 3D digital technologies to assist the restoration of a fragmented terracotta statue[END_REF] [START_REF] Powell | An ancient statue, re-created[END_REF]. These examples target one or several particular objects and use different types of equipment that can be expensive. These aspects can put off uninitiated archaeologists. So it would be appropriate to see if these techniques could be generalized to a whole class of geometrically simple and common artifacts, such as ceramics.

From these observations, associated to ceramics specialists with fragments of broken ceramics, we aimed at arranging different tools and methods, including photogrammetry, to explore opportunities for a cheap and attainable reconstruction methodology and its possible applications. Our first objective was to establish a protocol for scanning fragments with photogrammetry, and for reconstruction of original ceramics. We used the digital reconstitutions of the ceramics we got following our process to calculate some metrics and to design and 3D print a display for the remaining fragments of one pottery.

Archaeological context

The archaeological material studied in this work consists in domestic ceramics from an Iron Age site in Rezé (Loire-Atlantique -France) and from a Bronze Age site in Lannion Penn An Alé (Côte d'Armor, France, excavated by S. Blanchet, Inrap).

EQUIPMENT AND SOFTWARES

Equipment

As the study is based on a reconstitution of archaeological ceramics with photogrammetry, it was essential to use a high quality digital single-lens reflex camera with a good lens and good image resolution. Indeed, this allows to obtain very good quality photos and better results during the search for common points with the photogrammetry software. We used a Nikon D60 Digital SLR to take photos of the various ceramic fragments [START_REF] Labracherie | Test du reflex Nikon D60, Focus Numérique[END_REF]. Its technical features include a 10-megapixel CCD sensor, an 18-55mm optical AF-S with a 40mm focal distance and a Nikon autofocus module Multi-Cam 530 (3 AF points).

The environment used to take photos was composed of two light tables, a white background with its support and a halogen lamp (fig 1). This equipment allowed getting a light from below and above to minimize the shadow games. Homogeneity of brightness was difficult to obtain with only one single halogen lamp. This was a major problem because we needed a uniform and constant light regardless of the orientation of the shooting. Indeed, if the light intensity varies too much, then common points determined on a photo by photogrammetry software will not be the same from one photo to another. For 3D restitution, we used the photogrammetry software Agisoft PhotoScan version 0.8.3 beta 64bits. This software has been chosen for its ease of use and its ability to generate highly detailed 3D models. 3D reconstruction from scanned fragments was performed via Blender software. The scale and volume calculation were processed with MeshLab software.

3D DIGITIZATION

We present now the protocol followed to perform the 3D digitization of ceramic fragments by photogrammetry, the different steps to obtain realistic and satisfactory results, the technique used to calculate the volume of ceramics, and the problems we encountered during the work.

Shooting

The goal was here to shoot the object rotating on itself in order to get every viewing angle [START_REF] Kraus | Manuel de photogrammétrie -principes et procédés fondamentaux[END_REF]. We took a first series of photos, describing a semicircle, with a distance camera-object as regular as possible. Then, we rotated the object through 180° in order to complete the series and to get a first coverage of the object. Other series of photographs were taken from different angles to cover a maximum area (Fig. 2).

Point cloud

Photos of the object were then imported into PhotoScan. Before starting the search for common points, a mask constraint (clipping of the object) must be included on each photo to focus the calculation only on the object and not all the environment. The search for common points and pictures alignment is then performed with a "strong" accuracy. Scaling the object is usually performed by placing two markers whose distance is known, on a photo of the object. The distance between the markers is set in Photoscan as a reference distance that enables the software to generate a 3D model with a correct scale. This functionality requires a sufficiently uniform light, which was not the case in our installation. According to this limitation, we chose to handle the scaling of the model in another further step. The maximum number of common points was set to 60000 instead of 40000 (default) to get a better rendering. The point cloud generated allowed to verify the coherency of the digital output with respect to the object (fig 3).

Mesh and scale

The generation of geometries and textures was processed in PhotoScan with "arbitrary" option for the object type, "high" for the quality and "sharp" for geometric precision. The number of expected surfaces was set to 250000. The color of the mesh was computed from the color information contained in the points cloud. Photoscan computed the texture from several pieces of images extracted from the pictures, and applied it on the surface of the mesh (fig 4). The options for this step were a "generic" mapping mode and "average" fusion mode. The resulting mesh was then exported in the ".obj" format and imported in Meshlab to perform the scaling. A measure of distance was computed with the tool "rule" between two points whose real distance is known. The ratio between the real distance and the measured distance on the mesh provided the proportionality factor between these two distances. Thanks to the "Transform: scale" MeshLab function, distances on the X, Y and Z axis were multiplied by this factor achieving the correct scaling of the mesh. / textured meshs (bottom)

RECONSTRUCTIONS AND APPLICATIONS

Rezé ceramics

In the case of the Rezé ceramics, the available fragments cover only a partial part of each object. Therefore, the reconstruction was theoretical and performed by extrapolating the different shapes from the fragments in Blender. Most of the fragments contain a part of the lip, the neck, the edge, the body or the bottom. Therefore, if the fragment has a piece of ceramic bottom, it can be aligned with the horizontal plane XY for an easier handling. The lip of the ceramic fragment has a curvature. Using this curvature, a circle (circular mesh) was created and the perimeter of this circle adjusted to coincide perfectly with the contour of the lip. This provided the theoretical circumference of the upper ceramic part. The same operation was repeated several times by shifting the other circles on the Z-axis which were then adjusted on the other parts of the ceramic (kept and visible on the fragment). Therefore, the assembly of these circles on the fragment provided the "skeleton" of the ceramic. Once this frame done, it was imported in MeshLab again to calculate the volume of the ceramic. The "convex hull" function, in the "remeshing, simplification and reconstruction" tab allowed to create a mesh from the different circles of the frame (fig 5). This mesh, constituting a completely closed shape, was used to calculate a volume with the "compute geometric measures" function. The table 1 gathers the characteristics of the reconstitutions of the different Rezé ceramics.

Lannion Penn An Alé ceramic

This Bronze Age ceramic is "almost" full, so we tried to assemble the different fragments in Blender. Fragments of ceramic having a part of the lip were first assembled and wedged on a circle perfectly matching its shape. The other pieces were then imported and positioned one by one. One difficulty, however, was met with several small fragments. For them, only the inner or outer part of the fragment was reconstituted. The reasons for this failure are many and difficult to determine. Indeed, this can be induced by low quality of light, insufficient image quality or a too thin thickness of the ceramic fragments that corrupted the detection of common points on the edges. To solve this problem, the inner and outer parts were rescanned separately and then scaled in MeshLab. The two parts were then positioned parallel to the horizontal plane XY, superimposed on each other by playing with a rotation on the Z axis, and then merged by sliding the inner part relative to the Z axis.

The reconstruction approach was here different. After deleting the handles, we used the function "Surface reconstruction: Poisson" in Meshlab allowing to extrapolate the initial ceramic shape from the assembly of fragments (fig 6). This shape was used to compute the volume of the ceramic whose result is 2742,9 cm3. A last step of the work consisted in using the generated shape to build a support for the fragments, in order to get a 3D tangible representation of the pottery. We used a 3D printer MarkerBot Replicator 2x, which is an ABS printer, to produce the support display. The mesh required some corrections made by a graphic designer under 3DSmax (fig 7), because the Poisson reconstruction introduced some approximation at the frontiers between the real fragments and the computed shape.

DISCUSSION

The digitization process has highlighted difficulties often masked by the apparent simplicity of photogrammetry. However, we believe that the goal was reached regarding the digitization of fragments, and for a very low cost. The different approaches about reconstruction discussed here provided practical useful data for the analysis work performed by archaeologists. The estimate of their volumes, the 3D visualization of the shape, and the tangible 3D presentation represented concrete help in the study of the ceramics.

The 3D printer used for this work is a relatively inexpensive one, recently bought by the computer science laboratory involved in the project. The generation of a support display for the Lannion ceramic was a proof of concept for presenting fragments physically held around the content. However, problems have arisen. Because the Bronze Age ceramic is too large for this printer, especially because of the size of the handles, we had to print the support in two parts. In addition, the Poisson reconstruction of the internal volume of the ceramic generated approximations, which implied that the fragments did not exactly fit the display. Modifications on the 3d model of the display to match the shape of the fragments were required.

The result of this experimental work convinced us that the production such display could well be generalized within the museum community due to its low cost. A first reflection on its use has led us to believe that it would be wise to add on top of the display (not covered by a fragment) signs indicating metadata in the print-ready 3D model (fig 8).

CONCLUSION AND FUTURE WORKS

We proposed a simple protocol to digitize fragments of ceramics. This protocol is accurate enough to provide a digital reconstitution of ceramics with heterogeneous morphologies and compute basic metrics. Furthermore, from the digital reconstitution of the Rezé ceramics, we designed and 3D printed a display for the remaining fragments. This kind of display can be used for exhibition of artefacts in museums.

The next step of this work is to validate these methods on more ceramics. To do this, we must solve the problems of accuracy in digitization and printing of large volumes. We also wish to extend these methods to other kinds of less symmetrical archaeological fragmented artifacts. Finally, it would be interesting to integrate museologists to our future works on 3D printed displays.

Figure 1 .

 1 Figure 1. Shooting equipment and ceramic fragments 2.2 SoftwareFor 3D restitution, we used the photogrammetry software Agisoft PhotoScan version 0.8.3 beta 64bits. This software has been chosen for its ease of use and its ability to generate highly detailed 3D models. 3D reconstruction from scanned fragments was performed via Blender software. The scale and volume calculation were processed with MeshLab software.

Figure 2 .

 2 Figure 2. The different positions to take the series of photos. On the right, front view of the position of the series, and on the left, top view of the position of the camera in the 2 nd , 3 rd and 4 th series

Figure 3 .

 3 Figure 3. Piece of the Bronze Age ceramic with connection points (blue) determined by PhotoScan / point cloud get from the determination of connection points and photos alignment 3.3 Mesh and scaleThe generation of geometries and textures was processed in PhotoScan with "arbitrary" option for the object type, "high" for the quality and "sharp" for geometric precision. The number of expected surfaces was set to 250000. The color of the mesh was computed from the color information contained in the points cloud. Photoscan computed the texture from several pieces of images extracted from the pictures, and applied it on the surface of the mesh (fig4). The options for this step were a "generic" mapping mode and "average" fusion mode. The resulting mesh was then exported in the ".obj" format and imported in Meshlab to perform the scaling. A measure of distance was computed with the tool "rule" between two points whose

Figure 4 .

 4 Figure 4. Raw mesh (top left) / smoothed mesh (top middle) / mesh with colors (top right) / textured meshs (bottom)

Figure 5 .

 5 Figure 5. Fragments of eating ceramics reconstituted by photogrammetry / Frame of circles aligned along the Z axis and dimensioned according to the curvature of the fragments / "convex hull" function on the Gr8B ceramic

Figure 6 .

 6 Figure 6. Inner part (top) and outer part (bottom) of the bottom / reconstruction step / reconstruction with the function "Surface reconstruction: Poisson"

Figure 7 .

 7 Figure 7. 3D model of the display/ empty display / display with fragments

Figure 8 .

 8 Figure 8. Separated part of the display / metadata inscription on the 3d model

Table 1 .

 1 Summary of the characteristics of each reconstituted ceramic fragments

	Ceramic	Number of	Number of	Number	Number	Calculation	Volume
	name	photos	points (point	of	of	time	(cm3)
			cloud)	surfaces	vertices
	Gr1C	36	24101	249 999	130234	40m	536,4
	Gr5E	78	55470	246768	128077	1h14m	1722,8
	MGr2E	89	31204	249999	139674	52m	88,9
	MGr11D	75	18354	250000	129272	43m	138,8
	Gr1E	44	29396	17911	11036	39m	338,7
	Gr7D	62	45250	154286	77517	41m	211,1
	Gr8B	111	59623	250000	149449	1h15m	1938,7


Christophe Steiner

Michel Mehrenberger

Nicolas Crouseilles

Virginie Grandgirard

Guillaume Latu

Gyroaverage operator for a polar mesh

Keywords:

. In particular, we show that in the context of a drift-kinetic simulation, the proposed method has similar computational cost as the standard method and its precision is independent of the radius.

Introduction

In strongly magnetized plasma, when collision effects are negligible, one has to deal with kinetic models since fluid models, which assume that the distribution function is close to an equilibrium, are not suitable. However, the numerical solution of Vlasov type models is challenging since this model involves six dimensions in the phase space. Moreover, multi-scaled phenomena make the problem very difficult. Gyrokinetic theory enables to get rid of one of these constraints since the explicit dependence on the phase angle of the Vlasov equation is removed through gyrophase averaging while gyroradius effects are retained. The soobtained five-dimensional function is coupled with the Poisson equation (or its asymptotic counterpart, the so-called quasi-neutrality equation) which is defined on the particle coordinates. Thus, solving the gyrokinetic Vlasov-Poisson system requires an operator that transforms the gyro-center phase space in the particles phase space. This operator is the so-called gyroaverage operator. We refer to an abundant literature around this subject (see [START_REF] Brizard | [END_REF][START_REF] Hahm | Nonlinear Gyrokinetic Equations for Tokamaks Microturbulence[END_REF][START_REF] Lin | Method for solving the gyrokinetic Poisson equation in general geometry[END_REF] and references therein). The present work is devoted to the numerical computation of the gyroaverage operator for a polar mesh, following a a A part of this work was carried out using the HELIOS supercomputer system at Computational Simulation Centre of International Fusion Energy Research Centre (IFERC-CSC), Aomori, Japan, under the Broader Approach collaboration between Euratom and Japan, implemented by Fusion for Energy and JAEA.

previous work which was designed for cartesian geometry [START_REF] Crouseilles | Numerical solution of the gyroaverage operator for the finite gyroradius guiding-center model[END_REF]. We propose an alternative to the classical Padé approximation, that is especially employed in the gyrokinetic code GYSELA [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF][START_REF] Grandgirard | Computing ITG turbulence with a full-f semi-Lagrangian code[END_REF] and which is known to be only valid for small Larmor radius. The method is based on direct integration and interpolation. Close to the method used already in the GENE code (see [START_REF] Görler | The global version of the gyrokinetic turbulence code GENE[END_REF][START_REF] Görler | Multiscale effects in plasma microturbulence[END_REF]), it is applied in the context of 2D interpolation, as we consider here a fully global simulation.

The paper is organized as follows. In Section 2, we introduce the gyroaverage operator. The classical Padé approximation for the gyroaverage computation is described in Section 3. In Section 4, we present the method based on interpolation. A numerical comparison with analytical solutions is performed in Section 5. Finally, these gyroaverage operators will be applied to gyrokinetic simulations in Section 6.

Gyroaverage operator definition

Let ρ be the gyro-radius which is transverse to b = B/B (with B the magnetic field) and which depends on the gyrophase α ∈ [0, 2π], i.e ρ = ρ(cos(α)e ⊥1 + sin(α)e ⊥2 )

Here e ⊥1 and e ⊥2 are the unit vectors of a cartesian basis in the plane perpendicular to the magnetic field direction b. Let x G be the guiding-center radial coordinate and x the position of the particle in the real space. These two quantities differ by a Larmor radius ρ, i.e x = x G +ρ. The gyroaverage J ρ (f ) of any function f : (r, θ) → g(r cos(θ), r sin(θ)) depending on the spatial coordinates is defined as

J ρ (f )(r, θ) = 1 2π 2π 0 g(x + ρ)dα.
This gyro-average process consists in computing an average on the Larmor circle. It tends to damp any fluctuation which develops at sub-Larmor scales.

Introducing f (k) the Fourier transform of f , with k = k(cos(θ), sin(θ)) the wave vector, then the operation of gyro-average reads

J ρ (f )(x G ) = 2π 0 dα 2π +∞ -∞ d 3 k (2π) 3 f (k) exp{ik • (x G + ρ)} = +∞ -∞ d 3 k (2π) 3 2π 0 dα 2π exp(ik ⊥ ρ cos α) × f (k) exp(ik • x G )
where k ⊥ is the norm of the transverse component of the wave vector

k ⊥ = k -(b.k)b. Let consider J n the Bessel function of the first kind and of order n ∈ N, i.e ∀z ∈ C, J n (z) = i -n π π 0 exp(iz cos θ) cos(nθ) dθ.
Therefore, the previous gyro-average operation can be expressed as a function of the Bessel function of first order J 0 as

J ρ (f )(x G ) = +∞ -∞ d 3 k (2π) 3 J 0 (k ⊥ ρ) f (k)e ik•x G (1)
Let consider a uniform polar mesh (r, θ) ∈ [r min , r max ] × [0, 2π[ with N r × N θ cells, our goal is to approximate the operator

(f j,k ) ∈ R (Nr+1)×N θ → (J ρ (f ) j,k ) ∈ R (Nr+1)×N θ .
Considering expression [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF], in Fourier space the gyro-average reduces to the multiplication by the Bessel function of argument k ⊥ ρ. Indeed, the Fourier transform of J ρ (f ) can be written as

J ρ (f )(k) = 1 2π R 2 2π 0 f (x + ρ)dα e -ix•k dx = 1 2π 2π 0 R 2 f (x + ρ)e -i(x+ρ)•k dx e iρ•k dα = 1 2π 2π 0 e ikρ cos(α-θ) dα f (k)
which leads to

J ρ (f )(k) = J 0 (kρ) f (k) ( 2
)
This operation is straightforward in simple geometry with periodic boundary conditions, such as in local codes. Conversely, in the case of global codes, the use of Fourier transform is not applicable for two main reasons: (i) radial boundary conditions are non periodic, and (ii) the radial dependence of the Larmor radius has to be accounted for. Several approaches have been developed to overcome this difficulty. The most widespread method for this gyro-averaging process is to use a quadrature formula. In this context, the integral over the gyro-ring is usually approximated by a sum over four points on the gyro-ring [START_REF] Lee | Gyrokinetic approach in particle simulation[END_REF]. This is rigorously equivalent to considering the Taylor expansion of the Bessel function at order two in the small argument limit, namely

J 0 (k ⊥ ρ) ≃ 1-(k ⊥ ρ) 2 /4
, and equivalent to computing the transverse Laplacian at second order using finite differences. This method has been extended to achieve accuracy for large Larmor radius [START_REF] Hatzky | Energy conservation in a nonlinear gyrokinetic particle-in-cell code for ion-temperature-gradient-driven modes in theta-pinch geometry[END_REF],

i.e the number of points (starting with four) is linearly increased with the gyro-radius to guarantee the same number of points per arclength on the gyro-ring. In this approach -used e.g. in [START_REF] Idomura | Gyrokinetic theory of drift waves in negative shear tokamaks[END_REF] and [START_REF] Jolliet | A global collisionless PIC code in magnetic coordinates[END_REF]-the points that are equidistantly distributed over the ring are rotated for each particle (or marker) by a random angle calculated every time step. This is performed on a finite element formalism and enables therefore high order accuracy by keeping the matricial formulation.

In [START_REF] Crouseilles | Numerical solution of the gyroaverage operator for the finite gyroradius guiding-center model[END_REF] the influence of the interpolation operator (which is of great importance when the quadrature points do not coincide with the grid points) has been studied and has shown that the cubic splines are a good candidate. Some techniques used in [START_REF] Crouseilles | Numerical solution of the gyroaverage operator for the finite gyroradius guiding-center model[END_REF] taking advantages of the cartesian coordinates properties are no more valid in (r, θ) polar geometry. In this paper we present a method based on a direct integration of the gyroaverage operator which is directly applicable for global gyrokinetic code in toroidal geometry as GYSELA code. This new approach is tested for two different interpolation methods, one based on cubic splines and the other one on Hermite polynomials. Both are compared to the Padé approximation.

Method based on Padé approximation

One possible solution to compute the gyroaverage operator is to approximate the Bessel function with a Padé expansion J Pade (kρ) = 1/ 1 + (kρ) 2 /4 (e.g. see [START_REF] Sarazin | Kinetic features of interchange turbulence[END_REF]). As described in the following, such a Padé representation then requires the inversion of the Laplacian operator ∇ 2 ⊥ in real space. Indeed, using this approximation, the relation ( 2)

reads 1 + (kρ) 2 4 J ρ (f )(k) = f (k)
which corresponds in real space to

1 - ρ 2 4 ∆ J ρ (f )(x) = f (x).
We project f and J(f ) in the Fourier basis :

f (r, θ) ≈ N θ -1 n=0 A n (r)e inθ , J ρ (f )(r, θ) ≈ N θ -1 n=0 B n (r)e inθ .
The Laplacian in polar coordinates is expressed as

∆ = ∂ 2 ∂r 2 + 1 r ∂ ∂r + 1 r 2 ∂ 2 ∂θ 2
then we get for n = 0...N θ -1

- ρ 2 4 B ′′ n (r) - ρ 2 4r B ′ n (r) + 1 + ρ 2 n 2 4r 2 B n (r) = A n (r)
that can be solved by finite differences. For the results presented in the following, we use finite differences of second order which leads to a tridiagonal system. This Padé approximation gives the correct limit in the large wavelengths limit kρ ≪ 1, while keeping J Pade finite in the opposite limit kρ → ∞. The drawback is an overdamping of small scales: in the limit of large arguments x → ∞, J Pade (x) → 4/x 2 , whereas J 0 → (2/πx) 1/2 cos(xπ/4) (see figure 2). The method proposed in the next section which is no more based on an approximation of the Bessel function but on the direct calculation of the integral on a circle of radius ρ has been developed to overcome this drawback. Representation in the real space (figures from [START_REF] Sarazin | Kinetic features of interchange turbulence[END_REF]).

Method based on interpolation

We put N uniformly distributed points on the circle of integration and we approximate the function value at these points by interpolation. The gyroaverage is then obtained by the rectangle quadrature formula on these points. More precisely, for a given point (r j , θ k ), the gyroaverage at this point is approximated by

J ρ (f ) j,k ≃ 1 2π N -1 ℓ=0 P(f )(r j cos θ k + ρ cos α ℓ , r j sin θ k + ρ sin α ℓ )∆α,
where α ℓ = ℓ∆α and ∆α = 2π/N . Since the quadrature points do not coincide with grid points, we introduce an interpolation operator P which can be -Hermite interpolation, -Cubic splines interpolation.

As detailed in [START_REF] Mehrenberger | Vlasov on GPU[END_REF], the interpolation can be reformulated into a matrix-vector product

P(f )(r j cos θ k + ρ cos α ℓ , r j sin θ k + ρ sin α ℓ ) = (A ℓ c) j,k ,
where c denotes the splines coefficient (cubic splines method) or the function values (Hermite method) so that the gyroaverage can be itself viewed as a matrix-vector product

J ρ (f ) j,k = 1 2π N -1 ℓ=0 (A ℓ f ) j,k ∆α = (A ρ c) j,k .
As a consequence, for a given Larmor radius ρ, the matrix A ρ can be stored once for all.

•

For each method, 2 versions are implemented :

a basic version a version with precomputation where we first compute the matrix A ρ such that

(J ρ (f ) j,k ) = A ρ c,
where c are spline coefficients (vector of size (N r + 1)N θ ) or function values and derivatives in the case of Hermite interpolation (size is 4(N r + 1)N θ ). Remark 1 Note that for the Hermite interpolation, we need first to compute the derivatives at each cell interface. These derivatives are reconstructed by centered finite differences of arbitrary even order (we could also use odd order by having only a C 0 reconstruction, but then the size of c would increase to 9(N r + 1)N θ ). In the numerical results, we will take the order 4.

Remark 2 Time comparison are given in Tables 1 and2.

The use of the precomputation version is here quite efficient, as the Larmor radius ρ is fixed and the matrix is the same for each value of θ, which implies that the storage is reduced; but the basic version permits to give a rough indication of the time that would be used for more general situations (where for example the storage would be an issue), that are not considered for the moment. Optimization strategies, in a parallel environment, may be the subject of further extensions.

Remark 3

We may wonder about the numerical cost of the Hermite method in comparison to the Padé approximation, especially for a large radius. Luckily, the picture will change in the 4D case, since the gyroaverage is applied only in 3D. As an example, we have obtained the following times in the case of a drift kinetic simulation (see subsection 6): 42s for PADE, 43s for Hermite interpolation with precomputation (using 1024 quadrature points). Without precompution, the time for Hermite is 47s with 16 quadrature points and 270s with 1024 quadrature points. Computations are made on a local cluster of the University of Strasbourg using 16 processors (grid 32 × 32 × 32 × 64, 100 iterations). Further time measures will be detailed in subsection 6.

Proposition. Let z ∈ C. The gyroaverage of f (r, θ) = C m (zr)e imθ reads J ρ (f )(r 0 , θ 0 ) = J 0 (zρ)C m (zr 0 )e imθ0 .
Proof : By definition,

J ρ (f )(r, θ) = 1 2π 2π 0 +∞ 0 2π 0
C m (zr 0 )e imθ0 δ {x0=x+ρ} dαdr 0 dθ 0 where x 0 = r 0 (cos(θ 0 ), sin(θ 0 )), x = r 0 (cos(θ), sin(θ)) and ρ = ρ(cos(α), sin(α)).

The additivity theorem of Graf for Bessel functions (see [START_REF] Abramowitz | Handbook of Mathematical Functions[END_REF]) states that if u, v and w are the lengths of a triangle and α, γ the angles as shown in the following figure :

u w v α γ
then for all integers m and all complex numbers z,

C m (zw)e imγ = ∞ k=-∞ C m+k (zu)J k (zv)e ikα .
Then, we obtain, with v = ρ, w = r 0 , u = r, γ = θ 0 -θ and α = α,

J ρ (f )(r, θ) = e imθ 1 2π 2π 0 +∞ 0 2π 0 C m (zr 0 )e im(θ0-θ) × δ {x0=x+ρ} dαdr 0 dθ 0 = e imθ 1 2π 2π 0 +∞ 0 2π 0 ∞ k=-∞ C m+k (zr)J k (zρ)e ikα × δ {x0=x+ρ} dαdr 0 dθ 0 = e imθ ∞ k=-∞ C m+k (zr)J k (zρ) × 1 2π 2π 0 +∞ 0 2π 0 e ikα δ {x0=x+ρ} dαdr 0 dθ 0 = e imθ ∞ k=-∞ C m+k (zr)J k (zρ) × 1 2π 2π 0 e ikα dα.
We use the fact that

1 2π 2π 0 e ikα dα = δ k,0
in order to conclude that

J ρ (f )(r, θ) = J 0 (zρ)C m (zr)e imθ .
In the following, we give some examples of these test functions depending on the boundary conditions we want to test.

Examples

1. r min = 0 and homogeneous Dirichlet condition on r max .

Here we consider a disk and the function

f 1 (r, θ) = J m r j m,ℓ r max e imθ
where j m,ℓ is the ℓ th zero of J m verifies

f 1 (r max , θ) = 0, 0 ≤ θ < 2π,
and its gyroaverage reads

J ρ (f 1 )(r 0 , θ 0 ) = J 0 ρ j m,ℓ r max f 1 (r 0 , θ 0 ).
2. Homogeneous Dirichlet condition on r min > 0 and r max .

The following function is defined on an annulus

f 2 (r, θ) = J m (γ m,ℓ )Y m r γ m,ℓ r max - Y m (γ m,ℓ )J m r γ m,ℓ r max e imθ
where γ m,ℓ is the ℓ th zero of

y → J m (y)Y m y r min r max -Y m (y)J m y r min r max verifies f 2 (r min , θ) = 0, f 2 (r max , θ) = 0, 0 ≤ θ < 2π,
and its gyroaverage reads

J ρ (f 2 )(r 0 , θ 0 ) = J 0 ρ γ m,ℓ r max f 2 (r 0 , θ 0 ).
3. Homogeneous Neumann condition on r min > 0 and r max .

The following function is defined on an annulus

f 3 (r, θ) = J ′ m (η m,ℓ )Y m r η m,ℓ r max - Y ′ m (η m,ℓ )J m r η m,ℓ r max e imθ
where η m,ℓ is the ℓ th zero of

y → J ′ m (y)Y ′ m y r min r max -Y ′ m (y)J ′ m y r min r max verifies ∂ r f 3 (r min , θ) = 0, ∂ r f 3 (r max , θ) = 0, 0 ≤ θ < 2π,
and its gyroaverage reads

J ρ (f 3 )(r 0 , θ 0 ) = J 0 ρ η m,ℓ r max f 3 (r 0 , θ 0 ).
We have here used the fact that for C n = J n or Y n , the derivative reads :

C ′ n (r) = -C n+1 (r) + nC n (r) r .
We show in Fig. 3 the real and imaginary parts of the function

(r, θ) ∈ [0, 5] × [0, 2π] → J 1 (r)e iθ .

Numerical results

In this part, the different numerical methods are compared in the case of the second test case (Homogenous Dirichlet condition on r min and r max ) with r min = 0.1, r max = 0.9, N r = N θ = 512, ℓ = 1 and m = 1, 5, 20. In Tables 3 and4, we give here the L 2 -norm error of the gyroaverage function for m = 1 whereas Tables 5,6 refer to m = 5 and Tables 7,8 refer to m = 20. In Table 8, we use various orders of interpolation for Hermite. Note that for each method, the error is computed over the domain [r min + ρ, r max -ρ]. Remark 4 1. For ρ = 0, all the methods are exact.

2. The method based on interpolation gives almost the same results with Hermite interpolation or cubic splines interpolation.

The method based on Padé approximation gives very

good results for small values of ρ. 4. For large values of ρ, the method based on interpolation gives better results even with a relatively small number of points on the circle (for example with ρ = 0.1 and N ≥ 8).

Application to gyrokinetic simulations

The computational effort to numerically solve the 6 dimensional Vlasov-Maxwell systems which describes plasma turbulence in tokamak plasmas still remains out of reach for present day supercomputers. All the numerical simulations performed until now in this domain take care of the gyrokinetic ordering to reduce this problem of one dimension. This ordering take into account the fact that (i) electromagnetic fluctuations occur on time scales much longer than charged particle gyration period (ω/Ω c ≪ 1 with ω the fluctuation frequency and Ω c the cyclotron frequency), and (ii) the wavelength of these fluctuations is much smaller than the characteristic scale length of the gradients of magnetic field, density and temperature. See [START_REF] Garbet | Gyrokinetic simulations of turbulent transport[END_REF] for a detailed review on the gyrokinetic framework and simulations to compute turbulent transport in fusion plasmas. Within this gyro-ordering, the so-called gyrokinetic model can be derived (see [START_REF] Lee | Gyrokinetic approach in particle simulation[END_REF]) by averaging on the fast gyration of charged particles around the magnetic field lines. The magnetic toroidal configuration considered in this paper is simplified. Indeed, magnetic flux surfaces are assumed concentric torii with circular cross-sections. The gyroaverage operator described in section 2 occurs in this reduction from 6 to 5 dimensions. The new 5D set of coordinates corresponds to: (i) 3D toroidal spatial coordinates (r, θ, ϕ) (with r the radial direction, θ and ϕ the poloidal (resp. toroidal) angle), and (ii) 2D in velocity space with v the velocity parallel to the magnetic field line and µ = mv 2 ⊥ /(2B) the magnetic moment where v ⊥ represents the velocity in the plane orthogonal to the magnetic field. It is important to note that in this ordering µ is an adiabatic invariant, so it plays the role of a parameter in the 5D gyrokinetic Vlasov equation.

In the following, the 4D problem which is treated in section 6.1 corresponds to the case were we consider a unique value of µ, i.e the same Larmor radius is taken for all particles. The 5D problem, we deal with in section 6.2, several values of µ considered to take into account the dependence on the Larmor radius with v ⊥ .

The time evolution of the gyrocenter distribution function F is given by the gyrokinetic conservative equation (see also Eqs ( 17)-( 20) in [START_REF] Garbet | Gyrokinetic simulations of turbulent transport[END_REF]):

B * ∂ F ∂t + ∇ ∇ ∇ • B * dx G dt F + ∂: ∂v G B * dv G dt F = 0 (3)
where x G and v G are respectively the space coordinates and the parallel velocity of the guiding centers. In the electrostatic limit, for a particle of mass m and charge q the motion equations of the guiding centers are given by

dx G dt = v G b * + v E×B + v D (4)
m dv G dt = -µ∇ ∇ ∇ * B -q∇ ∇ ∇ * Φ + mv G v E×B • ∇ ∇ ∇B B ( 5
)
where ∇ ∇ ∇ * ≡ b * • ∇ ∇ ∇, while b * and B * are defined by:

b * ≡ B B * + mv G qB * B ∇ ∇ ∇ × B (6)
B * ≡ B + mv G qB b • (∇ ∇ ∇ × B) (7)
The 'E × B' drift is equal to

v E×B = (1/B * )b ×∇ ∇ ∇ Φ while curvature drift is defined as v D = mv 2 G +µB qB * b × ∇ ∇ ∇B B .
We focus on the turbulent transport driven by collisionless ITG instability so electrons are assumed adiabatic. In this limit, the gyro-averaged electrostatic potential Φ (equivalent to J √ 2µ Φ notation) is solution of the self-consistently coupled 3D quasi-neutrality equation:

-∂ 2 r Φ + 1 r + ∂ r n 0 (r) n 0 (r) ∂ r Φ + 1 r 2 ∂ 2 θ Φ + 1 T e (r) (Φ -λ Φ ) = 1 n 0 (r) J √ 2µ f -f eq dv , (8)
Φ = 1 L L 0 Φ(r, θ, z)dz, ( 9
)
where T e and n 0 are electron temperature and density profiles which will be defined latter.

The new gyroaverage methods presented before have been tested with two codes: (i) the SELALIB plateform [START_REF] Selalib | [END_REF] for the 4D simplified case and (ii) the GYSELA code [START_REF] Grandgirard | Computing ITG turbulence with a full-f semi-Lagrangian code[END_REF] for the benchmark with the classical Cyclone DIID 5D case. Both are based on a classical Backward semi-Lagrangian scheme (BSL) with cubic splines interpolation, predictor-corrector method.

In the following numerical solutions are computed using normalized equations. The temperature is normalized to T e0 , where T e0 is defined by the initial temperature profile such that T e (r p )/T e0 = 1. The time is normalized to the inverse of the ion cyclotron frequency ω c = e i B 0 /m i . Velocities, including the parallel velocity, are expressed in units of the ion speed v T 0 = T e0 /m i , the electric potential is normalized to T e0 /e i and the magnetic field is normalized to B 0 . Consequently, lengths are normalized to the Larmor radius ρ = m i v T 0 /e i B 0 and the magnetic moment µ to T e0 /B 0 .

Simplified 4D SLAB case

In this section, we consider a simplified model of the system of equations ( 3)- [START_REF] Grandgirard | Computing ITG turbulence with a full-f semi-Lagrangian code[END_REF]. A periodic cylindrical plasma of radius a and length 2πR (with R the major radius) is considered as a limit case of a stretched torus. The plasma is confined by a strong magnetic which is uniform B = Be z where e z stand for the unit vector in the toroidal direction z. With theses assumptions the velocity drifts are reduced to the E × B drift. This SLAB 4D case is equivalent to the one treated in [START_REF] Klein | Gyrokinetic water-bag modeling of a plasma column : Magnetic moment distribution and finite Larmor radius effects[END_REF] or [START_REF] Grandgirard | A drift-kinetic Semi-Lagrangian 4D code for ion turbulence simulation[END_REF].

The equation satisfied by the distribution function of ions f (t, r, θ, z, v) following the guiding center movement reads :

∂ t f - ∂ θ J √ 2µ Φ r ∂ r f + ∂ r J √ 2µ Φ r ∂ θ f + v∂ z f -∂ z J √ 2µ Φ ∂ v f = 0. (10) for (r, θ, z, v) ∈ [r min , r max ]×[0, 2π]×[0, L]×[-v max , v max ].
To deals with this equation system, we have used the SELALIB platform [START_REF] Selalib | [END_REF] with a classical semi-Lagrangian method with cubic splines interpolation, predictor corrector method and Verlet algorithm for the characteristics (see [START_REF] Crouseilles | Semi-Lagrangian simulations on polar grids: from diocotron instability to ITG turbulence[END_REF] for details). The platform has been improved by adding the one fixed µ capability and by implementing the three different gyroaverage operators described in sections 3 and 4. In our case the MPI parallelization is based on transpositions between (r, θ, v) domain decomposition and z domain decomposition. In this section the numerical instability growth rates are compared to the one deduced from the dispersion relation obtained by linearizing the self-consistent equation system ( 9)-(10).

Dispersion relation

In order to validate the linear part of the numerical results, we compute the dispersion relation with the gyroaverage. We make the following expansions:

f = f 0 + εf 1 + O(ε 2 ), φ = φ 0 + εφ 1 + O(ε 2 ) with f 0 (r, v) = f eq (r, v) = n 0 (r) exp -v 2 2Ti(r) (2πT i (r)) 1/2 , φ 0 = 0.
Then we obtain

J √ 2µ (f ) = J √ 2µ (f 0 + εf 1 ) + O(ε 2 ), and J √ 2µ (φ) = J √ 2µ ( φ0 + ε φ1 ) + O(ε 2
). Substituting the above relations into [START_REF] Grandgirard | Computing ITG turbulence with a full-f semi-Lagrangian code[END_REF], we obtain

∂ t f 1 - ∂ θ J √ 2µ (φ 1 ) r ∂ r f 0 + v∂ z f 1 -∂ z J √ 2µ (φ 1 )∂ v f 0 = O(ε).
Similarly, the equation

-∂ 2 r φ + 1 r + ∂ r n 0 (r) n 0 (r) ∂ r φ + 1 r 2 ∂ 2 θ φ + 1 T e (r) (φ -λ φ ) = 1 n 0 (r) J √ 2µ f -f eq dv becomes -∂ 2 r φ 1 + 1 r + ∂ r n 0 (r) n 0 (r) ∂ r φ 1 + 1 r 2 ∂ 2 θ φ 1 + 1 T e (r) (φ 1 -λ φ 1 ) = 1 n 0 (r) J √ 2µ f 1 dv + O(ε).( 11
)
We assume that the solutions have the form :

f 1 = f m,n,ω (r, v)e i(mθ+kz-ωt) , φ 1 = φ m,n,ω (r)e i(mθ+kz-ωt)
J √ 2µ (f 1 ) = fm,n,ω (r, v)e i(mθ+kz-ωt) , J √ 2µ (φ 1 ) = φm,n,ω (r)e i(mθ+kz-ωt)

with k = 2πn L . Then, we obtain

(-ω + kv)f m,n,ω = m r ∂ r f 0 + k∂ v f 0 φm,n,ω (12)
and the relation [START_REF] Görler | Multiscale effects in plasma microturbulence[END_REF] becomes

-∂ 2 r φ m,n,ω + 1 r + ∂ r n 0 (r) n 0 (r) ∂ r φ m,n,ω - m 2 r 2 φ m,n,ω + 1 T e (r) (φ m,n,ω -λδ 0 n φ m,0,ω ) = 1 n 0 (r) fm,n,ω dv.
If we assume that m = 0 and n = 0, the last relation and the equation ( 12) lead to :

-∂ 2 r φ m,n,ω + 1 r + ∂ r n 0 (r) n 0 (r) ∂ r φ m,n,ω - m 2 r 2 φ m,n,ω + 1 T e (r) φ m,n,ω = 1 n 0 (r) φm,n,ω fm,n,ω f m,n,ω m r ∂ r f 0 + k∂ v f 0 kv -ω dv and then - ∂ 2 r φ m,n,ω φ m,n,ω + 1 r + ∂ r n 0 (r) n 0 (r) ∂ r φ m,n,ω φ m,n,ω - m 2 r 2 + 1 T e (r) = 1 n 0 (r) φm,n,ω φ m,n,ω fm,n,ω f m,n,ω m r ∂ r f 0 + k∂ v f 0 kv -ω dv.
We make the approximations :

Φm,n,ω Φ m,n,ω ≈ J 0 (κ 2µ), fm,n,ω f m,n,ω ≈ J 0 (κ 2µ)
where κ ∈ R + . Rigorously the previous approximations are true when f and Φ are Fourier-Bessel functions in (r, θ), i.e. when

f m,n,ω (r, v) = J m (κr) × g(v)
and Φ m,n,ω (r) = J m (κr) (see Proposition in Section 3). In general, we are not in this case and this will explain the differences we observe in Figure 7. Then, by considering the previous approximations, we obtain :

- ∂ 2 r φ m,n,ω φ m,n,ω + 1 r + ∂ r n 0 (r) n 0 (r) ∂ r φ m,n,ω φ m,n,ω - m 2 r 2 + 1 T e (r) = J 0 (κ 2µ) 2 1 n 0 (r) m r ∂ r f 0 + k∂ v f 0 kv -ω dv.
By setting

I = m r ∂ r f 0 + k∂ v f 0 kv -ω dv
and by using the expression of f 0 , we have

I = -v Ti + m kr ( ∂rn0 n0 -∂rTi 2Ti + v 2 ∂rTi 2T 2 i ) v -ω k f 0 dv.
Now, we introduce for n ∈ N :

I n = 1 n 0 v n f 0 v -ω k f 0 dv
and we obtain the relations :

I 1 = 1 + ω k I 0 , I 2 = ω k 1 + ω k I 0 .
By using the change of variable v = (2T i (r)) 1/2 w and the expression of f 0 , we have by setting k * = (2T i ) 1/2 k.

I 0 = exp -v 2 2Ti (2πT i ) 1/2 (v -ω k )dv = exp(-ω) π 1/2 ((2T i (r)) 1/2 w -ω k )dv = 1 (2T i ) 1/2 Z ω k * with Z(z) = 1 √ π exp(-x 2 ) x -z dx = i √ π exp(-z 2 )(1 -erf(-iz)), erf(x) = 2 √ π x 0 exp(-t 2 )dt.
Finally, for z = ω/k * :

J 0 (κ √ 2µ) 2 n 0 (r) I = - 1 T i (1 + zZ(z)) + m k * r Z(z) ∂ r n 0 n 0 - ∂ r T i 2T i + z(1 + zZ(z)) ∂ r T i T i . (13)

Instability growth rate comparisons

We adapt a code available in SELALIB, that computes the zeros for the dispersion relation [START_REF] Hahm | Nonlinear Gyrokinetic Equations for Tokamaks Microturbulence[END_REF], as in [START_REF] Coulette | Numerical comparisons of gyrokinetic multi-water-bag models[END_REF], by adding the gyroaverage term. Figure 7 presents the instability rates as a function of µ. We obtain the two first curves by solving the dispersion relation with J 0 ( √ 2µ) 2

(curve in red) or by substituting J 0 ( √ 2µ) by its Padé approximation (curve in green). We have chosen κ = 1, in [START_REF] Hahm | Nonlinear Gyrokinetic Equations for Tokamaks Microturbulence[END_REF]. The two remaining curves are obtain numerically with the Padé method for the gyroaverage operator (curve in blue) or the method with Hermite interpolation (curve in magenta). It appears that the slope decreases faster with the Hermite interpolation method than with the Padé method. The slopes obtained with the numerical Padé are different from these obtained with the dispersion relation and the Padé approximation because the functions that we consider here are not Fourier-Bessel functions.

In the simulations, we will take λ = 0 (no zonal flow case). The initial distribution function reads :

f (0, r, θ, z, v) = f eq (r, v) × 1 + ε exp - (r -r p ) 2 δr cos 2πn L z + mθ
where the equilibrium function f eq is

f eq (r, v) = n 0 (r) exp -v 2 2Ti(r) (2πT i (r)) 1/2 .
The profiles T i , T e and n 0 are given by : P Numerical results are given in Fig. 4 -8. We will consider here N = 1024 for Hermite, wiht precomputation. For small µ, the differences between Hermite and Padé are very small and they are comparable to the case without gyroaverage. When µ increases, the differences between the two methods increase and the structures become coarser. The gyroaverage tends to reduce the instability rate; the more µ is large, the more this rate is small (Fig. 7). Time results (Tab. 10 and 11) show that the choice of the gyroaverage operator is not very influential in the total time. Indeed, the computation of the gyroaverage appears to be a 3D problem in a 4D environment.

Benchmark with the classical 5D Cyclone DIII-D case

For this part, the gyroaverage operator based on cubic spline and Hermite interpolation have been implemented in the GYSELA code [START_REF] Grandgirard | Computing ITG turbulence with a full-f semi-Lagrangian code[END_REF] and compared to the existing Padé approximation. The 5D Vlasov-Poisson problem considered is the one described by equations ( 3)-( 9). The numerical comparisons have been performed on a linear benchmark based on the classical cyclone DIII-D case [START_REF] Dimits | Comparisons and physics basis of tokamak transport models and turbulence simulations[END_REF]. This kind of typical benchmark had already been performed several years ago to validate the GYSELA code [START_REF] Grandgirard | Computing ITG turbulence with a full-f semi-Lagrangian code[END_REF]. For the present tests, the same parameters than in section 4 of [START_REF] Grandgirard | Computing ITG turbulence with a full-f semi-Lagrangian code[END_REF] have been used. The comparison between the Padé operator and the gyroaverage operator based on Hermite interpolation is presented in the following tabular for four unstable modes (m, n) with m the poloidal mode number and n the toroidal one. The results obtained with cubic spline interpolation are not detailed because very similar to the results obtained with the Hermite interpolation. For the simulations, we have used N = 32 quadrature points, with the Hermite method and without precomputations. We observe the same behaviour as in the previous testcase: the instability growth rate is smaller (except for the first point), with Hermite, and the values are in the same range.

Conclusion

We approximation, considering on the one hand analytical cases, for which we know the exact solution, and on the other hand some basic gyrokinetic simulations: a 4D driftkinetic model with one Larmor radius and the classical linear DIII-D benchmark case. We find that the introduction of the gyroaverage operation tends to diminish the growth rate of the instability and this is amplified, when considering the direct (right) gyroaverage operator, instead of the Padé approximation. Linear analysis predicts a similar behaviour, when we compare Padé approximation and the J 0 Bessel function. This is coherent as the Padé approximation is above the Bessel function for kρ relatively small (see on Figure 2). Note that the result remains at the qualitative level, as here, in polar geometry, the multiplication by the J 0 Bessel function is not the exact solution (expect for Fourier-Bessel functions), and this differs from the cartesian geometry.

Fig. 1 .Fig. 2 .

 12 Fig. 1. The zero-th order Bessel function J0(kρ) compare to its Padé approximation 1/ 1 + (kρ) 2 /4 .

5

  Numerical comparison with analytical solutions 5.1 Definition of a class of analytical solutions depending on boundary conditions First, we give a family of functions whose gyroaverage is analytically known. For these functions, we obtain the gyroaverage just by multiplying them by the Bessel function. Let m ≥ 0 be an integer and C m be the Bessel function of the first kind (denoted by J m ) or the Bessel function of the second kind (denoted by Y m ). The following proposition gives the analytic expression of the gyroaverage of Fourier-Bessel type functions.

8 Table 3 . 8 Table 4 .

 8384 Comparison Padé with Hermite (m = 1). ρ Padé N = 4 N = 8 N = 16 N Comparison Padé with splines (m = 1). ρ Padé N = 4 N = 8 N = 16 N =

Fig. 3 .

 3 Fig. 3. Real and imaginary parts of the test function (r, θ) → J1(r) exp(iθ).

dr

  r) = C P exp -κ P δr P tanh r -r p δr P where P ∈ {T i , T e , n 0 }, C Ti = C Te = 1 and C n0 = r max -r min rmin rmax exp -κ n0 δr n0 tanh r-rp δrn 0 We consider the parameters of [3] [Medium case] : r min = 0.1, r max = 14.5, v max = 7.32, κ n0 = 0.055, κ Ti = κ Te = 0.27586, δr Ti = δr Ti = δr n0 2 = 1.45, ε = 10 -6 , n = 1, m = 5, L = 1506.759067, r p = r min + r max 2 , δ r = 4δr n0 δr Ti .

Fig. 4 .

 4 Fig. 4. Poloidal cut f (r, θ, 0, 0) at time T = 7000 for 128 × 256×128×128, ∆t = 2, successively from top to bottom µ = 0.5 with Hermite and then Padé; µ = 1, with Hermite and then Padé.

Fig. 5 .

 5 Fig. 5. Poloidal cut f (r, θ, 0, 0) at time T = 5000 for 128 × 128×128×128, ∆t = 1, successively from top to bottom µ = 0.1 with Hermite and then Padé; µ = 0, with Hermite and then Padé.

Fig. 6 .

 6 Fig. 6. Poloidal cut f (r, θ, 0, 0) at time T = 5000 for 64 × 64 × 32 × 64, ∆t = 5 (left, Hermite; right, Padé) for µ = 0.1, . . . , 0.8 (from top to bottom).

Fig. 7 .Fig. 8 .

 78 Fig. 7. Time evolution of rmax r min 2π 0 Φ(r, θ, 0)rdrdθ: Hermite (top), Padé (middle). Bottom: instability rates as a function of µ; comparison between solution of dispersion relation (13) (using J0( √ 2µ) or its Padé approximation) and numerical results.

Table 1 .

 1 Time (in sec) of Hermite precompute and Padé as a function of ρ with different orders for Hermite interpolation. Parameters : rmin = 0.1, rmax = 0.9, Nr = N θ = 512, N = 1024, 100 iterations of the gyroaverage.

	ρ	Hermite (4) (6) (10) (18) splines Padé
	0		6	7	9	14	8	0.6
	0.001		10	11	13	18	12	2
	0.01		40	41	43	48	53	2
	0.1		446	442 448	453	594	2
		ρ	Hermite basic Hermite precompute
		0		301		0.6
		0.001		316		0.8
		0.01		312		1
		0.1		304		5

Table 2 .

 2 Time

(in sec) as a function of ρ for Hermite without and with precomputing. Parameters : rmin = 0.1, rmax = 0.9, Nr = N θ = 128, N = 1024, order of interpolation : 4, 100 iterations of the gyroaverage.

Table 5 .

 5 Comparison Padé with Hermite (m = 5).

	ρ	Padé	N = 4 N = 8 N = 16 N = 1024
	0	10 -17	10 -17	10 -17	10 -17	10 -17
	0.001	1.10 -8	5.10 -6 5.10 -6	5.10 -6	5.10 -6
	0.01	2.10 -6	2.10 -6 1.10 -6	8.10 -7	8.10 -7
	0.1	1.10 -3	1.10 -3 3.10 -5	2.10 -7	1.10 -7

Table 6 .

 6 Comparison Padé with splines (m = 5).

	ρ	Padé	N = 4 N = 8 N = 16 N = 1024
	0	10 -18	10 -17	10 -17	10 -17	10 -17
	0.001	1.10 -8	5.10 -6 4.10 -6	4.10 -6	4.10 -6
	0.01	6.10 -6	1.10 -6 1.10 -6	7.10 -7	7.10 -7
	0.1	9.10 -3	3.10 -3 1.10 -5	9.10 -8	6.10 -8

Table 7 .

 7 Comparison Padé with Hermite (m = 20).

	ρ	Padé	N = 4 N = 8 N = 16 N = 1024
	0	10 -18	10 -17	10 -17	10 -17	10 -16
	0.001	1.10 -8	4.10 -6 4.10 -6	4.10 -6	4.10 -6
	0.01	6.10 -6	1.10 -6 1.10 -6	8.10 -7	8.10 -7
	0.1	9.10 -3	3.10 -3 1.10 -5	1.10 -7	1.10 -7

Table 8 .

 8 Comparison Padé with splines (m = 20).

	ρ	Hermite(4)	(6)	(10)	(18)
	0	10 -17	10 -17	10 -17	10 -17
	0.001	3.10 -6	3.10 -6 3.10 -6 3.10 -6
	0.01	5.10 -7	6.10 -7 6.10 -7 6.10 -7
	0.1	5.10 -8	8.10 -8 1.10 -7 1.10 -7

Table 9 .

 9 Comparison

Hermite with different interpolation order. Parameters : N = 1024, m = 1.

Table 10 .

 10 Time (in min.) on HPC, mesocenter of the University of Strasbourg, with 32 processors (4 nodes). Nr × N θ × Nz × Nv = 64 × 64 × 32 × 64, ∆t = 5, 1600 iterations. For Hermite, with precomputation : N = 1024.

	µ	0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9
	Herm.	15	16	15	16	16	16	16	16	15
	Padé	15	15	14	15	17	15	16	16	15
			µ		Hermite	Padé
			0.5		12524	13326
			1		12556	13454

Table 11 .

 11 Final time (in sec.) reached for a 24 hours simulation with Nr × N θ × Nz × Nv = 128 × 256 × 128 × 128, ∆t = 2. On Helios Computational Simulation Centre, International Fusion Energy Research Centre of the ITER Broader Approac. Supercomputer with 128 processors (8 nodes; each node having 16 threads). For Hermite, with precomputation : N = 1024.

Table 12 .

 12 have validated the gyroaverage computation on polar geometry. Comparisons are made with classical Padé Herm. 1.145e -4 5.89e -4 8.175e -4 8.276e -4 5.4e -4 Linear mode growth rates for the Cyclone DIID-D base case.

	(m, n)	(4, 3)	(10, 7)	(14, 10)	(17, 12)	(21, 15)
	Padé	9.6e -5	5.9e -4	8.34e -4	8.63e -4	6.29e -4