Marcin Pietrzak
email: marcin.pietrzak@ncbj.gov.pl

Heidi Nettelbeck

Yann Perrot

Carmen Villagrasa

Aleksandr Bancer

Marion Bug

Sebastien Incerti

Intercomparison of nanodosimetric distributions in nitrogen simulated with Geant4 and PTra track structure codes

Keywords: nanodosimetry, Geant4-DNA, particle track structure, ionising radiation, Monte Carlo simulation, electron cross sections

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

Introduction

Fundamental knowledge of the mechanisms driving ionising radiation interactions with biological targets is the key to answer many questions that persist in the field of radiation protection, particularly for medical applications. Being able to adapt current approaches for dosimetry to novel radiation therapy treatments, such as hadrontherapy, BNCT (boron neutron capture therapy) or the use of metal nanoparticles as radiosensitisers [START_REF] Rabus | Intercomparison of Monte Carlo calculated dose enhancement ratios for gold nanoparticles irradiated by X-rays: Assessing the uncertainty and correct methodology for extended beams[END_REF], as well as addressing questions related to exposure at very low doses (where only stochastic effects are expected), are just some of the research challenges that would benefit from such knowledge.

Whilst achieving a complete picture of how biological effects are induced by ionising radiation might be seen as a long-term goal, significant experimental and computational research conducted during the last decades have advanced our understanding of the different processes leading to radiation-induced effects. In this context, Monte Carlo (MC) track structure codes [START_REF] Nikjoo | Track-structure codes in radiation research[END_REF][START_REF] Tang | Evaluation, from nanodosimetric modeling, of the influence of chromatin compaction on early radiation-induced effects and extension to late effects (DNA damage repair and cell death)[END_REF] are a powerful means to simulate, on an event-by-event basis, the detailed pattern of radiation interactions (on the nanometer scale) with individual intracellular molecules, which are at the origin of these effects. For that purpose, track structure codes include the interaction cross sections needed to transport electrons down to energies of a few eV. While this method may be controversial for charged particles with kinetic energies lower than 1 keV [START_REF] Thomson | On the Monte Carlo simulation of electron transport in the sub-1 keV energy range[END_REF], it has been argued [START_REF] Liljequist | On the validity of trajectory methods for calculating the transport of very low energy ( < 1 keV) electrons in liquids and amorphous media[END_REF] that statistical results obtained with MC track structure codes do provide a good approximation of multiple quantum scattering down to electron energies in the order of 10 eV.

Generally, track structure codes aim to include cross sections for calculating the probability of interactions that lead to energy deposition in condensed matter. In the case of low-energy charged particles, the interaction probability not only depends on the atomic composition of the target, but also on its aggregation state (i.e. gas, liquid or solid). Hence, there is an additional complexity to obtain interaction cross sections experimentally. For liquid water, which is considered one of the most abundant components of biological materials, few experimental cross sections exist [START_REF] Heller | Collective oscillation in liquid water[END_REF][START_REF] Hayashi | The complete optical spectrum of liquid water measured by inelastic x-ray scattering[END_REF]. Nevertheless, most track structure codes for radiobiological modelling use liquid water as a surrogate for all biological target materials, such as the cell nucleus or DNA [START_REF] Emfietzoglou | Accurate electron inelastic cross sections and stopping powers for liquid water over the 0.1-10 keV range based on an improved dielectric description of the Bethe surface[END_REF][START_REF] Incerti | Comparison of GEANT4 very low energy cross section models with experimental data in water[END_REF].

Existing experimental approaches, however, cannot provide any measurable quantities describing particle track structure on the nanometric scale in liquid water. Even the use of low-pressure water vapour in experimental setups is challenging. Thus, the most common experimental approach is to use more inert gases like nitrogen or propane. For example, microdosimeters often use tissue-equivalent gases comprising a mixture of propane, nitrogen and carbon dioxide, while nanodosimeters are usually operated with nitrogen or propane [START_REF] Bantsar | State of the art of instrumentation in experimental nanodosimetry[END_REF].

MC tools capable of simulating such experiments are not only scarce, but also challenging to develop and implement. In the case of nitrogen, one of the most tested MC track structure codes is the PTra code developed at PTB [START_REF] Grosswendt | Transport of low energy electrons in nitrogen and air[END_REF][START_REF] Grosswendt | The track structure of α-particles from the point of view of ionization-cluster formation in "nanometric" volumes of nitrogen[END_REF][START_REF] Bug | Ionization cross section data of nitrogen, methane, and propane for light ions and electrons and their suitability for use in track structure simulations[END_REF]. The code has been successfully applied to simulate experiments, in particular, those performed using three different nanodosimeters using nitrogen or propane as a target gas [START_REF] Grosswendt | The track structure of α-particles from the point of view of ionization-cluster formation in "nanometric" volumes of nitrogen[END_REF][START_REF] Grosswendt | Experimental equivalent cluster-size distributions in nanometric volumes of liquid water[END_REF][START_REF] Bashkirov | Ion-counting nanodosemeter with particle tracking capabilities[END_REF][START_REF] Bug | Proton-impact ionisation cross sections for nanodosimetric track structure simulations[END_REF][START_REF] Bug | Nanodosimetric characterization of ion beams[END_REF][START_REF] Hilgers | Comparison of measuredand monte carlo simulated track structure parameters in nanometric volumes[END_REF][START_REF] Conte | Nanodosimetric track structure in homogeneous extended beams[END_REF].

Written in FORTRAN, the PTra code provides a rudimentary interface regarding the target and beam geometry. Only simple geometries can be easily implemented and anything more complicated demands careful changes in many of the subroutines, which is prone to human error. One solution is to implement appropriate cross section models for these gases in a more contemporary and open source MC toolkit such as Geant4 [START_REF] Agostinelli | Geant4 -a simulation toolkit[END_REF], which provides a more versatile interface for defining beam and detector geometries. Furthermore, the widespread use of Geant4 compared to the PTra code would enable the simulation tools based on such cross-section data to be available to more scientific users.

The Geant4-DNA collaboration therefore initiated a new activity to implement interaction cross section data of materials other than liquid water for track structure simulations. This began with four DNA-substitute materials tetrahydrofuran (THF), trimethyl phosphate, (TMP), purine (Pu) and pyrimidine (Py), which were already implemented in PTra [START_REF] Bug | Ionization cross section data of nitrogen, methane, and propane for light ions and electrons and their suitability for use in track structure simulations[END_REF], and has been furthered in this work with the inclusion of nitrogen. In future releases, it is planned to also include cross section data for propane. For this purpose, specific model classes were designed and implemented in the Geant4-DNA extension of the Geant4 code (i.e. in the Geant4DNAPTB model classes) [START_REF] Incerti | Comparison of GEANT4 very low energy cross section models with experimental data in water[END_REF][START_REF] Bernal | Track structure modeling in liquid water: A review of the Geant4-DNA very low energy extension of the Geant4 Monte Carlo simulation toolkit[END_REF][START_REF] Incerti | The Geant4-DNA project[END_REF][START_REF] Incerti | Geant4-DNA example applications for track structure simulations in liquid water: A report from the Geant4-DNA Project[END_REF]. Currently, these publicly available models can be used to perform track structure simulations in the four DNAsubstitute materials mentioned above, as demonstrated in the Geant4 "ICSD" user example [START_REF] Incerti | Comparison of GEANT4 very low energy cross section models with experimental data in water[END_REF][START_REF] Bug | An electron-impact cross section data set (10 eV-1 keV) of DNA constituents based on consistent experimental data: A requisite for Monte Carlo simulations[END_REF].

