AMS subjects:
###### Abstract

This paper is devoted to numerical simulation of a charged particle beam submitted to a strong oscillating electric field. For that, we consider a two-scale numerical approach as follows: we first recall the two-scale model which is obtained by using two-scale convergence techniques; then, we numerically solve this limit model by using a backward semi-lagrangian method and we propose a new mesh of the phase space which allows us to simplify the solution of the Poisson’s equation. Finally, we present some numerical results which have been obtained by the new method, and we validate its efficiency through long time simulations.

Two-scale semi-lagrangian simulation of a charged particle beam in a periodic focusing channel

Alexandre Mouton

IRMA, Université Louis Pasteur, F-67084 Strasbourg

#### AMS subjects:

35B27, 76X05, 65M25 (65Y20).

#### Keywords:

two-scale convergence, semi-lagrangian method, Vlasov-Poisson model.

## 1 Introduction

Recent papers have proved that the two-scale convergence theory developed by Allaire in [2] and Nguetseng in [17] can be used successfully in order to develop numerical methods for solving ODEs or PDEs with oscillatory singular perturbations: for example, Frénod, Salvarani and Sonnendrücker have developed a two-scale PIC method in [13] for simulations of charged particle beams in a periodic focusing channel, and Frénod, Mouton and Sonnendrücker have developed a two-scale finite volume method in [9] for solving the weakly compressible 1D Euler equations. We can also cite the work of Ailliot, Frénod and Monbet in [1] about the simulation of tide oscillation for long term drift forecast of objects in the ocean.
On the other hand, many papers devoted to the numerical simulation of Vlasov-type problems involve Particle-In-Cell methods (see Birdsall and Langdon [3]) or Eulerian methods like semi-lagrangian schemes (see Sonnendrücker et al.[19], Grandgirard et al.[14], Filbet and Sonnendrücker [7, 8], Cheng and Knorr [5]). Since papers like [9] or [13] can be viewed as parts of a work programme which goal is the development of two-scale numerical methods for simulations of magnetic confinement fusion, it can be interesting to couple two-scale convergence results on Vlasov models such as two-scale models obtained in [12] with a semi-lagrangian scheme.
However, it is preferable in a first time to develop such a method on a more simple problem in order to study its behavior especially in the context of non-smooth solutions. The context of charged particle beams in a periodic focusing external field described in [13] offers a relatively simple framework for answering these questions.
We recall that such a phenomenon can be successfully represented by the 3D Vlasov-Maxwell system. In the same spirit of [13], we will consider non-relativistic long and thin beams, so we can replace the full three-dimensional Vlasov-Maxwell system by its paraxial approximation. To obtain this approximation, we do the following assumptions:

• the beam has already reached its stationary state,

• the beam is long and thin,

• the beam is propagating at constant velocity along the longitudinal axis ,

• the beam is sufficiently long so we can neglect the longitudinal self-consistent forces,

• the external electric field is supposed to be -periodic in and independent of the time,

• the beam is axisymmetric,

• the initial distribution is concentrated in angular momentum.

Under all these assumptions, the 3D Vlasov-Maxwell system reduces itself to a 2D Vlasov-Poisson system where the longitudinal position coordinate plays the role of time, and even to a 1D axisymmetric Vlasov-Poisson system of the form

 (1.1)

where is the radial component of the position vector in the transverse plane to the propagation direction, is the projection of the transverse velocity in the transverse plan to the propagation direction, is the ratio between the characteristic transverse radius of the beam and the characteristic longitudinal length of the beam, is the distribution function of the particles, is the radial part of the transverse self-consistent electric field, is the radial part of the transverse external electric field defined with the dimensionless real constant and the tension functions and which are respectively constant and periodic. All these quantities and variables are dimensionless. This system is naturally defined for but we can extend it to by using the conventions , , and .

The aim the paper is to simulate the system (1.1) with a two-scale semi-lagrangian scheme when . Inspired by Frénod, Salvarani and Sonnendrücker [13], and Frénod, Mouton and Sonnendrücker [9], we derive the model (1.1) by using the two-scale convergence theory developed in Allaire [2] and Nguetseng [17] and we obtain a new model which is independent of . Then, instead of discretizing the model (1.1) with a classical semi-lagrangian scheme, we discretize the new model with a semi-lagrangian method in order to obtain an approximation of a function , where is the second time scale, which verifies . Proceeding in such a way presents the advantage that there is no longer -oscillations in the limit model, so a very small time step is no longer required in order to simulate these oscillations. However, since semi-lagrangian schemes are based on interpolation on a phase space mesh, we have to pay attention to the contributions of the second time scale in the limit model. For this reason, in the same spirit of the work of Lang et al.[16], we introduce a -dependent moving mesh adapted for this two-scale numerical method the goal of which is to reduce the number of interpolations during the simulation.

