Giglio Eric

CPU-efficient numerical code for charged particle transport through insulating straight capillaries

Keywords:

est destinée au dépôt

I. INTRODUCTION

Since the pioneering works of Ikeda [2] and Stolterfoht [1], the guiding phenomenon of charged particles through nano or macro capillaries has attracted much attentions [3,4] and is still an active research field, as shown by recent publications [START_REF] Nguyen | Electron guiding in macroscopic borosilicate capillaries with large bending angles[END_REF][6][7]. On the theoretical side, simulations of beam particles (electrons, positrons or ions) through insulating capillaries provided quite early fruitful theoretical support to experimental results and are still improving the quality of their predictions. The numerical codes used by the authors for the simulations are usually home-made. The codes are typically segmented into two parts. One part calculates the trajectories of the inserted beam particles through the capillary until they escape or hit the inner capillary surface, depositing some charge at the impact point. The second part calculates the time evolution of the deposited charges in the capillary wall. Both parts are coupled by the electric field generated by the accumulated charges in the capillary.

There exist different approaches for the charge dynamics in insulating capillaries. On one hand, the charge that accumulates in the capillary is described by a set of N point charges. Depending on the model, the point charges are field driven along the surface or even into the bulk or simply sit at their impact point. The point charges then decay in time with a decay rate proportional to the bulk conductivity of the insulator. The electric potential at a given position is evaluated by summing over all (non-zero) point charges, V (⃗ r) = i 2 1+ϵr qi | ⃗ ri-⃗ r| , where a screening factor accounting for surface polarization is usually added. Such approaches were for example used in [8][9][10][11][12][13][14][15].

On the other hand, the accumulated charges are described by a surface charges represented on axis symmetric 2D grids at the inner and outer surfaces of the capillary. The charges that are deposited at the impact point are smeared over nearest grid points. The time evolution of both surface charge distributions is described by two coupled continuity equations. The electric field that drives the surface charges and the trajectory of the inserted beam particles is deduced from the free and induced polarization surface charges at the inner and outer surface of the capillary.

For the charge dynamics, the grid approach has some advantages over the point charge approach. First, boundary conditions at dielectric and conducting interfaces can easily be imposed on the electric potential and electric field, which allows calculating explicitly polarization charges that appear at the vacuum-dielectric interfaces and free charges that appear at the vacuum-metal interfaces. The electric field is evaluated self-consistently including the correct screening of the free charges at the inner surface. As a result, the screening term 2 1+ϵr , which anyway is a crude approximation here, is no longer necessary.

Second, the charge relaxation is controlled by the bulk and surface conductivity of the insulator, which may even depend non-linearly on the electric field. The charges that migrate from the inner to the outer surface accumulate at the outer surface and still contribute to the electric field. Absorbing boundary conditions at the outer surface can be added to permit the accumulated charges to decay, simulating a leakage current that allows to reproduce the measured charge decay time. The grid method allows thus to describe quite reliably the charge dynamics in insulating capillaries, which is crucial in order to get accurate beam particle trajectories and consequently accurate charge deposition. Reversely, as the transmitted particles monitor the time evolution of the electric field inside the capillary, one has also access to some aspects of the charge dynamics in insulators. If the comparison between experimental results and simulations are reliable enough, capillaries may thus also set to become a formidable tool to study the charge dynamics in insulators.

Last but not least, we consider the CPU cost for the charge dynamics. In the point charge approach, when the N point charges are field driven, the electric field needs to be evaluated at the position of each point charge, requiring thus N 2 evaluations. With increasing number of accumulated point charges, the charge dynamics may quickly became a bottleneck in the simulation. In the grid approach, the CPU cost for the charge dynamics is constant as the charge dynamics does not depend on the number of accumulated charges in the insulator but only on the number of grid points. The grid approach is thus well suited for studies where the insulator accumulates a large number of charges.