In this work, we validate the inclusion of appropriate cross section models for nitrogen in Geant4-DNA prior to its public release. Using these newly implemented data, the ability of the modified Geant4-DNA toolkit to simulate certain physical track-structure quantities is benchmarked against the PTra code in which the cross sections were first implemented. Physical track structure quantities relevant to nanodosimetry, such as the ionisation cluster size distribution (ICSD), were simulated with both Geant4-DNA and PTra using identical cross section data. In the course of comparing the transport methods and models adopted in both codes, some minor bugs were corrected in the PTra code. A validity check of the simulation results (obtained with the new cross section data) was carried out by benchmarking the simulated electronic stopping power and range in nitrogen with corresponding ICRU and NIST data [START_REF]ESTAR database[END_REF] for energies greater than 10 keV. For energies less than 10 keV, the comparison was made against a model presented by Gümüş [START_REF] Gümüş | Simple stopping power formula for low and intermediate energy electrons[END_REF]. In this low energy range, simulated ICSDs were also compared to those obtained experimentally (in nitrogen) by Bancer (aka Bantsar) et al. using the Jet Counter nanodosimeter at NCBJ [START_REF] Bantsar | Single track nanodosimetry of low energy electrons[END_REF]. The good agreement between Geant4-DNA and PTra codes for simulated quantities such as ICSDs as well as experimental results and published data on stopping power and range values satisfies the requirement for their release in the Geant4 code. Furthermore, it facilitates the next stage of this work, which is the implementation in the same model classes of interaction cross section data for propane.

In a broader context, this work is aligned with the Geant4 developers' initiative to provide a well-documented and tested open source MC tool for a variety of medical physics applications, from radiation quality studies through radiobiology to patient treatment with radiotherapy [START_REF] Incerti | Advances in Geant4 applications in medicine[END_REF].

Materials and methods

Physical processes and models for nitrogen electron interactions

With the aim of validating the implementation of nitrogen cross sections in Geant4-DNA classes, this study compares the performance of the MC track structure codes Geant4 (version 10.06.p01) and PTra (versions V10/08 and V22/04) using identical interaction cross section data for nitrogen. The key features of the cross section models are extensively described by Grosswendt and Waibel [START_REF] Grosswendt | Transport of low energy electrons in nitrogen and air[END_REF] with further amendments in [START_REF] Grosswendt | The track structure of α-particles from the point of view of ionization-cluster formation in "nanometric" volumes of nitrogen[END_REF][START_REF] Bug | Ionization cross section data of nitrogen, methane, and propane for light ions and electrons and their suitability for use in track structure simulations[END_REF]. An overview of these models is given in Table 1 and a brief description is provided in the subsequent sections.

In the case of Geant4-DNA, the cross sections were implemented in the form of lookup tables so that they can be easily replaced, if needed, in a future revised version of the code. For example, it is planned to introduce an alternative elastic scattering model for nitrogen based on recent Table 1: Summary of the models used to calculate interaction cross sections in nitrogen, which were obtained from PTra and then implemented in Geant4-DNA.

Interaction type

Cross section type Models

Impact ionisation

Total (including partial cross sections for subshells)

Binary-Encounter-Bethe (BEB) model [START_REF] Hwang | New model for electronimpact ionization cross sections of molecules[END_REF], except for K-shell ionisation which uses one proposed by Casnati et al. [START_REF] Casnati | An empirical approach to K-shell ionisation cross section by electrons[END_REF].

Differential (production energy of secondaries)

Breit-Wigner formula [START_REF] Green | Ionization cross sections and secondary electron distributions[END_REF].

Doubly-differential (production angle for given production energy)

Berger's kinematic model with an extension for low energies [START_REF] Grosswendt | Transport of low energy electrons in nitrogen and air[END_REF] Elastic scattering Total Fitted experimental data [START_REF] Grosswendt | Transport of low energy electrons in nitrogen and air[END_REF].

Differential (scattering angle)

Modified Rutherford model with atomic screening corrections [START_REF] Grosswendt | Transport of low energy electrons in nitrogen and air[END_REF].

Electronic excitation

Total (including 29 partial cross sections)

Formulas and cross section parameters based on Porter et al. [START_REF] Porter | Efficiencies for production of atomic nitrogen and oxygen by relativistic proton impact in air[END_REF].

experimental data for electrons of energies ranging from 30 to 1000 eV once the experimental data is published. In all simulations and calculations, an energy cut of 15.58 eV (i.e. the ionisation threshold of nitrogen) was used for electron transport.

While interaction cross section models presented in Table 1 are implemented in the PTra code, they are often not used directly. Rather, lookup tables are generated in the initialisation phase which are then used during the simulation of track structure. The Geant4-DNA classes were therefore designed to use pre-generated lookup tables with an identical format as those used in the internal tables of PTra. A different approach, however, was used to calculate the production energy of secondaries (in the ionisation model) and the calculation of the scattering angle (in the elastic scattering model). In both cases, PTra uses formulas that directly describe particular models, while Geant4-DNA uses pre-generated lookup tables based on these formulas. This difference in code design allows the necessary generalisation of Geant4-DNA to simulate a variety of materials using different sources of cross section data, often providing alternatives for the same material. This dissimilarity between PTra and Geant4-DNA can lead to residual differences in the simulation results obtained with these codes.

Impact ionisation model

The Binary-Encounter-Bethe (BEB) model [START_REF] Kim | Binary-encounter-dipole model for electron-impact ionization[END_REF] with the molecular orbital data from Hwang et al. [START_REF] Hwang | New model for electronimpact ionization cross sections of molecules[END_REF] is used for describing the electron impact ionisation processes in nitrogen. The model provides partial ionisation cross sections for the four subshells with binding energies of 15.58 eV, 17.07 eV, 21.00 eV and 41.72 eV. With a binding energy of 409.5 eV, the K-shell is treated separately according to the empirical formula proposed by Casnati et al. [START_REF] Casnati | An empirical approach to K-shell ionisation cross section by electrons[END_REF]. The energy of a secondary electron emitted after impact ioni-sation is determined by the Breit-Wigner formula. [START_REF] Green | Ionization cross sections and secondary electron distributions[END_REF]. For any given energy, the total ionisation cross section should be the sum of the cross sections for all five subshells. However, if the auto-ionisation process is activated, as is the default behaviour in both codes, then the total impact ionisation cross section is reduced by an amount that is then added to the total excitation cross section. This auto-ionisation process corresponds to an emission of a secondary electron after the excitation of one or more nitrogen shells. Therefore, the total emission of secondary electrons produced either by ionisation or auto-ionisation must be preserved in accordance with the total ionisation yield. Compared to when it is inactivated, activation of the auto-ionisation option can lead to small variations in energy and angular distribution of secondary electrons.

Excitation model

In this model, there is no scattering (change in momentum direction) of incident electrons. The treatment of electronic excitation processes is based on the formulas and cross section parameters of Porter et al. [START_REF] Porter | Efficiencies for production of atomic nitrogen and oxygen by relativistic proton impact in air[END_REF]. To obtain a better agreement with the experimental total scattering cross sections, a few corrections and extensions have been applied [START_REF] Grosswendt | The track structure of α-particles from the point of view of ionization-cluster formation in "nanometric" volumes of nitrogen[END_REF].

The auto-ionisation process is part of the excitation model in both MC codes, where a probability of 50% autoionisation is assumed if the excitation energy of a Rydberg state is greater than the 15.58 eV ionisation threshold for nitrogen. The secondary electrons produced by autoionisation are emitted isotropically and their kinetic energy is calculated from the given excitation energy minus the ionisation threshold.