The paper is organized as follows: in section 2, we recall the procedure leading from the paraxial approximation of the complete 3D Vlasov-Maxwell system to the model (1.1), and we will recall some two-scale convergence results about this model. In section 3, we build the semi-lagrangian method on the limit model and we see how to simplify it by considering a particular mesh. Section 4 is devoted to numerical results obtained with the two-scale semi-lagrangian method: on one hand, we compare them to some results obtained from a classical semi-lagrangian method on the system (1.1) in terms of quality of results and CPU time, and on the other hand, we see in this section the consequences of the use of the new mesh in the same terms.

## 2 The two-scale model

Firstly, we recall in this section the way to obtain the system (1.1) from the paraxial approximation of the 3D Vlasov-Maxwell equations. Then, we present a theorem about the two-scale convergence of the solutions of (1.1) which has been proved by Frénod, Salvarani and Sonnendrücker in [13].

### 2.1 Scaling of the paraxial model

By applying the hypothesis about the considered beam which are mentioned in the introduction, the full three-dimensional Vlasov-Maxwell is reduced to

 vz∂zf+v⋅∇xf+qm(E+Ξ)⋅∇vf=0,f(x,z=0,v)=f0(x,v),−∇xϕ=E,−Δxϕ=ρε0=qε0∫R2fdv, (2.1)

where is the longitudinal position coordinate, is the transverse position vector, is the constant longitudinal speed, is the transverse speed vector, is the distribution function of particles whose charge is and mass is , is the potential linked with the transverse self-consistent electric field , and is the transverse external electric field. More details about this derivation can be found in [6] and [8]. In this model, we assume that the external electric field is given by

 Ξ(x,z)=−H0x+H1(ω1zl)x, (2.2)

where is a positive constant tension, is a -periodic tension and is a dimensionless real constant.

In the same spirit of Frénod, Salvarani and Sonnendrücker [13], we build a dimensionless version of the system (2.1)-(2.2) by introducing the dimensionless variables defined by

 x=λx′,z=Lz′,v=vzv′, (2.3)

where is the characteristic transverse radius of the beam and is the characteristic length of the beam. Moreover, we define the dimensionless quantities by

 f(λx′,Lz′,vzv′)=mε0q2λLf′(x′,z′,v′),E(λx′,Lz′)=mv2zqLE′(x′,z′),f0(λx′,vzv′)=mε0q2λLf′0(x′,v′),H0=¯¯¯¯¯¯¯H0H′0,ϕ(λx′,Lz′)=mλv2zqLϕ′(x′,z′),H1(τ)=¯¯¯¯¯¯¯H1H′1(τ). (2.4)

With these new variables, the system (2.1)-(2.2) can be rewritten under the form

 ∂z′f′+Lλv′⋅∇x′f′+[E′−¯¯¯¯¯¯¯H0λqLv2zmH′0x′+¯¯¯¯¯¯¯H1λqLv2zmH′1(ω1Lz′l)x′]⋅∇v′f′=0,f′(x′,z′=0,v′)=f′0(x′,v′),−∇x′ϕ′=E′,−Δx′ϕ′=∫R2f′dv′. (2.5)

Since the beam is supposed to be long and thin, it is natural to take the ratio

 λL=ϵ. (2.6)

Furthermore, as we want to simulate the beam over a large number of periods of the external electric field, we also consider the ratio

 lL=ϵ. (2.7)

Finally, we suppose that the external electric field is much stronger than the self-consistent electric field and that its oscillations in direction are of the same order as , so we consider that

 ¯¯¯¯¯¯¯H0λqLv2zm=1ϵ,¯¯¯¯¯¯¯H1λqLv2zm=1. (2.8)

Under all these hypothesis the system (2.5) reduces to

 ∂tfϵ+1ϵv⋅∇xfϵ+(Eϵ+Ξϵ)⋅∇vfϵ=0,fϵ(x,v,t=0)=f0(x,v),−∇xϕϵ=Eϵ,−Δxϕϵ=∫R2fϵdv,Ξϵ(x,t)=−1ϵH0x+H1(ω1tϵ)x, (2.9)

