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
IRMA, Université Louis Pasteur, F-67084 Strasbourg
35B27, 76X05, 65M25 (65Y20).
two-scale convergence, semi-lagrangian method, Vlasov-Poisson model.
Recent papers have proved that the two-scale convergence theory developed by Allaire in  and Nguetseng in  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  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  for solving the weakly compressible 1D Euler equations. We can also cite the work of Ailliot, Frénod and Monbet in  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 ) or Eulerian methods like semi-lagrangian schemes (see Sonnendrücker et al., Grandgirard et al., Filbet and Sonnendrücker [7, 8], Cheng and Knorr ). Since papers like  or  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  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  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 , 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
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 , and Frénod, Mouton and Sonnendrücker , we derive the model (1.1) by using the two-scale convergence theory developed in Allaire  and Nguetseng  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., 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 .
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
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  and . In this model, we assume that the external electric field is given by
where is a positive constant tension, is a -periodic tension and is a dimensionless real constant.
where is the characteristic transverse radius of the beam and is the characteristic length of the beam. Moreover, we define the dimensionless quantities by
Since the beam is supposed to be long and thin, it is natural to take the ratio
Furthermore, as we want to simulate the beam over a large number of periods of the external electric field, we also consider the ratio
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
Under all these hypothesis the system (2.5) reduces to
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
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  and are recalled in the theorem below.
We assume that and that the initial distribution of (1.1) satisfies the following properties:
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
and is solution of
where is equal to 1 if , and 0 otherwise.
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 , 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
with a semi-lagrangian method. For this, we have to consider the characteristics of (3.1), which are the solutions of
It is an easy game to remark that is constant along the characteristics, i.e.
so we can write
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
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  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:
where is solution of
In many cases, is only known at points . Then we have to replace (3.7) by
where is an interpolation operator on points . Assuming that the polynomial function is regular enough, we write the following expansion of (3.8):
because, in the expansion of in , we get a term, which is a . Then, we obtain a second order accurate approximation of given by
Considering this approach, the semi-lagrangian 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.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
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
where is a 1D cubic spline interpolation operator on the points .
We do an advection of in direction and we define by
where is a 1D cubic spline interpolation operator on the points .
We compute by discretizing the formula
with the trapezoidal rule on the points .
We compute by doing a last advection of in direction:
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
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
where and are defined by
As in the paragraph 3.1, we remark that the solution of (2.12) is constant along the characteristics, so we can write:
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:
Assuming that we know the value of and on the mesh , we compute with the formula
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
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
where is a cubic spline interpolation operator on the points .
We compute the the shifts by using the following formula:
where the matrix is defined by
We compute by interpolating at the points :
We save the approximation of at time given by
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.
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 :
where , and . Considering the function defined by
we define and by
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 , 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 ,  or  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
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:
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
Assuming that the support of is included in