Elastic scattering model

The treatment of elastic electron scattering is based on integrated cross sections σ el (T) obtained by experiments and on Rutherford's differential cross sections (dσ/dΩ) el with respect to the solid angle, which were modified to take into account atomic screening effects [START_REF] Grosswendt | Transport of low energy electrons in nitrogen and air[END_REF].

Auger cascades and auto-ionisation

As mentioned above, PTra has the capability of simulating auto-ionisation following an excitation as well as Auger cascades after a K-shell ionisation of the target molecule. While it is also possible to simulate Auger cascades in Geant4-DNA, both codes have different implementations of the process, which may lead to differences in the simulated results. Since the Auger cascade process is independent of the cross section implementation to be validated in this work, the Auger cascade model was inactivated in both PTra and Geant4-DNA. The autoionisation process, on the other hand, was incorporated in the G4DNAPTBExcitationModel class in a similar manner to that in PTra, and thus is included in the comparison of simulated results obtained with both codes.

Modifications to the PTra code

Prior to implementing the physical models in Geant4-DNA, the PTra code was updated to version V22/04. This included some minor changes in an older version (V18/08) of PTra that were necessary to perform simulations for the intercomparison for specific target and beam geometries. It was during this process that a few minor bugs in PTra were identified and subsequently corrected, namely some outdated physics constants (minor differences compared to contemporary values), miscalculation of a few percent for the total ionisation cross section due to omission of the K-shell contribution in the summation (up to 100 keV of electron energy, the K-shell contribution is less than 5%), a mistake in the formula for production angles of secondary electrons with energies below 200 eV. (Owing to the limited range and high elastic-scattering cross section of these electrons, the exact angular spectrum during production is of little importance.) The overall impact of these minor bugs is only apparent during sensitive tests and, thus, in most typical cases has little effect on the simulation results. Consequently, simulated results were compared for the following versions of the codes:

1. PTra V22/04 (updated version of PTra V18/08); 2. Geant4-DNA (10.06.p01): including all cross section data derived from models implemented in PTra V22/04; 3. PTra V18/08: initial version of the PTra code used in previous studies. Simulations were performed with this version to determine the impact of the aforementioned code upgrade on the results.

Simulated geometries and energies

To compare and validate the results of both codes, different combinations of simulation geometries and energies were investigated, both with nitrogen target volumes of density equal to 1 g/cm 3 (i.e. liquid water density). The world volumes in which these targets were placed were also filled with nitrogen of density equal to 1 g/cm 3 . This density was chosen to enable the simulations to be performed in target volumes of sizes equivalent to those of biological targets, such as DNA constituents. In the first geometry, a cylindrical target 2.3 nm in diameter and 3.4 nm in height (representing the typical dimensions of a DNA segment comprising 10 base pairs) was placed in the centre of a 10×10×10 nm 3 cubic world volume. The target was irradiated with a pencil beam of monoenergetic electrons originating on the side of the cylinder at half of its height. This geometry is identical to the one used in the Geant4 "ICSD" example, which is already available for simulations with DNA-substitute materials (i.e. THF, TMP, Pu and Py). The objective of this work is to extend this ICSD example (in a future Geant4 release) to include nitrogen as a potential material.

In the second geometry, a spherical target of 8 nm in diameter (representing the volume of a nucleosome) was placed at the centre of a 20×20×20 nm 3 cubic world. Both the target and world were filled with nitrogen of density equal to 1 g/cm 3 . The target was irradiated by an isotropic source of monoenergetic electrons originating at the centre of the sphere. The use of these two targets allows the investigation of different ICSDs in terms of mean values and variances to cover all possible discrepancies that may arise in the results obtained with the two codes. For both geometrical setups, at least 10 7 primary electrons were simulated for each initial kinetic energy: 20 eV, 30 eV, 50 eV, 70 eV, 100 eV, 150 eV, 200 eV, 300 eV, 400 eV, 500 eV, 600 eV, 800 eV, 1 keV, 2 keV, 3 keV, 5 keV, 10 keV, 20 keV, 50 keV, 100 keV, 200 keV, 500 keV and 1 MeV, which led to results with statistical uncertainties of less than 1%.

Comparison with experimental data

Simulated ICSDs were also compared with those obtained experimentally in order to benchmark the cross section data implemented in both Geant4-DNA and PTra codes. Measured ICSDs for low-energy monoenergetic electron beams traversing a "nanometric" nitrogen target are reported by Bancer et al. [START_REF] Bantsar | Single track nanodosimetry of low energy electrons[END_REF][START_REF] Bantsar | Ionization Cluster Size Distributions Created by Low Energy Electrons and Alpha Particles in Nanometric Track Segment in Gases[END_REF]. The study includes a detailed description of the Jet Counter nanodosimeter and experimental setup used to carry out the measurements.

The most important aspects of the experimental setup can be summarised as follows. The target in the Jet Counter is created dynamically by a pulse-operated piezoelectric valve that injects a jet of nitrogen gas to a interaction chamber (IC). The IC volume is an open (bottomless) cylinder, 10 mm in diameter and 10 mm in height. Following passage of the projectile, ions created within the IC escape with the gas jet to a much larger vacuum chamber. While neutral molecules dissipates within the chamber, an electric field guides all ions to a single ion detector which then counts them one-by-one. Thus, the entire IC region is a sensitive volume of the detector, in which all ions created inside are measured with a known detection efficiency.

For each event, the target was irradiated with a single pulse from the electron gun. The average number of electrons per beam pulse was approximately one, but varied slightly from experiment to experiment (see Table 2). Despite the real number of electrons for any given event being random, it was verified to obey a Poisson process. This of course means that some events were blank shots (no electron was emitted from the electron gun). In order to compare these measurements with simulations, a Poisson probability distribution for the number of electrons per each event was taken into account.

The nitrogen gas density within the nanodosimeter was tuned to a value of 0.34 µg/cm 3 . Applying the scaling procedure described by Grosswendt [START_REF] Grosswendt | Recent advances of nanodosimetry[END_REF], this cylindrical nitrogen volume could be scaled to a cylindrical water target of unit density (roughly 2.3 nm in both diameter and height), which is equivalent in size to a short segment of DNA. In both PTra and Geant4-DNA simulations, however, the target's geometry and density were the same as those used in the experiment.

Reference quantities

For comparison purposes, ICSD was chosen as the basic physical reference quantity as it is usually regarded as one of the fundamental descriptors of the particle track structure. It is for this reason that many papers in the field of nanodosimetry are focused on this quantity and/or its derivatives. Sensitive to different aspects of the ICSD shape, the following three parameters of the ICSD were therefore chosen for direct comparison: the mean ionisation cluster size M 1 (i.e. first moment of the ICSD), which reflects the total ionisation yield; the variance σ 2 of the cluster size, which is most sensitive to the distribution's tail; and the probability P 1 to obtain a cluster size of one (i.e. only a single ionisation in the target), which for the geometries of interest was sensitive to both the production and behaviour of secondary electrons. The latter can lead to additional ionisations, thereby significantly lowering P 1 . Of these nanodosimetric quantities, M 1 is of most interest as the other two are sensitive to subtle differences in the shape of the ICSD.

The statistical uncertainty is always given as a standard deviation of a quantity. For M 1 and σ 2 , the uncertainty was calculated from the statistical properties of ICSD from which the quantity is derived. For P k (and P 1 in particular), the standard deviation is the square root of the number of counts of a given cluster size divided by the total number of events.