where the primed notations for the variables and the initial distribution have been eliminated, and where has been replaced by , by , by , by , and by .
We introduce the polar coordinates linked with by the relations

 x=rcosθ,vr=vxcosθ+vysinθ,y=rsinθ,vθ=vycosθ−vxsinθ. (2.10)

Since the beam is supposed to be axisymmetric, the system does not depend on . Furthermore, we assume that the initial distribution is concentrated in angular momentum. Then the system (2.9) reduces to (1.1).

### 2.2 Two-scale convergence results

Since the aim of the paper is to develop a two-scale numerical method in order to simulate the model (1.1), we have to establish that the solution of this model two-scale converges to a function in a certain Banach space, and we have to find a system of equations verified by . These results have been proved in [13] and are recalled in the theorem below.

###### Theorem 1.

We assume that and that the initial distribution of (1.1) satisfies the following properties:

• with ,

• for all ,

• .

Then, considering a sequence of solutions of (1.1) and extracting a subsequence from it, we can say that two-scale converges to and two-scale converges to . Furthermore, there exists a function such that

 F(r,vr,τ,t)=G(cos(τ)r−sin(τ)vr,sin(τ)r+cos(τ)vr,t), (2.11)

and is solution of

 ∂tG−[∫2π0sin(σ)Er(cos(σ)q+sin(σ)ur,σ,t)dσ]∂qG−[∫2π0sin(σ)IQ(ω1)2πH1(ω1σ)(cos(σ)q+sin(σ)ur)dσ]∂qG+[∫2π0cos(σ)Er(cos(σ)q+sin(σ)ur,σ,t)dσ]∂urG+[∫2π0cos(σ)IQ(ω1)2πH1(ω1σ)(cos(σ)q+sin(σ)ur)dσ]∂urG=0,G(q,ur,t=0)=12πf0(q,ur),1r∂r(rEr(r,τ,t))=∫RG(cos(τ)r−sin(τ)vr,sin(τ)r+cos(τ)vr,t)dvr, (2.12)

where is equal to 1 if , and 0 otherwise.

Of course, such a result exists for the solution of the system (2.9) and can be found in [13].

## 3 The two-scale semi-lagrangian method

In this section, we develop a two-scale semi-lagrangian method in order to approach the solution of (1.1), in the case where . As it has been suggested in the introduction, the strategy is to discretize the model (2.11)-(2.12) in order to obtain a good approximation of which can be used for approaching . As in the two-scale PIC-method developed by Frénod, Salvarani and Sonnendrücker in [13], there is an advantage by proceeding in such a way: since there is no longer -frequency oscillations in the system (2.12), we do not need a very small time step for good simulation. In a first time, we recall the basis of a semi-lagrangian method. Then we present a motivation for the development of a two-scale semi-lagrangian method through the description of a classical backward semi-lagrangian method on the model (1.1). Finally, we describe the two-scale numerical method itself and we suggest a new mesh in order to simplify it.

### 3.1 The semi-lagrangian method

In this paragraph, we recall the way to discretize the abstract model

 ∂tf(x,t)+U(x,t)⋅∇xf(x,t)=0 (3.1)

with a semi-lagrangian method. For this, we have to consider the characteristics of (3.1), which are the solutions of

 ∂tX(t)=U(X(t),t). (3.2)

It is an easy game to remark that is constant along the characteristics, i.e.

 ∂t(f(X(t),t))=0, (3.3)

so we can write

 f(X(t;x,s),t)=f(x,s), (3.4)

where is the solution (3.2) with the condition .
This property of is used in the semi-lagrangian method as follows: assuming that we know the value of at time on the mesh points , we use the property (3.4) to say that

 f(xi,tn+Δt)=f(X(tn−Δt;xi,tn+Δt),tn−Δt), (3.5)

so we have to compute the point first, then compute by interpolating on the points in order to obtain an approximation of . Sonnendrücker et al. have suggested in [19] a way to compute a second order approximation of : they discretize the equation (3.2) with a finite difference method in order to obtain the following approximation:

 X(tn−Δt;xi,tn+Δt)=xi−2di, (3.6)

where is solution of

 di=ΔtU(xi−di,tn). (3.7)

In many cases, is only known at points . Then we have to replace (3.7) by

 di=ΔtΠU(xi−di,tn), (3.8)

