Gridless particle technique for the Vlasov Poisson system in problems with high degree of symmetry

TL;DR

This paper introduces a gridless particle method using spherical shells and rings to efficiently simulate symmetric electrostatic plasmas without a computational grid.

Contribution

The paper presents a novel gridless particle technique tailored for symmetric Vlasov-Poisson problems, improving simulation flexibility in dynamic plasma regions.

Findings

Effective handling of rapidly changing plasma regions

Suitable for spherically and axially symmetric problems

Eliminates the need for a computational grid

Abstract

In the paper, gridless particle techniques are presented in order to solve problems involving electrostatic, collisionless plasmas. The method makes use of computational particles having the shape of spherical shells or of rings, and can be used to study cases in which the plasma has spherical or axial symmetry, respectively. As a computational grid is absent, the technique is particularly suitable when the plasma occupies a rapidly changing space region.

Equations

$${\mathbf{E}}_i=\left(\sum _{j=1}^{i-1}{q}_j+\frac{1}{2}{q}_i\right)\frac{{\mathbf{x}}_i}{r_i^3},$$

$$\{\begin{array}{l}\frac{\mathrm{d}{\mathbf{x}}_i}{\mathrm{d}t}=\frac{{\mathbf{p}}_i}{m_i},\\ \frac{\mathrm{d}{\mathbf{p}}_i}{\mathrm{d}t}=q_i{\mathbf{E}}_i({\mathbf{x}}_1,{\mathbf{x}}_2,...,{\mathbf{x}}_N),\end{array}$$

$$\{\begin{array}{l}{\mathbf{X}}_i=(r_i,0),i=1,2,...,N,\\ {\mathbf{P}}_i=\left({\mathbf{p}}_i\cdot \frac{{\mathbf{x}}_i}{r_i},\left|{\mathbf{p}}_i-\left({\mathbf{p}}_i\cdot \frac{{\mathbf{x}}_i}{r_i}\right)\frac{{\mathbf{x}}_i}{r_i}\right|\right).\end{array}$$

$${\mathbf{E}}_i=\left(\sum _{j=1}^{i-1}{q}_j+\frac{1}{2}{q}_i\right)\frac{{\mathbf{X}}_i}{R_i^3},$$

$$\{\begin{array}{l}\frac{\mathrm{d}{\mathbf{X}}_i}{\mathrm{d}t}=\frac{{\mathbf{P}}_i}{m_i},\\ \frac{\mathrm{d}{\mathbf{P}}_i}{\mathrm{d}t}=q_i{\mathbf{E}}_i({\mathbf{X}}_1,{\mathbf{X}}_2,...,{\mathbf{X}}_N).\end{array}$$

$$\mathcal{L}\left(r,\phi ,\dot{r},\dot{\phi },t\right)=\frac{m}{2}\left({\dot{r}}^2+r^2{\dot{\phi }}^2\right)-q\Phi \left(r,t\right),$$

$$\mathcal{H}\left(r,\phi ,p_r,p_{\phi },t\right)=\frac{1}{2m}\left(p_r^2+\frac{p_{\phi }^2}{r^2}\right)+q\Phi \left(r,t\right),$$

$$\{\begin{array}{l}\frac{\mathrm{d}r}{\mathrm{d}t}=\frac{p_r}{m},\frac{\mathrm{d}\phi }{\mathrm{d}t}=\frac{p_{\phi }}{m{r}^2},\\ \frac{\mathrm{d}{p}_r}{\mathrm{d}t}=-q\frac{\partial \Phi }{\partial r}+\frac{p_{\phi }^2}{m{r}^3},\frac{\mathrm{d}{p}_{\phi }}{\mathrm{d}t}=0.\end{array}$$

$$r_i=|{\mathbf{x}}_i|,p_{r,i}={\mathbf{p}}_i\cdot \frac{{\mathbf{x}}_i}{r_i},p_{\phi ,i}=r_i\left|{\mathbf{p}}_i-p_{r,i}\frac{{\mathbf{x}}_i}{r_i}\right|.$$

$$E_{r,i}=\left(\sum _{j=1}^{i-1}{q}_j+\frac{1}{2}{q}_i\right)\frac{1}{r_i^2}$$

$$\{\begin{array}{l}\frac{\mathrm{d}{r}_i}{\mathrm{d}t}=\frac{p_{r,i}}{m_i},\\ \frac{\mathrm{d}{p}_{r,i}}{\mathrm{d}t}=q_i{E}_{r,i}(r_1,r_2,...,r_N)+\frac{p_{\phi ,i}^2}{m_i{r}_i^3},\end{array}$$

$$E(r)=\left\{\begin{array}[]{cccccc}0,&r<r_{1},\\ \dfrac{q_{1}}{r^{2}},&r_{1}<r<r_{2},\\ \dfrac{q_{1}+q_{2}}{r^{2}},&r>r_{2},\end{array}\right.$$

$$U(r_1,r_2)={\int }_{{\mathbb{R}}^3}\frac{E^2}{8\pi }{\mathrm{d}}^3\mathbf{x}=\frac{q_1^2}{2r_1}+\frac{q_2^2+2q_1{q}_2}{2r_2}.$$

$$E_1=-\frac{1}{q_1}\frac{\partial U}{\partial r_1}=\frac{\frac{1}{2}{q}_1}{r_1^2}.$$

$$E_2=-\frac{1}{q_2}\frac{\partial U}{\partial r_2}=\frac{q_1+\frac{1}{2}{q}_2}{r_2^2}.$$

$$\begin{array}[]{cccccc}\dfrac{\mathrm{d}p_{1}}{\mathrm{d}t}=q_{1}\dfrac{\frac{1}{2}q_{1}}{r_{1}^{2}}&,&\dfrac{\mathrm{d}p_{2}}{\mathrm{d}t}=q_{2}\dfrac{q_{1}+\frac{1}{2}q_{2}}{r_{2}^{2}}.\end{array}$$

$$\{\begin{array}{l}\frac{\mathrm{d}{p}_1}{\mathrm{d}t}=-\frac{\partial }{\partial r_1}\left(\frac{\frac{1}{2}{q}_1^2}{r_1}\right),\\ \frac{\mathrm{d}{p}_2}{\mathrm{d}t}=-\frac{\partial }{\partial r_2}\left(\frac{q_1{q}_2+\frac{1}{2}{q}_2^2}{r_2}\right),\end{array}$$

$$\{\begin{array}{l}\frac{p_1^2}{2m_1}+\frac{\frac{1}{2}{q}_1^2}{r_1}=\text{Const},\\ \frac{p_2^2}{2m_2}+\frac{q_1{q}_2+\frac{1}{2}{q}_2^2}{r_2}=\text{Const}.\end{array}$$

$$\{\begin{array}{l}{\mathcal{E}}_1(+\infty )={\mathcal{E}}_1(0)+\frac{\frac{1}{2}{q}_1^2}{r_1(0)},\\ {\mathcal{E}}_2(+\infty )={\mathcal{E}}_2(0)+\frac{q_1{q}_2+\frac{1}{2}{q}_2^2}{r_2(0)}.\end{array}$$

$$\{\begin{array}{l}\frac{\mathrm{d}{p}_1}{\mathrm{d}t}=q_1\frac{q_2+\frac{1}{2}{q}_1}{r_1^2},\\ \frac{\mathrm{d}{p}_2}{\mathrm{d}t}=q_2\frac{\frac{1}{2}{q}_2}{r_2^2}\end{array}$$

$$\{\begin{array}{l}\frac{p_1^2}{2m_1}+\frac{q_1{q}_2+\frac{1}{2}{q}_1^2}{r_1}=\text{Const},\\ \frac{p_2^2}{2m_2}+\frac{\frac{1}{2}{q}_2^2}{r_2}=\text{Const}.\end{array}$$

$$\{\begin{array}{l}{\mathcal{E}}_1^{\prime }(t_c)={\mathcal{E}}_1(0)+\frac{\frac{1}{2}{q}_1^2}{r_1}-\frac{\frac{1}{2}{q}_1^2}{r_c}={\mathcal{E}}_1(+\infty )-\frac{\frac{1}{2}{q}_1^2}{r_c},\\ \\ {\mathcal{E}}_2^{\prime }(t_c)={\mathcal{E}}_2(0)+\frac{q_1{q}_2+\frac{1}{2}{q}_2^2}{r_2}-\frac{q_1{q}_2+\frac{1}{2}{q}_2^2}{r_c}={\mathcal{E}}_2(+\infty )-\frac{q_1{q}_2+\frac{1}{2}{q}_2^2}{r_c},\end{array}$$

$$\{\begin{array}{l}{\mathcal{E}}_1^{\prime }(+\infty )={\mathcal{E}}_1^{\prime }(t_c)+\frac{q_1{q}_2+\frac{1}{2}{q}_1^2}{r_c}={\mathcal{E}}_1(+\infty )+\frac{q_1{q}_2}{r_c},\\ {\mathcal{E}}_2^{\prime }(+\infty )={\mathcal{E}}_2^{\prime }(t_c)+\frac{\frac{1}{2}{q}_2^2}{r_c}={\mathcal{E}}_2(+\infty )-\frac{q_1{q}_2}{r_c}.\end{array}$$

$$U=\frac{1}{2}\sum _{i\ne j}{q}_i{q}_j{\phi }_{ring}(R_i,R_j,z_i-z_j)+\sum _{i=1}^N{q}_i^2{U}_{torus}(R_i,a),$$

$${\phi }_{ring}=\frac{1}{2\pi }{\int }_0^{2\pi }\frac{1}{(R^2+R^{\prime 2}+z^{\prime 2}-2R{R}^{\prime }\cos \theta {)}^{1/2}}d\theta ,$$

$${\phi }_{ring}=\frac{2q}{\pi }{\int }_0^{\pi /2}\frac{1}{{\left[(R+R^{\prime }{)}^2+z^{\prime 2}-4R{R}^{\prime }{\sin }^2\alpha \right]}^{1/2}}d\alpha ,$$

$$K[x]={\int }_0^{\pi /2}\frac{d\alpha }{(1-x{\sin }^2\alpha {)}^{1/2}},$$

$${\phi }_{ring}(R;R^{\prime },z^{\prime })=\frac{2K[\xi ]}{\pi s},$$

$$s=[(R+R^{\prime }{)}^2+z^{\prime 2}{]}^{1/2},\xi =\frac{4R{R}^{\prime }}{s^2}.$$

$$U_{torus}(R,a)\mathrel{\overset{}{\scalebox{2.0}[1.0]{$\sim$}}}-\frac{1}{2\pi R}\left[\log\left(\frac{a}{8R}\right)-\frac{1}{4}\right].$$

Full text

1 Introduction

The work investigates the possibility of using gridless particle techniques [1, 2] in the study of plasmas which are produced by laser-matter interaction with the purpose of accelerating positive ions. Avoiding to introduce a computational grid is useful in situations (as for plasma expansions and explosions), in which the physical domain occupied by the particles increases rapidly in time. In this framework, in general situations one could employ a set of computational particles and directly calculate the electric field acting on each of them, as the sum of the contribution of the other particles. This requires an extremely high computational effort, unless the problem under exam presents some symmetry. In the work, the cases of spherical and axial symmetry are considered. In the first case (Sect. 2), the problem is essentially one dimensional and computational particles are in the shape of spherical shells. By using the Gauss’s formula, the electric field is readily evaluated. For the second case (Sect. 3), particles are modeled as thin circular rings, which are characterized by their radii and their axial coordinates. In this case, the evolution of the force acting on each particle requires necessarily the calculation of the sum of contributions due to the other particles. Although some advantages which are present in the spherical case are lost, the technique here presented conserves interesting features also in this case. Results for both cases are shown and they are compared with exact calculations (when available) or with Particle-In-Cell simulations.

2 The shell method

This Section presents in a complete, rigorous way the method of the shells, which was already introduced and employed with different formulations by other Authors (in particular, in refs. [1, 2, 3]).

2.1 First formulation

In its simplest formulation, a set of $N$ computational particles is considered. After initializing their coordinates $\mathbf{x}_{i}$ and momenta $\mathbf{p}_{i}$, the particle are ordered according to their radial coordinates $r_{i}=|\mathbf{x}_{i}|$, so that $r_{j}>r_{i}$ if $j>i$. Then the radial electric field acting on each particle is evaluated simply as:

[TABLE]

by using the Gauss’s formula and taking advantage of the spherical symmetry of the problem. The presence of the factor $\tfrac{1}{2}$ multiplying $q_{i}$ can be explained in a simple way by considering that, for $r=r_{i}-\epsilon$ ($\epsilon\rightarrow 0^{+}$) $q_{i}$ does not contribute to the electric field, while for $r=r_{i}+\epsilon$ the total charge to be evaluated is $\sum_{j=1}^{i}q_{j}$. Thus, by supposing a linear behavior of $\mathbf{E}$ at the interface, the factor $\tfrac{1}{2}$ provides the correct value of the field (a rigorous proof of the formula is presented in Sect. 2.4). Finally, after evaluating $\mathbf{E}$ on each computational particle, the equations of motion:

[TABLE]

can be solved by using a suitable numerical technique (e.g., the leapfrog or the Runge-Kutta method), using a time step much smaller with respect to the inverse of the plasma frequency.

2.2 Second formulation

The technique described above is very simple (for example, a MATLAB code can be implemented in few lines of program), but it is excessively memory and time consuming, as it does not take fully advantage of the symmetry of the problem. In fact, in a central field of forces, the trajectory of each particle takes place on a plane. Therefore, the motion is essentially a two-dimensional problem. This fact suggests a new, simpler formulation of the method. After generating the initial 3D coordinates $\mathbf{x}_{i}$ and momenta $\mathbf{p}_{i}$, a set of 2D coordinates $\mathbf{X}_{i}$ and $\mathbf{P}_{i}$ is defined as

[TABLE]

After that, the method is completely identical to the previous formulation, but it uses only 2D vectors. More in detail, the particles are ordered according to the radial position $R_{i}=|\mathbf{X}_{i}|$, the electric field is evaluated as

[TABLE]

and the evolution of the system is governed by the equations

[TABLE]

2.3 Third formulation

Starting form the Lagrangian

[TABLE]

for a single particle in a central potential ($\Phi$ depends on $t$ due to the interaction with the other particles of the plasma), one can obtain the Hamiltonian

[TABLE]

and the equations of the motion

[TABLE]

In other terms, as it is well known, for a central potential there is a constant of the motion, $p_{\varphi}$, which corresponds to the axial angular momentum, and the motion in radial direction is essentially one-dimensional. This suggests a third way of studying the dynamics of these systems. Starting again from the set $\{\mathbf{x}_{i},\mathbf{p}_{i}\}$ one can calculate

[TABLE]

Then, the radial electric field is evaluated as

[TABLE]

(of course, particles must be sorted according to $r_{i}$), and the equations of the motion assume the form:

[TABLE]

in which the $p_{\varphi,i}$’s are constants of the motion and they are fixed by the initial conditions. This last formulation is the most convenient in terms of memory usage and computational effort. However, the presence of the term $p^{2}_{\varphi}/(mr^{3})$ in Eqs. (11) require a special care when $r\rightarrow 0$. All things considered, the second formulation represents a good compromise in terms of computational efficiency and simplicity.

2.4 Interaction between shells

Due to symmetry, each computational particle can be regarded as a spherical surface (a “shell”) on which the electric charge is distributed uniformly. The points on the surface move according to different trajectories, all sharing the same radial coordinate, $r(t)$, and the same angular momentum $p_{\varphi}$. For simplicity, a system made of only two shells (having charge $q_{1}$ and $q_{2}$ and radii $r_{1}$ and $r_{2}$, with $r_{1}<r_{2}$) is considered now. As the electric field is given by

[TABLE]

the electrostatic energy $U$ can be readily evaluated, as

[TABLE]

If $r_{1}$ is changed of $\delta r_{1}$, the change $-\delta U$ of the energy is equal to the work $qE_{1}\cdot\delta r_{1}$ of the field on the shell itself. In other terms, one has:

[TABLE]

Similarly, the field acting on the second shell can be calculated as

[TABLE]

In both cases, the value of the electric field is in agreement with the rule “${\sum\limits_{j=1}^{i-1}q_{j}+\tfrac{1}{2}q_{i}}$”, which was introduced previously.

Now the dynamics of the two shells is considered. If there is no crossing (i.e., no collisions) between shells, $r_{1}$ is always smaller than $r_{2}$ and one has

[TABLE]

Here only radial motion is considered for simplicity (i.e., $p_{\varphi}=0$ for both shells). The two equations (16) can be also written as

[TABLE]

from which one immediately obtains

[TABLE]

As the two shells continue to expand, the asymptotic kinetic energy for $t\rightarrow+\infty$, $\mathcal{E}(+\infty)$, of the two shells can be readily evaluated, as

[TABLE]

Now, the case of collision is considered. When $t=t_{c}$ one has $r_{1}(t_{c})=r_{2}(t_{c})=r_{c}$, and for $t>t_{c}$ the shell $\#$1 overtakes the shell $\#$2. Therefore, Eqs. (16-18) are valid only for $t<t_{c}$. For $t>t_{c}$, Eqs. (16) must be replaced by

[TABLE]

(they are obtained by simply exchanging indices 1 and 2), from which one finally obtains

[TABLE]

In the case of collision, in order to evaluate the new asymptotic energy, $\mathcal{E}^{\prime}(+\infty)$, both Eqs. 18 (for $t<t_{c}$) and Eqs. 21 must be considered:

[TABLE]

and

[TABLE]

In other terms, the collision produces an increase $\Delta\mathcal{E}=q_{1}q_{2}/r_{c}$ in the energy of the shell $\#$1, and a corresponding decrease $-\Delta\mathcal{E}$ for the shell $\#$2. In a typical plasma expansion, the energy $\mathcal{E}$ of a shell is of the order of $qQ/R$, being $Q$ the total charge and $R$ the initial plasma radius. Being $\Delta\mathcal{E}\sim q^{2}/R$ for a single collision, one can conclude that the “plasma parameter” $\Delta\mathcal{E}/\mathcal{E}$ for a set on $N$ shells will be of the order of $q/Q=1/N$. In practice, for typical values of the number of computational particles, the system can always be regarded as collisionless.

2.5 Results

Some typical results are reported in the following. In all the calculations, suitable normalization for the physical quantities has been used such that the total charge, the total mass of the plasma and the initial radius $R$ are all equal to 1. Three cases are considered: 1) the electron expansion in a spherical plasma [4]; 2) the expansion of a plasma made of a mixture of two ion species [5]; 3) the formation of shocks in Coulomb explosions [6]. Figures 1 and 2 refer to the early stage of the electron expansion in a spherical plasma. It is assumed that electrons and positive ions are initially distributed uniformly in a sphere of radius $R$. Initially, electrons have Maxwellian velocity distribution with temperature $T$ and positive ions are considered at rest during all the transient. Calculations have been performed both with a reduced ($N\simeq 10^{3}$) and with a high number of shells ($N\simeq 10^{6}$), in order to obtain reference results. The initial phase-space distribution of the electrons was generated by using random numbers, so for a small number of particles the results will depend on the particular choice of positions and velocities. For this reason, the same calculation has been repeated for 300 times (with different initial conditions, all corresponding to the same physical situation) in order to obtain the mean behavior and the distribution of the physical quantities (as performed in [7]). In Figs. 1 and 2, the time evolution of the number of electrons inside the ion sphere (i.e., with $r\leq R$) and of the fraction of trapped electrons (i.e., with total energy $\tfrac{p^{2}}{2m}-e\Phi(r)\leq 0$) are reported, respectively. As can be observed, the shell method provides excellent results, even with a reduced set of particles.

The second set of results (Figs. 3 and 4) refers to the acceleration of an ion plasma made of a mixture of two different species. In this case, analytic solutions for the problem exist [5] and can be used as a reference. The two species ($m_{1}/m_{2}=2/3,q_{1}=q_{2}$) are initially at rest and the ions are accelerated by electrostatic repulsion. In Fig. 3 the phase-space distribution for the two species, calculated with the shell method and using $10^{3}$ computational particles, is reported at different times and compared with analytic results. Figure 4 shows $\mathcal{E}(t\rightarrow+\infty)/m$ of the light ions as a function of their initial radial coordinate, $r=r_{0}$. This curve is important in order to determine the asymptotic energy spectrum, $\tfrac{\operatorname{d}N}{\operatorname{d}\mathcal{E}}$, of the ions (considering that $\Delta\mathcal{E}=\tfrac{\operatorname{d}\mathcal{E}}{\operatorname{d}r_{0}}\Delta r_{0}$ and $\Delta N=4\pi r_{0}^{2}n_{0}\Delta r_{0}$). The two figures show the excellent agreement between numerical and analytic results.

The third case here considered concerns the shock formation in a Coulomb explosion [4, 8]. The phenomenon arises when the initial ion distribution is not uniform, in particular if the inner density is larger respect to the outer one. In fact, in this case the electric field has a maximum inside the plasma region (while it depends linearly on $r$ if the ion density is constant) and consequently inner particles acquire higher kinetic energy with respect to the outer ones and can “overtake” them. In the situation considered in Figs. 5 and 6, an ion plasma made of only one species presents two regions with different density for $t=0$. Figure 5 reports the value of the radial coordinate $r(r_{0},t)$ of the ions as a function of their initial radius, $r_{0}$, for different times, while in Fig. 6 the phase-space distribution is plotted. The results here reported show the ability of the shell method to analyze cases in which the density, in theory, may become infinite in some point; in fact, results obtained with a relative low ($10^{4}$) and with a very large ($10^{6}$) number of shells are in perfect agreement.

3 The ring method

In the case of axial symmetry the fundamental “brick” for a $N$-body technique is a ring. More precisely, tori having circular cross section (of radius $a$) are considered here. The tori shares the same axis of symmetry (the $z$ axis) and are characterized by their radii, $R_{i}$, and axial coordinates, $z_{i}$ (as in Fig. 7). When $N$ tori are considered, the electrostatic energy of the system can be written as:

[TABLE]

where $\varphi_{ring}(R,R^{\prime},z^{\prime})$ is the potential generated by a unit charge distributed on a ring (i.e., a torus with $a=0$) of radius $R$ laying on the $xy$ plane in a point of polar coordinates $(R^{\prime},z^{\prime})$, while $U_{torus}(R,a)$ is the potential energy of a torus of unitary charge. The potential $\varphi_{ring}(R;R^{\prime},z^{\prime})$ can be evaluated111As a generic point of the ring has coordinates $(R\cos(\vartheta),R\sin(\vartheta),0)$ and the point where the potential has to be evaluated has coordinates $(R^{\prime},0,z^{\prime})$, the potential $\varphi_{ring}$ can be written as

$\varphi_{ring}=\frac{1}{2\pi}\int_{0}^{2\pi}\frac{1}{(R^{2}+R^{\prime 2}+z^{\prime 2}-2RR^{\prime}\cos\theta)^{1/2}}\operatorname{d}\theta,$

(25)

By introducing the new integration variable $\alpha=\frac{\theta}{2}-\frac{\pi}{2}$, the formula for $\varphi_{ring}$ becomes:

$\varphi_{ring}=\frac{2q}{\pi}\int_{0}^{\pi/2}\frac{1}{\left[(R+R^{\prime})^{2}+z^{\prime 2}-4RR^{\prime}\sin^{2}\alpha\right]^{1/2}}\operatorname{d}\alpha,$

(26)

from which Eq. (28) immediately follows. in terms of the complete elliptic integral of the first kind [9]:

[TABLE]

as

[TABLE]

being

[TABLE]

The calculation of $U_{torus}(R,a)$ is reported in detail in the Appendix. For the case of interest in which $a\ll R$, one has:

[TABLE]

From Eq. (30), it can be noticed that $U_{torus}$ diverges for $a\rightarrow 0$, and this is the reason why tori are considered and not simply rings. Instead, in calculating the interaction energy between tori, the value of $\varphi_{ring}$ is employed, as it is supposed that when $a\ll R$ the energy of two tori or two rings is essentially the same.

Now, the equations of the motion for the set of rings are derived. In order to write the Lagrangian of the system, the kinetic energy

[TABLE]

must be considered. By introducing the momenta $p_{R,i}$, $p_{z,i},$ and $p_{\varphi,i}$:

[TABLE]

one finally obtains the Hamiltonian $\mathscr{H}$ of the $N$ interacting rings as:

[TABLE]

and the equations of the motion:

[TABLE]

The angular momenta $p_{\varphi,\alpha}$ are constants of the motion. The partial derivatives of $\varphi_{ring}$ can be readily evaluated considering that:

[TABLE]

being $E[x]=\int_{0}^{\pi/2}(1-x\sin^{2}\alpha)^{1/2}\mathrm{d}\alpha$ the complete elliptic integral of the second kind [9]. Equations (34) have been deduced by considering only electrostatic interaction in non relativistic limit. In principle, the method can be readily extended to include relativistic particles and magnetic field (with axial symmetry). To test its accuracy, the ring method has been employed to simulate the expansion of an ion sphere of uniform density, for which a simple analytic solution exists. The same normalizaion of the physical quantities of Sect. 2.5 is used here. The initial ring distribution $\{R_{i},z_{i}\}$ has been generated in two different ways: 1) by dividing the initial $[R,z]$ domain (i.e., a half circle of radius $R_{0}$) into a number $N$ of small squares, each corresponding to the cross section of a ring; 2) by suitably taking a set of $\{R_{i},z_{i}\}$ in a random way in order to obtain a uniform charge density. The radius $a_{i}$ of the section of each ring has been chosen as proportional to $R_{i}$, i.e., $a_{i}=k\cdot R_{i}$. The constant $k$ has been determined by requiring the potential energy of the set of the rings to be equal to the exact value of the energy of the sphere. Figures 8, 9 and 10, 11 refer to method 1 and method 2, for ring loading, respectively. In Figs. 8 and 9 the time evolution of the phase-space distribution, as obtained with the ring method, is shown and it is compared with its analytical behavior. Figures 10 and 11 show the total kinetic energy of the ions, $\mathcal{E}=\sum_{i=1}^{N}\tfrac{m_{i}}{2}\mathbf{v}_{i}^{2}(t)$, as a function of $t$; moreover, the behavior of $[\mathcal{E}^{(t)}-\mathcal{E}_{r}(t)]/\mathcal{E}$, where $\mathcal{E}_{r}(t)=\sum_{i=1}^{N}\tfrac{m_{i}}{2}\left[\mathbf{v}_{i}^{2}(t)\cdot\mathbf{e}_{r,i}(t)\right]$ is the kinetic energy due to the motion in radial direction, is also presented. Obviously, in the exact solution $\mathcal{E}_{r}(t)\equiv\mathcal{E}(t)$, so a value of \bigl{\lvert}\tfrac{\mathcal{E}-\mathcal{E}_{r}(t)}{\mathcal{E}}\bigl{\lvert}\ll 1 is expected. All the numerical results presented in Figs. 8, 9, 10, 11 are in excellent agreement with the theory.

The second group of results here presented concerns the Coulomb explosion of an ion plasma having initially a cylindrical form. These are cases of practical interest, as they simulate the ion acceleration of the positive ions of a thin solid target after interaction with a ultra intense laser pulse. Two cases are considered, in which the cylinder has different aspect ratio. Figures 12 and 13 show the trajectories of the ions and the angular distribution of the kinetic energies for the first case. The same physical quantities are presented in Figs. 14 and 15 for the second case. In the Figures, the results of the ring method are compared with those obtained by using a PIC code developed by the Authors 222The code makes use of an $(R,z)$ uniform grid that is expanding in order to follow the motion of the particles. Moreover, the electrostatic potential is calculated at the border of the computational domain by summing the contributions due to all the rings; in this way, “exact” boundary conditions are provided for the solver of the Poisson’s equation.. The agreement between the two techniques is excellent.

4 Final considerations

The results presented in the paper and all the tests that have been performed prove the effectiveness of the numerical technique here proposed. The interaction between computational particles is not mediated by a grid and, as shown in Sects. 2 and 3, the method can be deduced by using a Hamiltonian approach. Consequently, all the physical quantities of interest (e.g., momentum, energy and angular momentum) are conserved exactly by the method, and the only errors are due to time discretization. This properly represents an important feature of the method. When the problem has the required degree of symmetry, the methods of shells and of rings can be usefully employed in two cases: 1) to obtain results making use of a simple, easy-to-implement code; 2) to have reference results to test more complex codes, in particular when the physical region occupied by the plasma grows dramatically during the simulation. For these reasons, in the Authors’ opinion the method can be regarded as a useful tool, in particular in the study of laser-plasma interaction.

Appendix A Electrostatic energy of a torus with $a\ll R$

With reference to Figure 16, the electrostatic energy of a torus can be calculated by dividing the cross section $S$ in a large number of subdomains. Each of them generates an electrostatic potential that can be approximated as the one of a ring. Indicating by $\Delta q_{i}$ the charge of the $i$-th subdomain and by $\varphi_{ring}(\textbf{x}_{i};\textbf{x}_{j})$ the potential in $\textbf{x}_{i}$ due to a unitary charge in $\textbf{x}_{j}$, the energy of the torus can be approximated by

[TABLE]

In the limit when the size of the subdomains tends to zero, one obtains

[TABLE]

where $\sigma(\textbf{x})$ is the charge density for a unit cross section. If the torus is uniformly charged and if $a\ll R$, one can assume

[TABLE]

In order to evaluate $\varphi_{ring}(\textbf{x}_{P},\textbf{x}_{Q})$, the parameters $s$ and $\xi$, defined in Eq. (29), must be evaluated. One has:

[TABLE]

It turns out useful to introduce the quantity $\eta=R+\tfrac{x_{1,P}+x_{1,Q}}{2}$, such that $R+x_{1,P}=\eta+\tfrac{x_{1,P}-x_{1,Q}}{2}$, $R+x_{1,Q}=\eta-\tfrac{x_{1,P}-x_{1,Q}}{2}$. In this way, $\xi$ can be written as:

[TABLE]

with $r_{PQ}^{2}=\left(\textbf{x}_{P}-\textbf{x}_{Q}\right)^{2}$. In fact, $\eta$ is much larger with respect to $|x_{2,P}-x_{2,Q}|\leq a$, so the approximation $\tfrac{1}{1+\epsilon}\simeq 1-\epsilon$ can be used; moreover, $\eta$ can be approximated by $R$. Making use of the asymptotic behavior of $K[\xi]$ for $\xi\rightarrow 1$:

[TABLE]

and assuming that $s\simeq 2R$, the following expression for $\varphi_{ring}(\textbf{x}_{P},\textbf{x}_{Q})$ is obtained:

[TABLE]

Equation (42) can be employed in Eq. (36), which can be rewritten as

[TABLE]

being

[TABLE]

For $\textbf{x}_{Q}=0$, $\varphi_{torus}$ is readily evaluated:

[TABLE]

To calculate $\varphi_{torus}$ for a generic $\textbf{x}_{Q}\in S$, one can start by noticing that $\log(r_{PQ})$ is proportional to the Green function for the two-dimensional Poisson’s equation:

[TABLE]

So, by applying the Laplacian operator $\nabla_{Q}^{2}$ to Eq. (44), one obtains

[TABLE]

Due to the symmetry of the problem, $\varphi_{torus}$ is a function of $r_{Q}=|\textbf{x}_{Q}|$, and the Laplacian operator can be written as ${\nabla_{Q}^{2}=\tfrac{1}{r_{Q}}\tfrac{\operatorname{d}}{\operatorname{d}{r}_{Q}}r_{Q}\tfrac{\operatorname{d}}{\operatorname{d}{r}_{Q}}}$. Therefore, Eq. (47) can be immediately solved, so obtaining

[TABLE]

Finally, the energy of the torus can be calculated by using Eq. (43):

[TABLE]

Formula (49) is very accurate for $a\ll R$. If compared with the value of $U$ obtained from numerical integration of Eq. (36), the relative error is less than 0.5$\%$ for $a/R<0.2$. A similar formula (without the term -1/4) has been deduced in a concise, brilliant way in [10] by using the technique of asymptotic matching.