Results and discussion

Figures 1 and2 show the energy dependence of each of the three above-mentioned quantities for PTra V22/04 and Geant4-DNA. These results should reflect the accuracy and suitability of the implemented physical models in both codes. For both of the investigated geometries, there is a good agreement between the simulated results for all three quantities (less than ∼1% and ∼3% difference for the first and second geometries, respectively). This not only verifies the correct implementation of nitrogen cross sections in Geant4-DNA, but also suggests that the transport in both track structure codes is equivalent. Any residual (minor) discrepancies can be attributed to differences in the interpolation method of cross section values used in each respective MC code.

In the current study, the cross section models obtained from PTra and implemented in Geant4-DNA were also validated by calculating the stopping power and range of electrons in nitrogen using equivalent implementation in both codes (i.e. PTra V22/04 and Geant4-DNA). The results of this comparison are shown in Figure 3 [START_REF]ESTAR database[END_REF] and Gümüş model [START_REF] Gümüş | Simple stopping power formula for low and intermediate energy electrons[END_REF] calculated using CSDA. The cutoff energy (i.e. when electron transport is ceased) in both simulations and calculations was assumed to be 15.58 eV (nitrogen ionisation threshold). The relative difference between Geant4-DNA simulation results and the other data is shown in the bottom panel of the figure. Statistical uncertainties are smaller than the symbol size in both panels.

data from the NIST ESTAR database [START_REF]ESTAR database[END_REF] (based on ICRU report 37 [START_REF]Stopping Powers for Electrons and Positrons[END_REF]) and Gümüş model [START_REF] Gümüş | Simple stopping power formula for low and intermediate energy electrons[END_REF]. The overestimation of the stopping power and corresponding underestimation of range is most likely related to the overestimation of impact ionisation cross sections in the BEB model (in relation to the values recommended by Itikawa [START_REF] Itikawa | Cross Sections for Electron Collisions with Nitrogen Molecules[END_REF], which was discussed by Bug et al. [START_REF] Bug | Ionization cross section data of nitrogen, methane, and propane for light ions and electrons and their suitability for use in track structure simulations[END_REF]). The substantial underestimation of stopping power beyond 200 keV has little impact 0.0 0.5 in most applications involving radiobiology or nano-and microdosimetry. For example, in the case of a 230 MeV proton beam, which is typically the maximum initial energy of clinical beams, less than 1% of secondary electrons are produced with energies greater than 100 keV with an average value of just a few keV. Finally, to serve as a reference for previous studies performed with earlier versions of the PTra code, we show the influence of the modifications made to PTra in this work. The overall impact of all bug fixes is shown in Figures 5 and6, where the most significant differences are present for high-energy electrons. This suggests that the bugs present in older versions of the PTra code have little impact on any previous results for heavy charged particles (protons and alpha particles) since the upper energy limit of secondary electrons produced by primary particles in these studies was typically less than 10 keV.

Experimental validation