where is an interpolation operator on points . Assuming that the polynomial function is regular enough, we write the following expansion of (3.8):

 di=ΔtU(xi,tn)−Δt∇x(ΠU)(xi,tn)di+O(Δt2), (3.9)

because, in the expansion of in , we get a term, which is a . Then, we obtain a second order accurate approximation of given by

 di=Δt(Id+Δt∇x(ΠU)(xi,tn))−1×U(xi,tn). (3.10)

Considering this approach, the semi-lagrangian method writes:

1. knowing at time and at time , we compute by using the relation (3.10) for each ,

2. we compute at time as follows:

 f(xi,tn+Δt)=Πf(xi−2di,tn−Δt). (3.11)

### 3.2 Implementation of the non-homogenized model

In this paragraph, we describe a classical semi-lagrangian method on the system (1.1). Since the electric fields and do not depend on , we can do a time-splitting on the first equation of the model, i.e. solving separately at each time step the equations

 ∂tfϵ+vrϵ∂rfϵ=0, (3.12)

and

 ∂tfϵ+(Eϵr+Ξϵr)∂vrfϵ=0, (3.13)

with a second order in time numerical scheme instead of solving the complete equation with the same scheme. As a consequence, we only do 1D interpolations instead of 2D interpolations. Then, denoting with the aproximation of and with the approximation of on the uniform mesh where is the size of one cell in direction and is the size of one cell in direction, an iteration of the method is organized as follows:

1. Knowing and , we do a backward advection of in direction and we define by

 fϵ∗(ri,vrj)=Πvrfϵn(ri,vrj−Δt2(Eϵn(ri)+Ξϵr(ri,tn))), (3.14)

where is a 1D cubic spline interpolation operator on the points .

2. We do an advection of in direction and we define by

 fϵ∗∗(ri,vrj)=Πrfϵ∗(ri−Δtϵvrj,vrj), (3.15)

where is a 1D cubic spline interpolation operator on the points .

3. We compute by discretizing the formula

 Eϵn+1(ri)=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩1ri∫ri0∫Rsfϵ∗∗(s,vr)dvrdsif ri≠0,0otherwise, (3.16)

with the trapezoidal rule on the points .

4. We compute by doing a last advection of in direction:

 fϵn+1(ri,vrj)=Πvrfϵ∗∗(ri,vrj−Δt2(Eϵn+1(ri)+Ξϵr(ri,tn+1))). (3.17)

If we use such a method, we have to guarantee the accuracy of the scheme, especially we consider the -frequency oscillations of the external electric field. A solution is to assume that the time step satisfies

 ⎧⎪ ⎪ ⎪⎨⎪ ⎪ ⎪⎩ri−Δr≤ri−Δtϵvrj≤ri+Δr,vrj−Δvr≤vrj−Δt2(Eϵn(ri)+Ξϵ(ri,tn))≤vrj+Δvr, (3.18)

for all . Then, in order to obtain good results with this method, has to be of the same order of which penalizes the method in terms of CPU time cost when we consider a very small .

### 3.3 Implementation of the two-scale model