The grid approach was used for example to simulate the self-organized radial focusing of a 2.5 keV ion beam by a conical capillary [16,[START_REF] Giglio | Stabilizing the ion-beam transmission through tapered glass capillaries[END_REF]. Further it was used to show in how far secondary electrons, generated by ions hitting the capillary inner surface, can avoid Coulomb blocking of the transmitted beam [18,[START_REF] Giglio | Effect of secondary electrons on the patch formation in insulating capillaries by ion beams[END_REF]. In both cases it was necessary to simulate a large number of trajectories (about 10 7 ) and accumulate a large amount of charges in the capillary to study the time evolution of the transmitted beam.

While the surface grid method has many advantages, it may not be implemented as straightforwardly as the point charge method, which may explain why the surface charge dynamics on an axial 2D grid was not widely adopted by others. However, in the special case where the capillary is a straight cylinder and where the bulk and surface conductivities can be considered constant, one can obtain an implementation of the charge dynamics which is robust, simple and fast. In addition, making intense use of Fast Fourier Transform, one obtains also a fast numerical scheme for evaluating the electric field at a beam particle's position, which is relevant for the computation of the trajectory of the beam particles.

In this work we will present such an implementation, based on the spectral method. Instead of using a cylindrical 2D grid to describe the charge dynamics, all the calculations are done using the moments of the distributions. We emphasize that the numerical schemes proposed in this work are suited only for straight capillaries. For example for conical capillaries, a grid method for the surface charges is still preferable.

Surface charge dynamics on a grid has not been adopted by other authors in the field and one reason may be algorithmic/numerical complexity of such an approach. This work presents the numerical schemes that we use for simulating the guiding of ions through insulating straight nano or macro capillaries. The paper may be regarded as a walk through our numerical code where we will detail how we simulate the inserted beam and control its emittance, how we compute the charge deposition and charge dynamics in the insulating capillary, how we evaluate the electric field by accounting for imposed boundary conditions, and how we evaluate efficiently the trajectories of the inserted charged particles.

II. THEORETICAL MODEL

The theoretical model, on which the charge dynamics is based, have already been presented elsewhere [START_REF] Giglio | Charge relaxation rates in insulating straight capillaries[END_REF]. In this section we will merely present those parts of the model that are relevant for understanding the numerical schemes of the code.

A. Parameters of the system

Dimensions of the capillary

A straight insulating capillary has a length H, an inner radius R 1 and an outer radius R 2 , as depicted in Fig. 1. We label S 1 the inner vacuum-dielectric interface and S 2 the outer dielectric-vacuum interface. We further suppose that the capillary is surrounded by a metallic cylindrical surface S 3 of radius R 3 ≥ R 2 and length H. The interface S 3 shares the same symmetry axis than the capillary and is electrically grounded. Surface S 3 is an important addition because it imposes Dirichlet boundary conditions on the electric potential outside the capillary, which allows determining the electric potential unambiguously. In experiments, a grounded metal cylinder should therefore always surround the capillary, if a reliable comparison between experiments and theory is desired.

In the particular case where the outer surface (S 2 ) of the capillary is covered by a conducting paint and grounded like in [22,24], it is sufficient to put R 3 = R 2 . On contrary, if no grounded cylinder is present in the experiment, we recommend to use R 3 = H, which should be large enough so that S 3 does not influence significantly the electric field inside the capillary. In experiments, the front circular section of the capillary is grounded or mounted behind a grounded collimator plate with an aperture centered on the capillary entrance.

Electric properties of the capillary

The model assumes that the deposited charges do not accumulate in the bulk but are localized at the inner and outer interface of the capillary. This assumption should be valid for insulators where the electric conductivity is dominated by ionic conductivity, such as glasses. Indeed, when beam particles hit the interface S 1 , electrons are ejected from the impact point, or equivalently, holes are injected at the impact point. Those holes are considered to have a negligible mobility compared to alkali ions and the holes are mainly relaxed by ionic conductivity. And because ionic charge carriers cannot pass the interfaces S 1 and S 2 , a positive charge builds up at the interfaces, which is well approximated by surface charge densities at S 1 and S 2 . Our model assumes that the electric bulk and surface conductivities are constant in time and space. Consequently, field dependent electric conductivities are not supported.

The electric properties of the capillary are described by three parameters, the relative permittivity ε r , the bulk conductivity κ b and the surface conductivity κ s . Note that the latter is mainly due to absorbed impurities and surface defects. The author recommends therefore to use, if possible, values for the surface conductivity which were measured in secondary vacuum. But in general the surface conductivity is not a well-known quantity and its contribution to the charge relaxation rate is still an open question. Note that the surface conductivity on S 1 and S 2 are not necessarily equal and the reader may use different values for both.

Symmetries and boundary conditions

We assume here that the capillary axis is aligned with the z-axis of the cylindrical coordinate system (r, θ, z) and that the beam axis lies in the xOz plane. The deposited charge distributions and consequently the electric potential will thus have xOz plane symmetry.

At the interfaces S 1 and S 2 and S 3 , boundary conditions are imposed. At S 1 and S 2 , the normal components of the electric displacement field are discontinuous, with their difference equal to the respective free surface charge densities. At the grounded interface S 3 , the potential is zero.

ε r ε 0 E r (R + 1 ) -ε 0 E r (R - 1 ) = σ (1) (1)
ε 0 E r (R + 2 ) -ε r ε 0 E r (R - 2 ) = σ (2) (2) V (R 3 ) = 0 (3)
Equations (1-3) add surface polarization charges at both dielectric interfaces and free charges at the conducting interface S 3 . To determine the electric potential unambiguously, we also have to impose boundary conditions at z = 0 and z = H. The front circular section of the capillary is supposedly covered by a conducting layer or a collimator plate and grounded. We put thus the electric potential at z = 0 equal to zero. Further we assume an ohmic contact between the grounded layer and the front circular section of the capillary so that electrons can flow from the grounded layer to the front circular section of the capillary and neutralize the excess charges that gather at z = 0. Consequently, we impose an absorbing boundary condition at z = 0 for the surface charge densities σ (1) and σ (2) at the inner and outer surface respectively,

V = 0 σ (1,2) = 0 z = 0 (4)
At the rear circular section of the capillary (z = H), the reader can choose between an Dirichlet type or Neumann type boundary boundary condition.

1. If the rear circular section of the capillary is covered by a conducting layer and grounded, which is usually the case for nano-capillaries, we will chose a Dirichlet type boundary condition. In that case, the electric potential V and the surface charge densities σ (1,2) are set to zero at z = H. Imposing the charge density to be zero at the boundary implies that the charge densities are absorbed at the grounded rear circular section.

2. If the rear circular section of the capillary is electrically isolated, which is often the case for macrocapillaries, we may choose a Neumann boundary condition. In that case the z-component of the electric field ∂V ∂z and the z-component of the surface charge current ∂σ ∂z are set to zero at z = H. The latter implies that the charges cannot pass the boundary (blocking boundary).

The Neumann condition seems to be the natural choice if the rear circular section is not explicitly grounded. But if a positively charged capillary is not screened from stray electrons then the positive potential of the rear circular section may attract stray electrons until the potential at z = H is close to zero [23]. In that case, a Dirichlet boundary condition at z = H is probably better suited.

B. Spectral method

Discrete Fourier transforms (DFT) expand a function f (⃗ r) on orthogonal circular basis functions. In this Fourier representation, differentiation is straightforward and a linear differential equation with constant coefficients is transformed into an easily solvable algebraic equation. The inverse DFT allows then to transform the result back into the ordinary spatial representation. Such an approach is called a spectral method. One advantage of the spectral method is that the circular basis functions used in the expansion can be chosen to automatically satisfy imposed symmetries and boundary conditions on f (⃗ r). In the following we define the quantities that are relevant for the charge dynamics.

Surface charge distributions

In cylindrical coordinates (r, θ, z), the surface charge densities at the inner σ (1) and outer surface σ (2) are expanded on the following circular basis functions

σ (1,2) (θ, z, t) = M m=0 N n=1 σ (1,2) mn (t) cos(mθ) sin(k n z) . (5)
The orthogonal basis functions {cos(mθ)} account for the xOz plane symmetry of the problem, while the {sin(k n z)} functions satisfy automatically the absorbing boundary condition at the capillary inlet z = 0. The σ

k n = n π H , ∀ n ≥ 1 , (6)
then sin(k n H) = 0, which imposes an absorbing boundary condition at the rear circular section of the capillary.

If we define k n by

k n = (n - 1 2 ) π H , ∀ n ≥ 1 , (7)
then ∂ z sin(k n z)| z=H = 0, which imposes a zero-flux (blocking) boundary condition at z = H. In both cases, the {sin(k n z)} functions form an orthogonal basis set on the interval [0, H]. Consequently, knowing the moments σ

(1,2) mn (t) is equivalent to knowing the densities σ (1,2) (θ, z, t) at the interfaces of the tube, with the addition that the imposed boundary conditions at z = 0 and z = H are automatically satisfied.

Source distributions

When beam particles hit the surface, a charge is deposited at the impact point. We label γ (1) and γ (2) the deposited charge per unit time and unit area at the inner and outer surface, respectively. In the following we will refer to γ (1,2) as source distributions. They are expanded similarly,

γ (1,2) (θ, z, t) = M m=0 N n=1
γ (1,2) mn (t) cos(mθ) sin(k n z) . (8) inheriting thus also the xOz plane symmetry and boundary conditions properties of the circular basis functions.

Electric potential inside the capillary

Inside the capillary, 0 ≤ r ≤ R 1 , there are no charges and the electric potential satisfies the Laplace equation. The potential is thus expanded on harmonic basis functions. Retaining only those basis functions which wellbehave at r = 0, satisfy the xOz plane symmetry and imposed boundary conditions at z = 0 and z =, one has,

V (r, θ, z, t) = M,N m=0 n=1 v mn (t)I m (k n r) cos(mθ) sin(k n z) (9)
where I m () are modified Bessel functions of order m. Depending on the chosen definition of k n , the potential or the z-component of the electric field are zero at z = H. The terms v mn (t) depend linearly on the time-dependent surface charge moments [START_REF] Giglio | Charge relaxation rates in insulating straight capillaries[END_REF],

v mn (t) = a mn σ (1) mn (t) + a ′ mn σ (2) mn (t) .

The coefficients a mn and a ′ mn , are constants and depend merely on the parameters R 1 , R 2 , R 3 , H, ε r of the system. They account for polarization charges that appear at the two vacuum-dielectric interfaces (S 1 ) and (S 2 ) and at the vacuum-metal interface (S 3 ). In appendix V A, we give a Python code that generates the a mn and a ′ mn coefficients. Once those coefficients given, the potential V and the electric field -⃗ ∇V inside the capillary tube are readily computed using (10) and (9).

Strictly speaking, imposing Dirichlet or Neumann boundary conditions on the potential V at z = 0 and z = H introduces smalls errors on the electric field near the inlet and outlet of the capillary. For a Dirichlet boundary condition at z = 0 or z = H, the error can be corrected by adding to the potential V () the potential generated by a disk of radius R 1 , placed at the boundary and carrying the surface charge distribution ε 0 ∂V ∂z . For the Neumann boundary condition at z = H, the error can be corrected by adding to the potential V the electric potential generated by a disk of radius R 1 , placed at the boundary and carrying the potential -V (r, θ, z = H)/2. Fortunately, the error is limited to the regions near the inlet and outlet of the capillary and, because of the large aspect ratio of the capillary, it does not affect sensibly the trajectories of the beam particles and can thus usually be ignored.

C. Time evolution of the surface charge moments

In the case where the charge dynamics takes place on the interfaces of a straight cylindrical capillary and where the surface and bulk conductivities are constants, the charge dynamics can actually be mapped onto a set of linear differential equations [START_REF] Giglio | Charge relaxation rates in insulating straight capillaries[END_REF],

  σ(1) mn (t) σ(2) mn (t)   = -   f (11) mn f (12) mn f (21) mn f (22) mn   Fmn   σ (1) mn (t) σ (2) mn (t)   ⃗ σmn(t) +   γ (1) mn (t) γ (2) mn (t)   ⃗ γmn(t) (11)
The time-independent coefficients f (ij) mn are rates and depend on the electric properties and dimensions of the capillary, R 1 , R 2 , H, κ b , κ s and ε r . The rates account for charge relaxation due to conductivity as well as for charge decay due to absorbing boundaries. Their expressions are given in [START_REF] Giglio | Effect of secondary electrons on the patch formation in insulating capillaries by ion beams[END_REF][START_REF] Giglio | Charge relaxation rates in insulating straight capillaries[END_REF]. Remarkably, there is no coupling between the different modes and ( 11) can be solved independently for each couple (m, n). The eigenvalues of F mn correspond to the two charge relaxation rates that characterize the straight capillary and are labeled respectively 1/τ (1) mn and 1/τ

(2) mn . Let δt be the time step in the numerical integration of the charge dynamics. We introduce the discrete time t i define by

t i = i × δt . ( 12
)
If we approximate the source term γ mn (t) in the time interval t i ≤ t ≤ t i+1 by its time-averaged value γ mn (t i ), see [START_REF] Giglio | Charge relaxation rates in insulating straight capillaries[END_REF], then the set of differential equations ( 11) can be solved analytically. For better readability we introduce the auxiliary quantities,

λ (i) mn = e -δt/τ (i) mn (13)
Λ (i) mn = τ (i) mn (1 -λ (i) mn ) (14)
where i = 1, 2. Note that λ mn and thus Λ mn depend on the time step δt. Further, we introduce the projectors P mn and Q mn , defined in [START_REF] Giglio | Charge relaxation rates in insulating straight capillaries[END_REF]. They are deduced from F mn and thus also time-independent. The density moments at time t i + δt are then given by the rather simple expression ⃗ σ mn (t i+1 ) = (λ (1) mn P mn + λ (2) mn Q mn ) ⃗ σ mn (t i ) + (Λ (1) mn P mn + Λ (2) mn Q mn ) ⃗ γmn (t i ) .

The time evolution of the density moments can thus be computed straightforwardly using a CPU inexpensive algorithm which is also unconditionally stable. The quantities λ

(1,2) mn , P mn and Q mn need to be computed just once for all. Given the parameters R 1 , R 2 , R 3 , H, κ b , κ s and ε r , as well as N and M , the python routine rates.py discussed in appendix V A, calculates the matrices F mn for each couple (mn) then diagonalizes each matrix F mn , extracts the two corresponding relaxation rates τ (1) mn and τ

(2) mn , and generates P mn and Q mn . With the latter given, implementing the charge dynamics is straightforward.

D. Computation of the source distributions

Inner surface

In the following, we calculate the time-averaged source terms γ(1,2) mn introduced in the previous section. We assume that all beam particles have the same charge q and that the inserted beam current I in is constant in time. When a particles hits the inner surface S 1 of the capillary at a grazing angle, a charge

q dep = (q + N se e) (16)
is deposited at the impact point (θ p , z p ), where N se is the average number of emitted secondary electrons per impact and e the elementary charge. We suppose that the deposited charge q dep is smeared over the interface according to a Gaussian distribution with a variance of ∆z = H/N along the z axis and a variance of ∆θ = π/M along the angle θ. The deposited point charge is thus replaced by a deposited surface charge,

σ dep θp,zp (θ, z) = q dep πR 1 ∆θ∆z e -(θ-θp) 2 ∆θ 2 e -(z-zp ) 2 ∆z 2 . ( 17
)
Let N δ ≥ 1 be the number of particles that are inserted into the capillary during the time step δt,

N δ = δt I in q with δt ≥ q I in . ( 18
)
Out of the N δ inserted particles, N p (t i ) hit the inner surface during the time t i and t i+1 . The deposited surface charge density per unit time, averaged over the time interval δt, reads γh (θ, z, t i ) = 1 δt

Np(ti) p=1 σ dep θp,zp (θ, z) (19)
where the bar symbol stands for the time average. Projecting [START_REF] Giglio | Effect of secondary electrons on the patch formation in insulating capillaries by ion beams[END_REF] on the basis functions {cos(mθ) sin(k n z)} yields the time-averaged moments γh mn (t i ), defined as

γh mn (t i ) = ∞ -∞ dθ ∞ -∞ γh (θ, z, t i ) cos(mθ) sin(k n z)dz = g mn × 1 N δ Np(ti) p=1 cos(m θ p ) sin(k n z p ) , (20)
where the factors

g mn = I in 1 + N se e q 4 πR 1 H e -1 4 (∆θ 2 m 2 +∆z 2 k 2 n ) (21)
are computed once for all for to save CPU time. Secondary electrons, emitted from the impact point, are field driven and eventually absorbed elsewhere at the the interface S 1 [START_REF] Giglio | Effect of secondary electrons on the patch formation in insulating capillaries by ion beams[END_REF]. The source term γse mn (t i ) accounts for the deposited secondary electrons distribution at time step t i . Summing both contributions yields the total source moments at the inner surface, γ( 1)

mn (t i ) = γh mn (t i ) + γse mn (t i ) ( 22
)

Outer surface

For the outer surface S 2 , we distinguish tree cases. If S 2 is covered by the conducting layer S 3 and grounded, R 2 = R 3 , no surface charges accumulate at the interface S 2 and we can safely put σ (2) mn (t i ) = 0 and γ (2) mn (t i ) = 0 (23) Equation ( 23) represents thus an absorbing boundary condition at S 2 .

If R 3 > R 2 and if S 2 is sufficiently screened from stray electrons, we can safely put γ( 2) mn (t i ) = 0. Otherwise, stray electrons may be attracted by the charged capillary, see [23], and a stray electron current per unit surface γ(2) (θ, z, t i ) must be defined at S 2 by the reader.

Computational cost

While ( 15) is unconditionally stable, the time step δt has nevertheless to satisfy the following condition. The source terms γ(t) should vary sufficiently slowly during the time step δt so that is can be approached by its timeaveraged value. That implies that the electric field inside the capillary should not vary significantly during the time step, which is ensured if δt is small compared to the characteristic relaxation times τ (1,2) mn and if the inserted charge per time step does not exceed a certain value Q δ ,

I in δt = qN δ ≤ Q δ , (24)
where the value of Q δ depends on the dimensions of the capillary. For nano-capillaries, we estimate that Q δ should not exceed several elementary charges. With inserted currents typically in the range of 10 -18 A, we take N δ = 1 and thus δt = q/I in .

For macro-capillaries, we estimate that Q δ should not exceed 10 -14 C. With inserted currents typically in the range of 10 -13 A, equation ( 24) yields a time step of the order of δt ∼ 0.1s and about N δ ∼ 10 5 inserted particles per time step δt. This is an unnecessarily large number of trajectories that need to be computed per time step δt to evaluate the time-averaged source term, see Eq. [START_REF] Giglio | Charge relaxation rates in insulating straight capillaries[END_REF]. In our simulations we therefore state that one calculated trajectory stands for y ≥ 1 physical particles. The number of trajectories N δ that need to be calculated per time step is thus reduced by the factor y

N δ = Q δ y q (25)
while equation ( 18) is replaced by

δt = N δ y q I in (26)
While larger values of y reduce the number of the trajectories to be computed per time step, too large values for y may introduce spurious effects due to low statistics in the evaluation of the time-averaged source terms γh mn (t i ) [START_REF] Giglio | Charge relaxation rates in insulating straight capillaries[END_REF]. We recommend to keep y below 5000. For inserted currents I in ≥ 10 -15 A, a thumb rule for y may be given by the relation

y = 10 -3 × I in q . (27)
The computational cost for the evaluation of the source moments γ(h) mn (t i ) per time step scales thus as N δ × (1 + N se ) × M × N , which is relatively inexpensive, especially that recursive relationships can be used to generate efficiently the M terms {cos(m θ p )} and N terms {sin(k n z p )} that appear in [START_REF] Giglio | Charge relaxation rates in insulating straight capillaries[END_REF], see Appendix V B. As a result, the computation of the moments γ(1,2) mn (t i ) is fast with the advantage of being simple, without the need of an additional grid mapping the interface S 1 and S 2 .

E. Computation of the electric field on a 3 D grid

Fast Fourier related transforms

Unlike charge densities and source terms, the electric field and electric potential will be be evaluated on a 3D cylindrical grid (r i , θ j , z k ) mapping the interior, r i ≤ R 1 , of the capillary. Along the angular and axial coordinates the grid point θ j and z k are uniformly distributed,

θ j = j 2π M , j ∈ [0, M -1] (28)
z k = k H N , k ∈ [0, N ] . (29)
This is motivated by the desire to use Fast Fourier related transforms to pass from the representation in the moment space to the grid points representation. Therefore N and M should be powers of 2. As an indicator, for capillaries with an aspect ratio H/R 1 of less the 50, the author uses N = 256 and M = 16 in his simulations. For larger aspect ratios, N may be doubled. Along the radial coordinate, no Fast Fourier related transform will be used and the radial grid is spaced nonuniformly, becoming denser with increasing r. This is motivated by the fact that the electric field varies faster when approaching the surface. Let L + 1 be the number of radial grid points. We propose to distribute the radial grid point according to the relation

r i = R 1 i L , i ∈ [0, L] (30)
The grid has thus Ω = (L + 1) × M × N grid points. The author typically uses for the radial grid L = 7, resulting in Ω ∼ 32000.

The potential and electric field on the (r i , θ j , z k ) grid points are computed by applying the following scheme.

In a first step we define for a given time t, the auxiliary moments ψ mn (r i ) and its radial derivative ψ ′ mn (r i )

ψ mn (r i ) = v mn (t)I m (k n r i ) (31) ψ ′ mn (r i ) = -v mn (t) I m (k n r) ∂r r=ri , (32)
where the v mn (t) are given by (10). In a second step, we pass from the moments of index n to the real space points z k , by applying a discrete sine (DST) or discrete cosine (DCT) transforms to the moments ψ and ψ ′ as indicated below.

ψ mn (r i ) DST ---→ ξ m (r i , z k ) ∀ i, m ψ ′ mn (r i ) DST ---→ ξ ′ m (r i , z k ) ∀ i, m k n × ψ mn (r i ) DCT ---→ χ m (r i , z k ) ∀ i, m
and store the results in auxiliary quantities ξ, ξ ′ and χ. In a third step, we pass from the angular moments (m) to the angular grid points (θ j ),

ξ m (r i , z k ) DCT ---→ V (r i , θ j , z k ) ∀ i, n (33)
ξ ′ m (r i , z k ) DCT ---→ E r (r i , θ j , z k ) ∀ i, n (34) m × ξ m (r i , z k ) DST ---→ E θ (r i , θ j , z k ) ∀ i, n (35)
-χ m (r i , z k ) DCT ---→ E z (r i , θ j , z k ) ∀ i, n (36)
Thank to the efficiency of the DST and DCT, the CPU cost for evaluating the potential V and the 3 component of the electric field (E r , E θ , E z ) on the 3D cylindrical grid scales like 7 × Ω(log(M ) + log(N )) .

Updating the electric field on the grid is relatively CPU expensive but fortunately it does not need to be carried out every time step. Only when a relative change in the moments v mn exceeds a certain amount ϵ > 0 does one need to re-evaluated the electric field at the grid points. Let t j be the last time when the electric field has been evaluated on the grid. If at a later t i > t j there is at least one couple of indexes (m, n) for which

|v mn (t i ) -v mn (t j )| > ϵ |v mn (t j )| (38)
then and only then will the electric field on the grid be updated with the new values taken at time t i . As a result, the frequency with which the field is updated on the grid decreases with increasing amount of accumulated charge in the capillary. The value of parameter ϵ is fixed by the reader. The author typically uses ϵ = 1%.

III. TRAJECTORIES A. Equation of motion

Particles of charge q and mass m, extracted by a potential V s have an initial velocity of u 0 = 2qV s /m. Here we assume that only the electric field generated by the deposited charge in the capillary wall acts on the charged particles. Other forces like the image charge force at the vacuum-dielectric interface may be added by the reader [7]. The trajectories of charged beam particles are calculated using classical Newtonian dynamics, which in its dimensionless form reads [START_REF] Giglio | Effect of secondary electrons on the patch formation in insulating capillaries by ion beams[END_REF]

d ⃗ ũ d t = -1 2 ∂ Ṽ ∂ ⃗ r d ⃗ r d t = ⃗ ũ , with ⃗ ũ = ⃗ u/u 0 ⃗ r = ⃗ r/H t = t (u 0 /H) Ṽ = (V /V s ) (39)
where the tilded quantities are dimensionless. Note that the equation of motion (EOM) is independent from the charge over mass ratio q/m of the beam particles. The trajectory of a beam particle is thus also independent from q/m and depends only on its initial position (x 0 , y 0 , 0) and the initial velocity vector ⃗ u. The initial conditions are sampled by a rejection method which will be discussed in the next section III B. Indeed, some care must taken when choosing the initial position and velocity vector of the beam particles, as the emittance of the injected beam was found to have in some cases a strong influence on the transmission rate [START_REF] Giglio | Stabilizing the ion-beam transmission through tapered glass capillaries[END_REF]. The initial conditions should therefore be sampled such as to reproduce at best the emittance of inserted beam in the experiments.

B. Initial phase-space conditions

The parameters of the beam that enter the simulation are:

• the injected intensity I in ,

• the extraction potential V s of the source,

• the mass m the particles,

• the charge q of the particles,

• the tilt angle β between the capillary axis and the beam axis • the divergence (opening angle) α of the beam,

• the radius R s of the virtual source and distance D of the virtual source from the capillary entrance.

With the above parameters, the sampling of the initial positions and velocities of the inserted beam particles is fully controlled. We remind that the entrance of the capillary is placed at z = 0. The components (x 0 , y 0 , 0) of the initial position and (u x , u y , u z ) of the initial velocity vector are sampled by a rejection method according to the following scheme, already discussed in [START_REF] Giglio | Stabilizing the ion-beam transmission through tapered glass capillaries[END_REF]. Let us assume that the profile of the beam in the plane z = 0 has a Gaussian shape, centered on the capillary aperture and with a full width half maximum (FWHM) labeled W . The charged particles are emitted from a virtual source, located at a distance D upstream from the capillary entrance. The ratio W/D = tan α ≃ α ≪ 1 may be regarded as the divergence of the beam.

The virtual source, which corresponds a the focal waist of the beam, is modeled by a disc of radius w s . We typically use w s = 1 mm. The positions (x s , y s ) of the particles on the virtual source disc are sampled uniformly. The u x and u y components are sampled according to a normal distribution P () having a width of ∆u = α u 0 ,

P (u i ) = √ π ∆u e -u 2 i /∆u 2 , i = x, y (40)
u z = u 0 (41)
Assuming that no electromagnetic field is present between the virtual source and the cappilary's entrance, the particles travel in a straight line until they reach the plane z = 0, defining initial coordinates

x 0 = x s + u x d u 0 (42)
y 0 = y s + u y d u 0 (43)
z 0 = 0 (44)
Only the trajectories that pass through the aperture,

x 2 0 + y 2 0 < R 2 1
, are retained and constitute the initial conditions for the equation of motion (39).

In the present model, the emittence of the beam that enters the capillary is thus controlled by three parameters, namely the distance D of the virtual source (focal waist) from the aperture, the radius w s of the virtual source disk and the divergence of the beam α, which are easily accessible in experiments. Eventually, the EOM is integrated numerically using a 4 th order explicit Beeman algorithm with adaptive step-size control, which was designed to allow for a high number of particles in molecular dynamics simulations [START_REF] Beeman | Some multistep methods for use in molecular dynamics calculations[END_REF].

C. Interpolation of the electric field at the particle's position

In order to compute the trajectory of a beam particle trough the capillary, one needs to evaluate the electric field at a large number of different positions. As the electric field is given at the Ω grid points of a 3D cylindrical grid mapping the inner space of a straight capillary, the electric field at given position can be obtained simply by a 3D interpolation. The simplest approach is the trilinear interpolation, where the 8 nearest points from a given position ⃗ r are used. Trilinear interpolation needs 7 calls to a one-linear interpolation function. We recommend however a tricubic interpolation scheme where the 64 nearest points are used. It needs 21 calls to a one-dimensional cubic interpolation function but will yield much better energy conservation along the trajectory than a trilinear scheme.

IV. CONCLUSION

We present a walk through our numerical code InCa4D where we show how to implement the surface charge dynamics at the inner and outer surface of a straight insulating capillary, without the use of cylindrical 2D grids. The advantages of such an approach with respect to point charges are highlighted. In particular, by satisfying the imposed boundary conditions, the electric field in InCa4D is evaluated accurately, yielding the correct charge relaxation and charge depletion rates of the deposited charges. The surface dynamics approach is suited for both macro and nano capillary but its CPU efficiency is mostly appreciable for macro capillaries where the number of deposited charges may be large.

We also show how to sample the initial positions and velocities of the beam particles as a function of the diameter of the virtual source disk, its distance to the capillary entrance and the divergence of the beam. Finally we show how to compute efficiently the trajectories of the beam particles through the capillary. The key idea is to map the inner space of the capillary by a cylindrical 3D grid, then use Fast Fourier related transforms to evaluate the electric field at the grid points. Then the field at the position of the particles is obtained using, for example, a tricubic interpolation scheme. This work hopefully shows that it is rather easy and advantageous to implement the charge dynamics at the inner and outer surface of a straight capillary using the moments of the distributions. In a future step, will extend the charge dynamics to the bulk, so the the charges can also accumulate in the bulk and not only at the interfaces, which should be relevant for those cases where the hole mobility cannot be neglected.

V. APENDIX

A. Numerical routine rates.py The Python3 routine rates.py calculates all the constant coefficients that are necessary to (i) simulate via equation (15) the charge dynamics in a straight insulating capillary and (ii) evaluate the electric potential inside the capillary using (9). The code is freely available on github.com.

Given the inputs R 1 , R 2 , R 3 , H, ε r , κ b , κ s , M, N and flag f bc , the routine rates.py calculates for all indexes 0 ≤ m < M and 1 ≤ n ≤ N

• the coefficients a mn and a ′ mn that appear in Eq. (10) and that allow to express the electric potential inside de capillary as a function of the moments of the charge distributions.

• the two characteristic relaxation times τ (1) mn and τ

(2) mn of the system.

• the projector P mn

FIG. 1 :

 1 FIG. 1: The three interfaces S1, S2, and S3 of the capillary setup. The conducting S3 interface is grounded.

  (1,2) mn (t) are the time dependent moments of the distributions, where m is the angular index and k n the wave number of index n. The definition of k n differs according to the imposed boundary condition at z = H. If we define k n by

Acknowledgments

The author would like to thank for financial support by the the French National Research Agency (ANR) in the P2N 2012 program (ANR-12-NANO-0008) for the PELIICAEN project.

The flag f bc indicates if blocking or absorbing boundary conditions are used at z = H . The projector Q mn is deduced from P mn via the relation

where I is the 2 × 2 identity matrix. The routine rates.py is based on the definitions found in [START_REF] Giglio | Charge relaxation rates in insulating straight capillaries[END_REF]. It creates an output file 'coefficients.txt' that contains the data ordered as follow:

-

a mn τ

( To generate efficiently the M terms {cos(m θ p )} that appear in [START_REF] Giglio | Charge relaxation rates in insulating straight capillaries[END_REF], the recursive relation

should be used, starting from cos(θ p ). For the computation of the N terms {sin(k n z p )}, we must distinguish between two cases: If k n = nπ/N then the recursive relation sin((n + 1)ξ) = 2 cos(ξ) sin(nξ) -sin((n -1)ξ) (49) with ξ = π N z p should be used. If k n = (n -1 2 )π/N , we first note that sin((n -

)ξ) = -cos(nξ) sin(ξ/2) + sin(nξ) cos(ξ/2) .

(50) We can now use the Moivre formula cos(nξ) + i sin(nξ) = (cos(ξ) + i sin(ξ)) n (51

to generate all the N terms {cos(nξ)} and {sin(nξ)} recursively, starting from sin(ξ) and cos(ξ).