Figure 7 shows a comparison of the measured [START_REF] Bantsar | Single track nanodosimetry of low energy electrons[END_REF] and simulated ICSDs in nitrogen. In the interest of readability, only the simulation results obtained with Geant4-DNA are shown. (All data points obtained with both PTra and Geant4-DNA codes were found to overlap.) For the same reason, only four energies are shown. Table 2 compares the values of M 1 , P 1 and σ 2 for all available experimental data.In the presented simulation results, the ionisation yield for 1 and 2 keV electron beams may be underestimated due to the absence of Auger electron cascade contributions. However, based on calculations with PTra, the impact on M 1 is in the order of 1% . 7: Comparison of the simulated (empty symbols) and measured (full symbols) ionisation cluster size distributions for different monoenergetic electron beams. The simulated cylindrical target geometry and density were identical to those adopted in the experiment (i.e. diameter and height of approximately 2.3 nm and density of 0.34 µg/cm 3 . Statistical uncertainty is marked only for experimental values as for simulated it is smaller then the symbol size. Table 2: Experimental (exp.) [START_REF] Bantsar | Single track nanodosimetry of low energy electrons[END_REF] and simulation (sim.) results for the ICSD parameters M 1 , P 1 and σ 2 obtained from monoenergetic electrons of energy E traversing a nitrogen target of 0.34 µg/cm 2 thickness (i.e. equivalent geometry to that of the Jet Counter nanodosimeter). The beam intensity I is defined as the average number of primary electrons per event. Statistical uncertainties of the last digits are given in parentheses.

Conclusion and further work

The results in this work indicate a successful implementation of interaction cross section models for nitrogen, namely elastic scattering, impact ionizations and electronic excitations, in the Geant4 toolkit. These cross sections are identical to those originally implemented in the PTra code, which has been used extensively over the last few decades to simulate nanodosimetric experiments. This was shown in the comparison of three physical quantities (M 1 , P 1 and σ 2 ) describing the ICSD. A comparison of the energy dependence of these quantities with two different target geometries using Geant4-DNA and PTra V22/04 (in Figures 1 and2) showed a good agreement between the codes. It can therefore be assumed that the transport in both MC track structure codes is similar, thus verifying a correct implementation of these cross section models for nitrogen in Geant4-DNA. Furthermore, the good agreement between the simulated stopping powers and ranges as well as ICSDs for monoenergetic electrons over the entire energy range also validates the use of these models in both codes.

Hence, these nitrogen cross sections will be included in the next release of the Geant4 code for public use. Further work is planned to extend this data set to include the implementation of appropriate cross section models for propane, which is another tissue-equivalent gas often chosen as the operating gas in micro-and nanodosimeters. This inclusion of nitrogen and propane cross section data in Geant4-DNA would not only improve the accuracy of radiation transport simulations in these media, but also increase the code's versatility for applications particularly within the micro-and nanodosimetry community.

Figure 1 :Figure 2 :

 12 Figure 1: Energy dependence of the three quantities M 1 , P 1 and σ 2 obtained with the cylindrical geometry for equivalent simulations performed using Geant4-DNA and PTra codes. The relative difference is shown in the bottom panel of the figure. Statistical uncertainties are smaller than the symbol size in both panels.

Figure 3 :

 3 Figure3: Energy dependence of the electronic stopping power in nitrogen for equivalent simulations performed using Geant4-DNA and PTra codes compared with data from NIST ESTAR database[START_REF]ESTAR database[END_REF] and Gümüş model[START_REF] Gümüş | Simple stopping power formula for low and intermediate energy electrons[END_REF]. The relative difference between Geant4-DNA simulation results and the other data is shown in the bottom panel of the figure. Statistical uncertainties are smaller than the symbol size in both panels.

Figure 4 :

 4 Figure4: Energy dependence of the electron range in nitrogen for equivalent simulations with Geant4-DNA and PTra codes compared with data from NIST ESTAR database[START_REF]ESTAR database[END_REF] and Gümüş model[START_REF] Gümüş | Simple stopping power formula for low and intermediate energy electrons[END_REF] calculated using CSDA. The cutoff energy (i.e. when electron transport is ceased) in both simulations and calculations was assumed to be 15.58 eV (nitrogen ionisation threshold). The relative difference between Geant4-DNA simulation results and the other data is shown in the bottom panel of the figure. Statistical uncertainties are smaller than the symbol size in both panels.

Figure 5 :Figure 6 :

 56 Figure 5: Energy dependence of three quantities M 1 , P 1 and σ 2 obtained in a cylindrical geometry from simulations with PTra with and without updates to the code ((version V22/04 and V18/08 respectively). The relative difference is shown in the bottom panel of the figure. Statistical uncertainties are smaller than the symbol size in both panels.

Figure

  Figure7: Comparison of the simulated (empty symbols) and measured (full symbols) ionisation cluster size distributions for different monoenergetic electron beams. The simulated cylindrical target geometry and density were identical to those adopted in the experiment (i.e. diameter and height of approximately 2.3 nm and density of 0.34 µg/cm 3 . Statistical uncertainty is marked only for experimental values as for simulated it is smaller then the symbol size.

Acknowledgements

The resources of the Świerk Computing Centre (CIŚ) financed by the EU and the Polish Ministry of Science and Higher Education (project No. POIG.02.03.00-00-013/09) were used to carry out this work.


Yu Yeh

Nelson Cisneros

Yongxin Wu

Kanty Rabenorosoa

Yann Le Gorrec

Modeling and Position Control of the HASEL Actuator via port-Hamiltonian Approach

Keywords: Soft robot, soft actuator, HASEL actuator, port-Hamiltonian approach, IDA-PBC method

This paper deals with the modeling and control problem of a Hydraulically Amplified Self-healing Electrostatic (HASEL) actuator based on the port-Hamiltonian framework. A nonlinear spring-damper system is used to approximate the mechanical deformation of the actuator due to the motion of the fluid while a nonlinear capacitance is used to approximate the electric behavior of the system. The actuator position control strategy is investigated based on the Interconnection Damping Assignment-Passivity Based Control (IDA-PBC) method, with further Integral Actions (IA) added to cope with load uncertainties. The proposed model and control laws are validated on an experimental benchmark. The experimental tests demonstrate that the proposed model is accurate up to 94% of fitness. The controllers allow assigning the actuator position with ramp and sinusoidal references with the relative error less than 5%. At last, the robustness of the proposed IDA-PBC controller with IA has been shown with the experimental result for the unknown load disturbance rejection.

I. INTRODUCTION

Robots based on soft actuators have drawn more and more attention of the researchers during the last decade due to their safety for humans and suitability to unknown environments. Recent studies have shown that robots based on soft actuators perform competitively with existing solid robots in terms of force and strain [START_REF] Rus | Design, fabrication and control of soft robots[END_REF]. Soft pneumatic actuators and dielectric elastomer actuators are two popular kind of soft actuators. In the recent years, the Hydraulically Amplified Self-Healing Electrostatic (HASEL) actuator has provided a new road-map for the design of soft actuators [START_REF] Acome | Hydraulically amplified self-healing electrostatic actuators with muscle-like performance[END_REF]. By combining the design concept of dielectric elastomer actuators and pneumatic soft actuators, more desirable features of HASEL actuators are indicated [START_REF] Rothemund | Hasel artificial muscles for a new generation of lifelike robots-recent progress and future opportunities[END_REF] including self-sensing and self-healing capabilities and great potential for many applications, such as: HASEL artificial muscles made from elastomers, elastomeric donuts, quadrant donuts and curling HASEL actuators [START_REF] Rothemund | Dynamics of electrohydraulic soft actuators[END_REF], soft grippers [START_REF] Acome | Hydraulically amplified self-healing electrostatic actuators with muscle-like performance[END_REF] and soft-actuated joints created based on the hydraulic mechanism used in spider legs [START_REF] Kellaris | Spider-inspired electrohydraulic actuators for fast, soft-actuated joints[END_REF]. One of the advantages of the HASEL actuator is that it is easy to manufacture, low cost, and open to many potential design based on its working principle [START_REF] Kellaris | Peano-hasel actuators: Muscle-mimetic, electrohydraulic transducers that linearly contract on activation[END_REF]. The authors of [START_REF] Mitchell | An easy-to-implement toolkit to create versatile and high-performance hasel actuators for untethered soft robots[END_REF] present an easy to implement toolkit to design and fabricate multiple HASEL This work has been supported by the EIPHI Graduate School (contract ANR-17-EURE-0002) and the ANR IMPACTS project (contract: ANR-21-CE48-0018).

Y. Yeh, N. Cisneros, Y. Wu, K. Rabenorosoa, Y. Le Gorrec are with the FEMTO-ST institute, UMR CNRS 6174, département AS2M, Université Bourgogne Franche-Comté , ENSMM, F-25000 Besanc ¸on, France (e-mail: yongxin.wu@femto-st.fr).

actuators incorporating electrostatic zipping mechanisms and reducing operating voltages. Several attempts to model and control the HASEL actuator have recently been proposed. In [START_REF] Kellaris | An analytical model for the design of peano-hasel actuators with drastically improved performance[END_REF], the geometrical analysis of HASEL is investigated. In [START_REF] Johnson | Identification and control of a nonlinear soft actuator and sensor system[END_REF], an identification of a nonlinear HASEL actuator model has been presented and the controller design using the self-sensing strain is addressed with a basic control approach. In [START_REF] Liu | Modeling of planar hydraulically amplified self-healing electrostatic actuators[END_REF], the planar HASEL actuator model is investigated to describe the relationship between the output force and the applied voltage. These works reveal that the modeling and control of complex mechanical structure with fluid structure interaction is very hard and challenging. Only few studies on controloriented modeling and control design has been proposed in the literature [START_REF] Johnson | Identification and control of a nonlinear soft actuator and sensor system[END_REF].

To overcome these constraints, this paper proposes a physically based dynamic model and a control strategy for HASEL actuators using the port-Hamiltonian (PH) approach. This approach [START_REF] Maschke | Port-controlled Hamiltonian systems: Modelling origins and systemtheoretic properties[END_REF] has proven to be powerful for the modeling and control of complex multi-physical systems [START_REF] Duindam | Modeling and Control of Complex Physical Systems: The Port-Hamiltonian Approach[END_REF], [START_REF] Dòria-Cerezo | Passivity-based control of a wound-rotor synchronous motor[END_REF]. PH modeling is based on the characterization of energy exchanges between the different components of the system. This approach is particularly adapted for the modular modeling of multiphysical systems. On the other hand, the PH approach is well suited passivity based control design with clear physical interpretation, such as energy shaping and control by interconnection and damping assignment (IDA-PBC) [START_REF] Ortega | Interconnection and damping assignment passivity-based control of port-controlled Hamiltonian systems[END_REF]. The PH approach and IDA-PBC control method have been investigated for the modeling and control of a flexible beam using the Ionic Polymer Metal Composite (IPMC) actuators in [START_REF] Zhou | Port-Hamiltonian modeling and IDA-PBC control of an IPMC-actuated flexible beam[END_REF]. In this work, the flexible beam has been considered under small deformation and modeled by a linear infinite dimensional system and a very simple RC circuit model has been considered for the IPMC actuator. The overall system is then governed by a linear coupled Partial Differential Equation (PDE) -Ordinary Differential Equation (ODE) model. The control method has been proposed for the discretized model without any investigation on the robustness of the controller.

A common design of HASEL actuators consists in clustering HASEL units as shown in the right figure of Fig 1 ,  and each unit consists of pouches filled with dielectric fluid and covered with pairs of electrodes. The pouch is formed by flexible but inextensible polymer films made bonding shell. While the voltage is applied, the corresponding Maxwell force zips the paired electrodes and displaces the dielectric liquid. Because of the inextensible properties of the polymer films and dielectric liquid, the displacement of liquid causes the deformation of the pouches, leading to the motion of endpoints of the actuator. The main contributions of this paper are the proposition of a PH model to describe the dynamics of the HASEL actuator and the position control strategy based on the IDA-PBC method. The geometry of the HASEL pouch is investigated and based on this analysis, a PH dynamic model is established using a non linear spring-damper system to approximate the mechanical behavior of the flowing liquid in the pouches. Different from to the work presented in [START_REF] Zhou | Port-Hamiltonian modeling and IDA-PBC control of an IPMC-actuated flexible beam[END_REF], a precise model is proposed considering both the electrical and mechanical parts of the actuator, leading to a non linear finite dimensional model in this paper. A control law based on IDA-PBC method is employed to achieve the desired equilibrium position. An integral action is added to the previous controller to cope with the uncertainties in the load mass which has not been considered in [START_REF] Zhou | Port-Hamiltonian modeling and IDA-PBC control of an IPMC-actuated flexible beam[END_REF]. The proposed PH model is identified and validated by an experimental setup and the proposed control method is also validated by the same experimental setup which is described in Fig. 1.

This paper is organized as follow: Section II presents the geometry analysis and the PH modeling of the HASEL actuator. The IDA-PBC based position control and integral action are proposed in Section III. In Section IV, the experimental setup and the parameters identification are presented. The proposed control laws is experimentally validated and the results are discussed in Section V. Some final remarks and future works are given in Section VI.

II. PORT HAMILTONIAN MODELING OF HASEL ACTUATOR

In the following sections, the geometric constraints of the system are introduced. All the geometric variables can be interpreted as function of the endpoint displacement q (the geometry of HASEL actuator is shown in Fig 2b), the partial derivative of these variables are naturally derived and which is essential for the PH model.

A. Geometry description of HASEL actuator

We first introduce the HASEL unit, which is a simplified geometric model to approximate the real actuator. Unlike an ellipse geometric shape presented in [START_REF] Rothemund | Dynamics of electrohydraulic soft actuators[END_REF], we use a diamond shape to describe the unit pouch of the HASEL actuator as shown in Fig. 2a.

The HASEL unit is composed of a pouches filled with dielectric liquid, and equiped with a pair of electrodes covering the two polymer films. The parameters of the actuator include: • L p the length of pouches;

• w the width of pouches and electrodes;

• L e the length of paired electrodes;

• L v the initial length of actuator;

• A the cross-sectional area of pouches;

• t the thickness of dielectric film. When we apply the voltage to the actuator, the Maxwell force zips the pair of electrodes which leads to the deformation of actuator. We define the geometric variables due to the deformation as follows:

• l e the length of the zipped electrodes;

• l h the height of the pouches;

• l v the length of the unzipped pouches;

• l c the distance between the endpoints of the paired electrodes. In this paper, we make two assumptions for the actuator during its deformation: 1) the cross-section area A is constant because the dielectric liquid is incompressible, and 2) the length of the pouch L p is constant because the membrane is inextensible. The objective is to find how the displacement q relates to the geometry of the HASEL unit. From the deformation of the actuator, we can write the geometrical relationship to achieve the different variables. The actuator endpoint displacement is defined as:

q = L v -(l e + l v ) = Ω(l e ) (1)
Since the cross-sectional area is assigned to a diamond shape, the following geometry relationship holds:

A = l v l h 2 (L p -l e ) 2 = l 2 v + l 2 h (2)
By substituting the first equation of ( 2) to the second one, we can present the height l h and the width l v as the function of l e as following:

l h = (L p -l e ) 2 2 - (L p -l e ) 4 4 -4A 2 (3)
l v = (L p -l e ) 2 2 + (L p -l e ) 4 4 -4A 2 (4)
Then we can substitute the above equation to Eq. ( 1) and the length of the zipped electrodes l e can be presented as a function of the actuator endpoint displacement q: l e (q) = Ω -1 (q)

Finally, the geometric constraints of the actuator can written as a function of the displacement q.

The more precise geometric analysis of the actuator has been made based on the ellipse shape in [START_REF] Kellaris | An analytical model for the design of peano-hasel actuators with drastically improved performance[END_REF], [START_REF] Moretti | Analysis of dielectric fluid transducers[END_REF]. In this paper we have chosen the the diamond shape for the dynamic modeling, because it is a good compromise between the accuracy of the model and its simplicity of the geometric analysis. All the geometric variables can be expressed in the function of q with the diamond shape assumption. However, in the ellipse shape case, it is not easy to present or compute all the geometric variables in the function of q.

B. PH Formulation of HASEL

In this part, we derive a dynamic model of the HASEL actuator by using the PH framework. The main idea is the use of spring damper systems to describe the mechanical properties of the dielectric liquid in pouch. The motivation is from the observation of the physical properties during the actuator deformation: 1) the displacement of the liquid and the constraints on the structure cause the resistance to the deformation, and 2) the actuator returns to its initial state (position) when the voltage is switched off. These properties imply that the mechanical properties of the flowing dielectric liquid can be described by a spring-damper system. In this paper, we propose two spring-damper systems places in the vertical and horizontal directions of the actuator, as shown in Fig. 3. To model the dynamics of the actuator, we first introduce its Hamiltonian function (energy function) which is composed by the mechanical energy and the electric energy. The total mechanical energy is the sum of the kinetic and the potential energy. The potential energy of the spring system is written as