In this paragraph, we adapt the semi-lagrangian method we have described in the paragraph 3.1 to the model (2.12). In this case, the characteristics of the system are the solutions of

 {∂tQ(t)=⟨E1⟩(Q(t),Ur(t),t),∂tUr(t)=⟨E2⟩(Q(t),Ur(t),t), (3.19)

where and are defined by

 ⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩⟨E1⟩(q,ur,t)=−∫2π0sin(σ)[Er(cos(σ)q+sin(σ)ur,σ,t)+IQ(ω1)2πH1(ω1σ)(cos(σ)q+sin(σ)ur)]dσ,⟨E2⟩(q,ur,t)=∫2π0cos(σ)[Er(cos(σ)q+sin(σ)ur,σ,t)+IQ(ω1)2πH1(ω1σ)(cos(σ)q+sin(σ)ur)]dσ. (3.20)

As in the paragraph 3.1, we remark that the solution of (2.12) is constant along the characteristics, so we can write:

 G(q,ur,tn+1)=G(Q(tn−1;q,tn+1),Ur(tn−1;ur,tn+1),tn−1), (3.21)

where , and where is the solution of (3.19) with the condition .

Firstly, we define the mesh in and directions by considering the points and (, ). We also consider the uniform mesh on (). Finally, we fix a time step for the entire simulation. Then, denoting with the approximation of at time , an iteration of the semi-lagrangian method is organized as follows:

1. Assuming that we know the value of and on the mesh , we compute with the formula

 Er(qi,τm)=⎧⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪⎩1qi∫qi0∫RsGn(cos(τm)s−sin(τm)vr,sin(τm)s+cos(τm)vr)dsdvrif i≠0, 0otherwise. (3.22)

Since is only known at points , we have to interpolate . Assuming that the support of is included in with large enough, we use the trapezoidal rule in order to approach the integral above. We obtain

 Enr(qi,τm)≈ΔqΔur2Pur∑j=−Pur(Π2Gn(cos(τm)qi−sin(τm)urj,sin(τm)qi+cos(τm)urj)+2ii−1∑k=1kΠ2Gn(cos(τm)qk−sin(τm)urj,sin(τm)qk+cos(τm)urj)), (3.23)

where is a cubic spline interpolation operator on the points .

2. We compute and at points : since is only known at points , we have to interpolate . By using the trapezoidal rule for approximating the integrals in (3.20), we obtain

 (3.24)

where is a cubic spline interpolation operator on the points .

3. We compute the the shifts by using the following formula:

 (3.25)

where the matrix is defined by

 Ai,j=Id+Δt⎛⎜ ⎜⎝∂q(Π2⟨E1n⟩)(qi,urj)∂ur(Π2⟨E1n⟩)(qi,urj)∂q(Π2⟨E2n⟩)(qi,urj)∂ur(Π2⟨E2n⟩)(qi,urj)⎞⎟ ⎟⎠. (3.26)
4. We compute by interpolating at the points :

 Gn+1(qi,urj)=Π2Gn−1(qi−2d1i,j,urj−2d2i,j). (3.27)
5. We save the approximation of at time given by

 fϵ(r,vr,tn+1)∼2πΠ2Gn+1(cos(tn+1ϵ)r−sin(tn+1ϵ)vr,sin(tn+1ϵ)r+cos(tn+1ϵ)vr). (3.28)

In order to initialize this two time step advance, we have to compute from which is given as an initial data: for this purpose, we perform a complete iteration such as decribed above where is replaced by and where we assume that .

### 3.4 The two-scale mesh

In most of test cases, we assume that the support of the initial distribution is compact and is included in some for , large enough. Then, if we follow the algorithm presented in the previous paragraph, the first thing we have to compute is by approximating the integral in (3.22): for numerical reasons, we have to reduce the integral on to an integral on a compact interval. Furthermore, since we can only say that the support of is included in as it is illustrated in Figure 1, the integral on in (3.22) is reduced to an integral on . As a consequence, we have to do all the simulation on instead of in order to avoid losing some data and we have to increase the number of mesh points in order to keep good interpolation results, even if the distribution function will be equal to 0 at many new points.

Support of and with .

In this paragraph, we present a different approach in order to avoid this extension of the simulation domain and its mesh. Before explaining the main idea of this new method, we consider the following meshes on and :

 (3.29)

where , and . Considering the function defined by

 γ:R2×[0,2π]⟶R2(r,vr,τ)⟼(cos(τ)r−sin(τ)vr,sin(τ)r+cos(τ)vr), (3.30)

we define and by

 (3.31)

Mesh and support of .

Mesh and support of .

Mesh and support of .

The main idea of our new method is to compute at points and to compute at points , whereas the function is computed at points for every . This approach is similar as the time-dependent moving grid described by Lang et al. in [16], even if, in our case, the mesh only depends on and , and is completely defined before the beginning of the simulation.
As a first consequence, considering that the support of is included in is equivalent to considering that the support of is included in for any such as illustrated in Figures 2,3, and 4, where the support of a Kapchinsky-Vladimirsky distribution is represented (see [15], [18] or [8] for more details about this distribution). Then we do not have to extend in order to avoid losing some data. Furthermore, the equation (2.11) reads

 F(r,vr,τ,t)=G(γ(r,vr,τ),t). (3.32)

Considering this approach, we fix a time step for all the simulation. Then, an iteration of this new semi-lagrangian method is organized as follows:

1. Assuming that and are known on the mesh for all , we compute at points . With the new notations, the equation (3.22) simplifies itself to

 Enr(ri,τm)=⎧⎪ ⎪ ⎪ ⎪⎨⎪ ⎪ ⎪ ⎪⎩1ri∫ri0∫RsGn(γ(s,vr,τm))dvrdsif i≠0,0otherwise. (3.33)

Assuming that the support of is included in