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 twoscale numerical approach as follows: we first recall the twoscale model which is obtained by using twoscale convergence techniques; then, we numerically solve this limit model by using a backward semilagrangian 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.
Twoscale semilagrangian simulation of a charged particle beam in a periodic focusing channel
Alexandre Mouton
IRMA, Université Louis Pasteur, F67084 Strasbourg
AMS subjects:
35B27, 76X05, 65M25 (65Y20).
Keywords:
twoscale convergence, semilagrangian method, VlasovPoisson model.
1 Introduction
Recent papers have proved that the twoscale 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 twoscale 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 twoscale 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 Vlasovtype problems involve ParticleInCell methods (see Birdsall and Langdon [3]) or Eulerian methods like semilagrangian 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 twoscale numerical methods for simulations of magnetic confinement fusion, it can be interesting to couple twoscale convergence results on Vlasov models such as twoscale models obtained in [12] with a semilagrangian 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 nonsmooth 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 VlasovMaxwell system. In the same spirit of [13], we will consider nonrelativistic long and thin beams, so we can replace the full threedimensional VlasovMaxwell 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 selfconsistent 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 VlasovMaxwell system reduces itself to a 2D VlasovPoisson system where the longitudinal position coordinate plays the role of time, and even to a 1D axisymmetric VlasovPoisson 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 selfconsistent 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 twoscale semilagrangian 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 twoscale 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 semilagrangian scheme, we discretize the new model with a semilagrangian 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 semilagrangian 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 twoscale 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 VlasovMaxwell system to the model (1.1), and we will recall some twoscale convergence results about this model. In section 3, we build the semilagrangian 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 twoscale semilagrangian method: on one hand, we compare them to some results obtained from a classical semilagrangian 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 twoscale model
Firstly, we recall in this section the way to obtain the system (1.1) from the paraxial approximation of the 3D VlasovMaxwell equations. Then, we present a theorem about the twoscale 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 threedimensional VlasovMaxwell is reduced to
(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 selfconsistent 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
(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
(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
(2.4) 
With these new variables, the system (2.1)(2.2) can be rewritten under the form
(2.5) 
Since the beam is supposed to be long and thin, it is natural to take the ratio
(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
(2.7) 
Finally, we suppose that the external electric field is much stronger than the selfconsistent electric field and that its oscillations in direction are of the same order as , so we consider that
(2.8) 
Under all these hypothesis the system (2.5) reduces to
(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
(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 Twoscale convergence results
Since the aim of the paper is to develop a twoscale numerical method in order to simulate the model (1.1), we have to establish that the solution of this model twoscale 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 twoscale converges to and twoscale converges to . Furthermore, there exists a function such that
(2.11) 
and is solution of
(2.12) 
where is equal to 1 if , and 0 otherwise.
3 The twoscale semilagrangian method
In this section, we develop a twoscale semilagrangian 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 twoscale PICmethod 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 semilagrangian method. Then we present a motivation for the development of a twoscale semilagrangian method through the description of a classical backward semilagrangian method on the model (1.1). Finally, we describe the twoscale numerical method itself and we suggest a new mesh in order to simplify it.
3.1 The semilagrangian method
In this paragraph, we recall the way to discretize the abstract model
(3.1) 
with a semilagrangian method. For this, we have to consider the characteristics of (3.1), which are the solutions of
(3.2) 
It is an easy game to remark that is constant along the characteristics, i.e.
(3.3) 
so we can write
(3.4) 
where is the solution (3.2) with the condition .
This property of is used in the semilagrangian 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
(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:
(3.6) 
where is solution of
(3.7) 
In many cases, is only known at points . Then we have to replace (3.7) by
(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):
(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
(3.10) 
Considering this approach, the semilagrangian method writes:

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

we compute at time as follows:
(3.11)
3.2 Implementation of the nonhomogenized model
In this paragraph, we describe a classical semilagrangian method on the system (1.1). Since the electric fields and do not depend on , we can do a timesplitting on the first equation of the model, i.e. solving separately at each time step the equations
(3.12) 
and
(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:

Knowing and , we do a backward advection of in direction and we define by
(3.14) where is a 1D cubic spline interpolation operator on the points .

We do an advection of in direction and we define by
(3.15) where is a 1D cubic spline interpolation operator on the points .

We compute by discretizing the formula
(3.16) with the trapezoidal rule on the points .

We compute by doing a last advection of in direction:
(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
(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 twoscale model
In this paragraph, we adapt the semilagrangian 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
(3.19) 
where and are defined by
(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:
(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 semilagrangian method is organized as follows:

Assuming that we know the value of and on the mesh , we compute with the formula
(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
(3.23) where is a cubic spline interpolation operator on the points .

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 .

We compute the the shifts by using the following formula:
(3.25) where the matrix is defined by
(3.26) 
We compute by interpolating at the points :
(3.27) 
We save the approximation of at time given by
(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 twoscale 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.
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
(3.30) 
we define and by
(3.31) 
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 timedependent 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 KapchinskyVladimirsky 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
(3.32) 
Considering this approach, we fix a time step for all the simulation. Then, an iteration of this new semilagrangian method is organized as follows:

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
(3.33) Assuming that the support of is included in