H s (q) = 1 2 k v (q)∆ξ v (q) 2 + 1 2 k h (q)∆ξ h (q) 2 (6)
where ∆ξ v (q), ∆ξ h (q) are the deformation of the springs in the vertical and horizontal directions, and k v (q), k h (q) are the coefficients of the nonlinear springs. The deformation of the springs are defined as

∆ξ v (q) = l v (q) -ξ v (m) ∆ξ h (q) = l h (q) -ξ h (m) (7)
where ξ v , ξ h are the original length of the springs. In this work, the original length of the springs are designed as nonlinear functions depending on the load mass m. Because the geometry constrains limit initial states of the actuator, no matter how heavy the load is applied to the actuator, the initial states should remain the same. From a physical point of view, it can be explained by the fact that the uncompressible dielectric liquid generates an equivalent inner pressure to resist to the external force from the load and maintains the actuator in the initial geometry.

The electric energy of the capacitance between the paired electrodes is

H c (q, Q) = 1 2 Q 2 C(q) ( 8
)
where Q is the charge of the capacitor and C is the capacitance value. While C depends on q, its computation can be split into two parts: the zipped region and unzipped region. The capacitance in the zipped region can be described by

C z (q) = r 0 w 2t l e (q) (9)
For the unzipped region, it is approximated by the capacitance between two unparalleled plates

C uz (q) = (Le-le) cos(α) 0 r 0 w tan(α) 2t + 2x tan(α) dx ( 10
)
As a result, the total capacitance is

C(q) = C z (q) + C uz (q) (11)
Then including the kinetic energy of the load and the potential energy related to the gravity, the Hamiltonian function, i.e. the total energy of the system can be described as:

H(q, p, Q) = H s (q) + H c (q, Q) + 1 2 p 2 m + mgq ( 12
)
where the p is the momentum of the load and g is the gravity constant. By choosing the energy variables as the state variables x = [q, p, Q] T . The co-energy variables e are:

e =    ∂H ∂q ∂H ∂p ∂H ∂Q    =    ∂Hc ∂q + ∂Hs ∂q + mg p m Q C(q)    ( 13
)
where ∂H ∂q stands for the total force working on the load in the q direction, resulting from the sum of the spring force working on the load ∂Hs ∂q , the force from capacitance ∂Hs ∂q and the gravity force mg.

With the above variables, the dynamics of the actuator can be presented by the following PH formulation:

  q ṗ Q  =   0 1 0 -1 -b(q) 0 0 0 -1/r   (J-R)    ∂H ∂q ∂H ∂p ∂H ∂Q    +   0 0 1/r   g V ( 14
)
where the interconnection matrix J = -J T represents the energy exchanges in the system, while the damping matrix R = R T = diag 0 b(q) 1/r ≥ 0 describes the internal dissipation of the system with b(q) the nonlinear damping coefficient depending on the vertical and horizontal damper and r the resistance of the actuator. Let choose the current, the power conjugate variable of the input, as the output y = g T ∂H ∂x = i, the system is passive, since the Hamiltonian is such that H > 0 and H(0) = 0, moreover its time derivative satisfies:

Ḣ = - ∂H ∂x T R ∂H ∂x + y T u ≤ y T u. (15)
III. POSITION CONTROL DESIGN In this work, the objective is to control the end position of the HASEL in closed loop. Based on the PH model proposed in the previous section, the IDA-PBC method will be investigated to control the position of the actuator. The IDA-PBC method provides several degrees of freedom to assign the closed-loop equilibrium position and performance with clear physical interpretation from energy point of view [START_REF] Ortega | Interconnection and damping assignment passivity-based control of port-controlled Hamiltonian systems[END_REF]. Furthermore, a structure preserving integral action will be added to improve the robustness of the control law while guaranteeing the closed-loop stability.

A. Basics of IDA-PBC Design

For the design procedure of the IDA-PBC method, let consider the open loop system [START_REF] Ortega | Interconnection and damping assignment passivity-based control of port-controlled Hamiltonian systems[END_REF]:

ẋ = (J -R) ∂H ∂x + g(x)u. (16)
Define an asymptotically stable PH target system

ẋ = (J d -R d ) ∂H d ∂x (17)
where the matrices

J d (x) = -J d (x) T , R d (x) = R T d (x)
≥ 0 and the function H d that satisfy the following matching equation:

g ⊥ (J d -R d ) ∂H d ∂x = g ⊥ (J d -R d ) ∂H ∂x (18)
with g ⊥ a full rank left annihilator of g, i.e g ⊥ g = 0. H d (x) is the desired Hamiltonian (closed-loop energy) function such that x * = argmin H d (x), with x * the equilibrium to be stabilized. The closed-loop system [START_REF] Ferguson | New results on disturbance rejection for energy-shaping controlled port-Hamiltonian systems[END_REF] with the feedback law u = β(x), where

β(x) = (g T g) -1 g T ((J d -R d ) ∂H d ∂x -(J -R) ∂H ∂x ) (19)
behaves as the target system with the equilibrium x * asymptotically stable.

B. Application to HASEL actuator position control

To apply the IDA-PBC design procedure on the HASEL actuator position control, we firstly consider the open loop system [START_REF] Moretti | Analysis of dielectric fluid transducers[END_REF] and define the desired closed-loop Hamiltonian function

H d (p, q, Q) = k p 2 (q -q * ) 2 + 1 2 p 2 m + k Q 2 (Q -Q * ) 2 (20)
where q * and Q * are the desired equilibrium position and charge of capacitance, and the asymptotically stable PH target system corresponds to [START_REF] Ferguson | New results on disturbance rejection for energy-shaping controlled port-Hamiltonian systems[END_REF] with

J d =   0 J 12 α 1 -J 12 0 α 2 -α 1 -α 2 0   , R d =   0 0 0 0 r 1 0 0 0 1 r  
where α 1 , α 2 are the frees variable to be designed. By defining the full rank left annihilator g ⊥ = -1 0 0 0 1 0 , the matching equation ( 18) leads to

J 12 (q, p, Q) = 1 -α 1 m p k Q (Q -Q * ) r 1 (q, p, Q) = b + m p ( ∂H ∂q -J 12 k p (q -q * )- -α 2 k Q (Q -Q * )) (21)
which are the variable depending on the state variables q, p, Q.

The desired equilibrium of the system is x * = q * , 0, Q * . To make sure the closed loop system is asymptotically stable, the parameter r 1 ≥ 0 which implies

Q * ≥ α 2 k Q Q -κ α 2 k Q (22)
where κ = b + m p ( ∂H ∂q -J 12 k d (q -q * )) According to equation (19), the control law of closed-loop system is

β(x) = -α 1 rk p (q -q * )-k Q (Q-Q * )+ Q C(q) -α 2 r p m (23)

C. Integral Action (IA) for robustness

The main drawback of the controller proposed in the previous section is the lack of robustness to overcome the steady state error due to the external disturbance or the unmatched mass of load. To improve the robustness, the integral action control [START_REF] Ferguson | New results on disturbance rejection for energy-shaping controlled port-Hamiltonian systems[END_REF] based on the IDA-PBC method is applied in this subsection.

The uncertainties of the load can be interpreted as external force applied on the system. Hence, the closed-loop system in the previous section can be expressed as

  Q ṗ q   = J d (x) -R d (x)    ∂H ∂Q ∂H ∂p ∂H ∂q    +   0 d 0   ( 24
)
where the matrices J d (x) and R d (x) are defined as

J d (x) := J aa (x) J au (x) -J T au (x) J uu (x) =   0 -α 1 -α 2 α 1 0 -J 12 α 2 J 12 0   R d (x) := R aa (x) R au (x) -R T au (x) R uu (x) =   1 r 0 0 0 0 0 0 0 r 1  
and d is a unknown constant external force (disturbance) to the system that depends on the unknown load. According to [START_REF] Ferguson | New results on disturbance rejection for energy-shaping controlled port-Hamiltonian systems[END_REF], if we choose the new closed-loop Hamiltonian function

H cl H cl = H d + K i 2 Q -x c 2 , ( 25
)
we can derive the new closed-loop Hamiltonian system

  ẋa ẋu ẋc   = J cl (x) -R cl (x)    ∂H cl ∂xa ∂H cl ∂xu ∂H cl ∂xc    +     0 d 0 0     . ( 26
)
Thus the control law of the integral action is

u i = -J aa + R aa + J c1 -R c1 -R c2 ∂H d ∂x a ( 27
)
+ J c1 -R c1 K i (x a -x c ) + 2R au ∂H d ∂x u ẋc = -R c2 ∂H d ∂x a + (J au + R au ) ∂H d ∂x u (28)
where x c is the state of integral action. In this paper, we choose the following design of matrices:

J c1 = 0, R c1 = 1 r , R c2 = 0 (29)
and the interconnection and damping matrices given by

J cl :=   0 J au (x) + R au (x) 0 -(J T au (x) + R T au (x)) J uu (x) 0 0 0 0   R cl :=   1/r 0 1/r 0 R uu (x) 0 1/r 0 1/r   .
(30) Then the following integral control law and the controller state can be derived:

u i = 1 r K i (Q -x c ) (31) ẋc = -α 1 k p (q -q * ) -α 2 p m ( 32
)
where K i is the gain of the integral action need to be designed.

Combining the IDA-PBC (23) and Integral action (31), we can implement the controller as the block diagram in Fig. 4.

HASEL actuator

Integral action Eq. ( 32) and (33) q,p,Q q* IDA-PBC β(x) Eq. ( 24

A. Experiment Setup

The experimental setup description is presented as shown in Fig. 1 in Introduction . The real-life experimental setup is shown in Fig. 5 and composed of the following elements:

• The actuator used in this study is composed of 5 C-5015 HASEL actuators that were manufactured by Artimus Robotic ® . A 3-D printed load container is attached to the bottom of the actuator to simplify the load change during the manipulation. • A T rek ® 610E high voltage amplifier is used to linearly amplify 0-10V input signal to 0-10KV output voltage with guaranteeing the current up to 2mA which provides sufficient driven power for the actuator motion. This amplifier can also return the current measurement in realtime during the manipulation. • A Kenyence ® LK-G152 laser sensor is used to measure the actuator displacement. The sensor is tuned to have 10kHz of bandwidth and ±4cm measurement range which are large enough to track the dynamics of the actuator and to validate the proposed control law. • A computer with M atlab Simulink ® is used to generate the reference signals, to implement the proposed controller and to get the measurement data. • A dSP ACE board (dS1104) serves the signals converters interface such as DAC (digital analogic converter) between the computer and the actuator and the ADC (analogic digital converter) between experimental sensor measurement (displacement and current) and the computer. The HASEL actuator has the self-sensing property for its strain measurement [START_REF] Acome | Hydraulically amplified self-healing electrostatic actuators with muscle-like performance[END_REF] which can be used for the controller implementation. The aim of this paper is to propose an accurate and reliable dynamic model and passivity based position control law for the HASEL actuator. In this case, we use a high precision laser displacement sensor instead of self-sensing measurement to avoid the error which can occur in the self-sensing.

B. Model Identification and Validation

For the model identification and the experimental validation, we conduct experiments with different loads (m = 150g, 200g, 250g) and different applied voltages (V = 4kV, 5kV, 6kV ) to obtain the experimental data of the displacement q and the electrical charge Q. The objective of identification is to find the proper coefficients of the spring-damper systems which approximate the dynamics of dielectric liquid of the HASEL actuator. From the system (14), we can obtain the mechanical force of the HASEL actuator as F m = ∂Hs ∂q + F d = ṗ -∂Hc ∂q -mg with F d = b(q) p m . Thus the identification problem becomes to find the optimal solution for the original length, the stiffness of the springs and the damper coefficients such that the difference between the mechanical force F m and the measured force F e is minimum. The coefficients are nonlinear and defined as the first order polynomial form:

ξ i (m) = φ j 1 m + φ j 0 k i (q) = θ j 1 q + θ j 0 b i (q) = λ j 1 l i (q) + λ j 0 (33)
where j = v, h stands for the vertical and horizontal directions. The length of the springs is the function of the load m with coefficients φ j n , the stiffness and the damper are the function of deformation with coefficients θ j n and λ j n . The identification problem is solved by the nonlinear model identification ('nlgreyest') and the trust-region-reflective algorithm ('fmincon') which are implemented in Matlab Identification Toolbox ® .

The experimental data of the mechanical force used for the parameters identification are shown by the blue curve in Fig. 6. Using the identification toolbox as mentioned before, one can get the identified parameters. Then the model simulation result with the identified parameters is compared to the experimental data in Fig. 6. The curve fitting has a very high percentage (fit = 94%) where fit = 100 1 -yexp-ysim yexp-mean(yexp) % with the experimental data y exp and the model simulation result y sim .

In order to verify the previous identified parameters, the actuator position response in the simulation is also compared to the experimental position data as shown in Fig. 7 with high fitness (88%). Meanwhile, the model with the identified parameters provides high fitness both in position and force response in the case of the different loads and the different applied voltages. These results suggest that the proposed model is able to describe the dynamics of the HASEL actuator.

V. CONTROL IMPLEMENTATION

In this section, we implement the proposed IDA-PBC controller (23) and the Integrate action (31) to the actuator with identified parameters to get the desired positions. Fig. From Fig. 8, the closed loop responses with and without IA controller both follow the reference in a satisfactory way and the the difference is not easy to see. The relative error (e r = absolute error reference ) of the closed-loop responses are 3% to 5% with and without IA controller. In this paper, the controller parameters tuning rules are not investigated. The investigation on PID-PBC controller parameters tuning for a class of linear (linearized) port Hamiltonian mechanical systems can be found [START_REF] Chan-Zheng | Tuning rules for a class of passivity-based controllers for mechanical systems[END_REF]. The IDA-PBC controller parameters tuning rules of proposed controllers in this paper shall be investigated in the future.

In Fig 9, we use a sinusoidal reference signal to define the desired position. The signal frequency is 0.2Hz and the signal range is [0.5mm, 2mm] (black dotted line). In this case, we only use the IDA-PBC + IA controller to control the actuator. The experimental result shows the closed loop response (red solid line) follows well the reference (upper figure). The maximal error of the signal tracking is around 0.1mm and with 4% relative error.

In Figure 10, we show the effectiveness of the proposed IDA-PBC+IA controller to the frequency position assignment. The desired position is defined by the sinusoidal reference with 1 Hz, 3 Hz, 5 Hz, 7 Hz. We can see the proposed controller can guarantee the position assignment with different frequencies.

From Fig. 8, one can observe that the IA does not significantly improve the closed loop performance. But as mentioned before, the purpose of using IA is to improve the robustness of the controller to the unknown load disturbance. In Fig. 11, we show the robustness improvement with the IA (red solid line) compared to the IDA-PBC controller without IA (blue dotted line). We assign the desired position to 2mm and 0.5mm. A 100g load mass is added in the load container when the actuator position at 2mm and removed when the position at 0.5mm. we can see that the closed loop position is disturbed when we add or remove the load mass for both controllers. The IDA-PBC can not reject this external disturbance while the closed loop response with IDA-PBC+IA controller goes back to the desired position as expected. In this work, we propose a dynamic model of an HASEL actuator via the port-Hamiltonian approach. The mechanical behavior of the actuator is approximated by two perpendicular nonlinear spring-damper systems. The proposed model reproduces the main dynamic behavior of the actuator. The experimental results show that the proposed model is sufficiently accurate to cope with the main system dynamics (up to 94% fitness). Based on the proposed model, the IDA-PBC method is investigated for the positioning control of the actuator. Further integral action is added to cope with load uncertainties. The experimental closed loop responses validate the proposed control laws to different desired position references with the relative error less than 5%. This work is the first attempt to the control-oriented modeling and control design of the HASEL actuator taking the advantage of the PH approach. In the future, the dynamic modeling and control design for more complex structures based on the HASEL actuator will be investigated. The high speed position control design and the controller parameters tuning problem will be considered in the future.

Fig. 1 :

 1 Fig. 1: Experimental setup description

Fig. 2 :

 2 Fig. 2: (a) The geometry of HASEL unit; (b) Deformation description of the actuator

Fig. 3 :

 3 Fig. 3: (a) Initial state when load is placed at the end; (b) The deformed state with zipped part of electrodes.

Fig. 4 :Fig. 5 :

 45 Fig. 4: Block diagram of closed-loop control with integral action

Fig. 6 :Fig. 7 :

 67 Fig. 6: Mechanical force identification (Fitness: 94%)

Fig. 8 :Fig. 9 :

 89 Fig. 8: Position control with the successive ramp reference signal (Upper figure); Control voltage (Middle figure); Relative error e r (Bottom figure)

Fig. 10 :Fig. 11 :

 1011 Fig. 10: Position responses of the sinusoidal signal with 1 Hz, 3 Hz, 5 Hz, 7 Hz, respectively. )

TABLE I :

 I HASEL actuator's parameters

	Parameter	Value	Units
	Lp	0.012	m
	Le	0.06	m
	w	0.05	m
	t	18 × 10 -6	m
	0	8.85 × 10 -12	F/m
	r	2.2	F/m
	r	8000	Ω