Gridless particle technique for the Vlasov Poisson system in problems with high degree of symmetry
GGridless particle technique for the Vlasov-Poissonsystem in problems with high degree of symmetry
E. Boella , G. Coppa , A. D’Angola , and B. Peiretti Paradisi Center for Mathematical Plasma Astrophysics, KU Leuven,Leuven, Belgium Dipartimento di Elettronica e Telecomunicazioni, Politecnico diTorino, Torino, Italy Scuola di Ingegneria, Universit`a della Basilicata, Potenza, Italy Dipartimento Energia, Politecnico di Torino, Torino, ItalyJanuary 24, 2019
Abstract
In the paper, gridless particle techniques are presented in order tosolve problems involving electrostatic, collisionless plasmas. The methodmakes use of computational particles having the shape of spherical shellsor of rings, and can be used to study cases in which the plasma hasspherical or axial symmetry, respectively. As a computational grid isabsent, the technique is particularly suitable when the plasma occupies arapidly changing space region.
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 thepurpose of accelerating positive ions. Avoiding to introduce a computationalgrid is useful in situations (as for plasma expansions and explosions), in whichthe physical domain occupied by the particles increases rapidly in time. Inthis framework, in general situations one could employ a set of computationalparticles and directly calculate the electric field acting on each of them, as thesum of the contribution of the other particles. This requires an extremely highcomputational effort, unless the problem under exam presents some symmetry.In the work, the cases of spherical and axial symmetry are considered. In thefirst case (Sect. 2), the problem is essentially one dimensional and computationalparticles are in the shape of spherical shells. By using the Gauss’s formula, theelectric field is readily evaluated. For the second case (Sect. 3), particles are1 a r X i v : . [ phy s i c s . p l a s m - ph ] J a n odeled as thin circular rings, which are characterized by their radii and theiraxial coordinates. In this case, the evolution of the force acting on each particlerequires necessarily the calculation of the sum of contributions due to the otherparticles. Although some advantages which are present in the spherical case arelost, 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. This Section presents in a complete, rigorous way the method of the shells,which was already introduced and employed with different formulations by otherAuthors (in particular, in refs. [1, 2, 3]).
In its simplest formulation, a set of N computational particles is considered.After initializing their coordinates x i and momenta p i , the particle are orderedaccording to their radial coordinates r i = | x i | , so that r j > r i if j > i . Thenthe radial electric field acting on each particle is evaluated simply as: E i = i − (cid:88) j =1 q j + q i x i r i , (1)by using the Gauss’s formula and taking advantage of the spherical symmetry ofthe problem. The presence of the factor multiplying q i can be explained in asimple way by considering that, for r = r i − (cid:15) ( (cid:15) → + ) q i does not contribute tothe electric field, while for r = r i + (cid:15) the total charge to be evaluated is (cid:80) ij =1 q j .Thus, by supposing a linear behavior of E at the interface, the factor providesthe correct value of the field (a rigorous proof of the formula is presented in Sect.2.4). Finally, after evaluating E on each computational particle, the equationsof motion: d x i d t = p i m i , d p i d t = q i E i ( x , x , ..., x N ) , (2)can be solved by using a suitable numerical technique (e.g., the leapfrog orthe Runge-Kutta method), using a time step much smaller with respect to theinverse of the plasma frequency. The technique described above is very simple (for example, a MATLAB code canbe implemented in few lines of program), but it is excessively memory and timeconsuming, as it does not take fully advantage of the symmetry of the problem.2n fact, in a central field of forces, the trajectory of each particle takes place ona plane. Therefore, the motion is essentially a two-dimensional problem. Thisfact suggests a new, simpler formulation of the method. After generating theinitial 3D coordinates x i and momenta p i , a set of 2D coordinates X i and P i is defined as (cid:40) X i = ( r i , , i = 1 , , ..., N, P i = (cid:16) p i · x i r i , (cid:12)(cid:12)(cid:12) p i − (cid:16) p i · x i r i (cid:17) x i r i (cid:12)(cid:12)(cid:12)(cid:17) . (3)After that, the method is completely identical to the previous formulation, butit uses only 2D vectors. More in detail, the particles are ordered according tothe radial position R i = | X i | , the electric field is evaluated as E i = i − (cid:88) j =1 q j + q i X i R i , (4)and the evolution of the system is governed by the equations (cid:40) d X i d t = P i m i , d P i d t = q i E i ( X , X , ..., X N ) . (5) Starting form the Lagrangian L ( r, ϕ, ˙ r, ˙ ϕ, t ) = m (cid:0) ˙ r + r ˙ ϕ (cid:1) − q Φ ( r, t ) , (6)for a single particle in a central potential (Φ depends on t due to the interactionwith the other particles of the plasma), one can obtain the Hamiltonian H ( r, ϕ, p r , p ϕ , t ) = 12 m (cid:32) p r + p ϕ r (cid:33) + q Φ ( r, t ) , (7)and the equations of the motion d r d t = p r m , d ϕ d t = p ϕ mr , d p r d t = − q ∂ Φ ∂r + p ϕ mr , d p ϕ d t = 0 . (8)In other terms, as it is well known, for a central potential there is a constantof the motion, p ϕ , which corresponds to the axial angular momentum, and themotion in radial direction is essentially one-dimensional. This suggests a thirdway of studying the dynamics of these systems. Starting again from the set { x i , p i } one can calculate r i = | x i | , p r,i = p i · x i r i , p ϕ,i = r i (cid:12)(cid:12)(cid:12)(cid:12) p i − p r,i x i r i (cid:12)(cid:12)(cid:12)(cid:12) . (9)3hen, the radial electric field is evaluated as E r,i = i − (cid:88) j =1 q j + 12 q i r i (10)(of course, particles must be sorted according to r i ), and the equations of themotion assume the form: d r i d t = p r,i m i , d p r,i d t = q i E r,i ( r , r , ..., r N ) + p ϕ,i m i r i , (11)in which the p ϕ,i ’s are constants of the motion and they are fixed by the initialconditions. This last formulation is the most convenient in terms of memoryusage and computational effort. However, the presence of the term p ϕ / ( mr ) inEqs. (11) require a special care when r →
0. All things considered, the secondformulation represents a good compromise in terms of computational efficiencyand simplicity.
Due to symmetry, each computational particle can be regarded as a sphericalsurface (a “shell”) on which the electric charge is distributed uniformly. Thepoints on the surface move according to different trajectories, all sharing thesame radial coordinate, r ( t ), and the same angular momentum p ϕ . For simplic-ity, a system made of only two shells (having charge q and q and radii r and r , with r < r ) is considered now. As the electric field is given by E ( r ) = , r < r ,q r , r < r < r ,q + q r , r > r , (12)the electrostatic energy U can be readily evaluated, as U ( r , r ) = (cid:90) R E π d x = q r + q + 2 q q r . (13)If r is changed of δr , the change − δU of the energy is equal to the work qE · δr of the field on the shell itself. In other terms, one has: E = − q ∂U∂r = q r . (14)Similarly, the field acting on the second shell can be calculated as E = − q ∂U∂r = q + q r . (15)4n both cases, the value of the electric field is in agreement with the rule“ i − (cid:80) j =1 q j + 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 is always smaller than r and one hasd p d t = q q r , d p d t = q q + q r . (16)Here only radial motion is considered for simplicity (i.e., p ϕ = 0 for both shells).The two equations (16) can be also written as d p d t = − ∂∂r (cid:18) q r (cid:19) , d p d t = − ∂∂r (cid:18) q q + q r (cid:19) , (17)from which one immediately obtains p m + q r = Const ,p m + q q + q r = Const . (18)As the two shells continue to expand, the asymptotic kinetic energy for t → + ∞ , E (+ ∞ ), of the two shells can be readily evaluated, as E (+ ∞ ) = E (0) + q r (0) , E (+ ∞ ) = E (0) + q q + q r (0) . (19)Now, the case of collision is considered. When t = t c one has r ( t c ) = r ( t c ) = r c , and for t > t c the shell t < t c . For t > t c , Eqs. (16) must be replaced by d p d t = q q + q r , d p d t = q q r (20)(they are obtained by simply exchanging indices 1 and 2), from which one finallyobtains p m + q q + q r = Const ,p m + q r = Const . (21)5n the case of collision, in order to evaluate the new asymptotic energy, E (cid:48) (+ ∞ ),both Eqs. 18 (for t < t c ) and Eqs. 21 must be considered: E (cid:48) ( t c ) = E (0) + q r − q r c = E (+ ∞ ) − q r c , E (cid:48) ( t c ) = E (0) + q q + q r − q q + q r c = E (+ ∞ ) − q q + q r c , (22)and E (cid:48) (+ ∞ ) = E (cid:48) ( t c ) + q q + q r c = E (+ ∞ ) + q q r c , E (cid:48) (+ ∞ ) = E (cid:48) ( t c ) + q r c = E (+ ∞ ) − q q r c . (23)In other terms, the collision produces an increase ∆ E = q q /r c in the energy ofthe shell − ∆ E for the shell E of a shell is of the order of qQ/R , being Q thetotal charge and R the initial plasma radius. Being ∆ E ∼ q /R for a singlecollision, one can conclude that the “plasma parameter” ∆ E / E for a set on N shells will be of the order of q/Q = 1 /N . In practice, for typical values ofthe number of computational particles, the system can always be regarded ascollisionless. t N ( r < R ) / N Figure 1: Time evolution of the fraction of electrons inside the ion sphere fortwo different normalized temperature, T = 0 . , . T , ensemble averages (full black line) and standard deviation ranges (dashedblack lines) are reported for N = 10 shells and 300 simulations with differentinitial conditions, together with reference results provided by a simulation with N = 10 shells (dashed red line). 6 t N ( E < ) / N Figure 2: Time evolution of the fraction of trapped electrons for the same caseof Fig. 1.
Some typical results are reported in the following. In all the calculations, suit-able normalization for the physical quantities has been used such that the totalcharge, 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) theformation of shocks in Coulomb explosions [6]. Figures 1 and 2 refer to the earlystage of the electron expansion in a spherical plasma. It is assumed that elec-trons 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. Calculationshave been performed both with a reduced ( N (cid:39) ) and with a high numberof shells ( N (cid:39) ), in order to obtain reference results. The initial phase-spacedistribution of the electrons was generated by using random numbers, so for asmall number of particles the results will depend on the particular choice of po-sitions and velocities. For this reason, the same calculation has been repeatedfor 300 times (with different initial conditions, all corresponding to the samephysical situation) in order to obtain the mean behavior and the distributionof the physical quantities (as performed in [7]). In Figs. 1 and 2, the timeevolution of the number of electrons inside the ion sphere (i.e., with r ≤ R ) andof the fraction of trapped electrons (i.e., with total energy p m − e Φ( r ) ≤
0) arereported, respectively. As can be observed, the shell method provides excellentresults, even with a reduced set of particles.The second set of results (Figs. 3 and 4) refers to the acceleration of an7
10 20 30 4000.20.40.60.811.21.4 r v Figure 3: Phase-space distributions of a mixture with m /m = 2 / q = q at different times ( t = 3 ÷ m /m = 2 / , q = q ) are initially at rest and the ions are acceleratedby electrostatic repulsion. In Fig. 3 the phase-space distribution for the twospecies, calculated with the shell method and using 10 computational parti-cles, is reported at different times and compared with analytic results. Figure4 shows E ( t → + ∞ ) /m of the light ions as a function of their initial radialcoordinate, r = r . This curve is important in order to determine the asymp-totic energy spectrum, d N d E , of the ions (considering that ∆ E = d E d r ∆ r and∆ N = 4 πr n ∆ r ). The two figures show the excellent agreement between nu-merical and analytic results.The third case here considered concerns the shock formation in a Coulomb ex-plosion [4, 8]. The phenomenon arises when the initial ion distribution is notuniform, 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 consequentlyinner particles acquire higher kinetic energy with respect to the outer ones andcan “overtake” them. In the situation considered in Figs. 5 and 6, an ion plasmamade of only one species presents two regions with different density for t = 0.Figure 5 reports the value of the radial coordinate r ( r , t ) of the ions as a func-tion of their initial radius, r , for different times, while in Fig. 6 the phase-space8 .6 0.7 0.8 0.9 10.70.80.911.11.21.3 r E ( + ∞ ) / m Figure 4: E ( t → + ∞ ) /m of the light ions as a function of their initial radialcoordinate, r , for the case of Fig. 3. Results obtained with the shell method(blue dots) are compared with the analytic solution (red line).distribution is plotted. The results here reported show the ability of the shellmethod to analyze cases in which the density, in theory, may become infinite insome point; in fact, results obtained with a relative low (10 ) and with a verylarge (10 ) number of shells are in perfect agreement. In the case of axial symmetry the fundamental “brick” for a N -body techniqueis a ring. More precisely, tori having circular cross section (of radius a ) areconsidered here. The tori shares the same axis of symmetry (the z axis) and arecharacterized 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: U = 12 (cid:88) i (cid:54) = j q i q j ϕ ring ( R i , R j , z i − z j ) + N (cid:88) i =1 q i U torus ( R i , a ) , (24)where ϕ ring ( R, R (cid:48) , z (cid:48) ) is the potential generated by a unit charge distributed ona ring (i.e., a torus with a = 0) of radius R laying on the xy plane in a point ofpolar coordinates ( R (cid:48) , z (cid:48) ), while U torus ( R, a ) is the potential energy of a torus9 r r Figure 5: Radial coordinate, r , at different times ( t = 0 ÷ .
47) as a functionof their initial position, r , for a single-species ion plasma with a non uniforminitial density distribution. In the simulations, n ( r,
0) = n when r < R/ n when r ∈ (cid:2) R , R (cid:3) , with n /n = 8. Results for 10 shells (blue dots) arecompared with those obtained with 10 shells (red line).of unitary charge. The potential ϕ ring ( R ; R (cid:48) , z (cid:48) ) can be evaluated in terms ofthe complete elliptic integral of the first kind [9]: K [ x ] = (cid:90) π/ d α (1 − x sin α ) / , (27)as ϕ ring ( R ; R (cid:48) , z (cid:48) ) = 2 K [ ξ ] πs , (28) As a generic point of the ring has coordinates ( R cos( ϑ ) , R sin( ϑ ) ,
0) and the point wherethe potential has to be evaluated has coordinates ( R (cid:48) , , z (cid:48) ), the potential ϕ ring can be writtenas ϕ ring = 12 π (cid:90) π R + R (cid:48) + z (cid:48) − RR (cid:48) cos θ ) / d θ, (25)By introducing the new integration variable α = θ − π , the formula for ϕ ring becomes: ϕ ring = 2 qπ (cid:90) π/ (cid:2) ( R + R (cid:48) ) + z (cid:48) − RR (cid:48) sin α (cid:3) / d α, (26)from which Eq. (28) immediately follows. r v Figure 6: Ion phase-space distribution at different times ( t = 0 ÷ .
47) for thesame case of Fig. 5. Results for 10 shells (blue dots) are compared with thoseobtained with 10 shells (red line).being s = [( R + R (cid:48) ) + z (cid:48) ] / , ξ = 4 RR (cid:48) s . (29)The calculation of U torus ( R, a ) is reported in detail in the Appendix. Forthe case of interest in which a (cid:28) R , one has: U torus ( R, a ) ∼ − πR (cid:20) log (cid:16) a R (cid:17) − (cid:21) . (30)From Eq. (30), it can be noticed that U torus diverges for a →
0, and this is thereason why tori are considered and not simply rings. Instead, in calculating theinteraction energy between tori, the value of ϕ ring is employed, as it is supposedthat when a (cid:28) 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 towrite the Lagrangian of the system, the kinetic energy N (cid:88) i =1 m i (cid:16) ˙ R i + ˙ z i + R i ˙ ϕ i (cid:17) (31)must be considered. By introducing the momenta p R,i , p z,i , and p ϕ,i : p R,i = m i ˙ R i , p z,i = m i ˙ z i , p ϕ,i = m i R i ˙ ϕ i , (32)11 x a z R i z = z i Figure 7: Scheme of a torus.one finally obtains the Hamiltonian H of the N interacting rings as: H = N (cid:88) i =1 m i (cid:32) p R,i + p z,i + p ϕ,i R i (cid:33) + 12 (cid:88) i (cid:54) = j q i q j ϕ ring ( R i , R j , z i − z j )+ N (cid:88) i =1 q i U torus ( R i , a ) , (33)and the equations of the motion: d R α d t = p R,α m α , d z α d t = p z,α m α , d p R,α d t = p ϕ,α m α R α − (cid:80) β (cid:54) = α q α q β ∂∂R α ϕ ring ( R α , R β , z α − z β ) − q α ∂∂R α U torus ( R α ; a ) , d p z,α d t = − (cid:80) β (cid:54) = α q α q β ∂∂z α ϕ ring ( R α , R β , z α − z β ) . (34)The angular momenta p ϕ,α are constants of the motion. The partial derivativesof ϕ ring can be readily evaluated considering that:d K [ x ]d x = E [ x ] − (1 − x ) K [ x ]2 x (1 − x ) , (35)being E [ x ] = (cid:82) π/ (1 − x sin α ) / d α the complete elliptic integral of the secondkind [9]. Equations (34) have been deduced by considering only electrostaticinteraction in non relativistic limit. In principle, the method can be readilyextended to include relativistic particles and magnetic field (with axial symme-try). To test its accuracy, the ring method has been employed to simulate the12
10 20 30 40 5000.511.5 r v Figure 8: Phase-space distribution at different times ( t = 4 ÷
36) of a sphericalion plasma in the case of ring loading with method 1. Results obtained with thering method (blue dots) are compared with the analytic solution (red lines).expansion of an ion sphere of uniform density, for which a simple analytic solu-tion exists. The same normalizaion of the physical quantities of Sect. 2.5 is usedhere. The initial ring distribution { R i , z i } has been generated in two differentways: 1) by dividing the initial [ R, z ] domain (i.e., a half circle of radius R )into a number N of small squares, each corresponding to the cross section of aring; 2) by suitably taking a set of { R i , z i } in a random way in order to obtain auniform charge density. The radius a i of the section of each ring has been chosenas proportional to R i , i.e., a i = k · R i . The constant k has been determined byrequiring the potential energy of the set of the rings to be equal to the exactvalue of the energy of the sphere. Figures 8, 9 and 10, 11 refer to method 1 andmethod 2, for ring loading, respectively. In Figs. 8 and 9 the time evolution ofthe phase-space distribution, as obtained with the ring method, is shown andit is compared with its analytical behavior. Figures 10 and 11 show the totalkinetic energy of the ions, E = (cid:80) Ni =1 m i v i ( t ), as a function of t ; moreover, thebehavior of [ E ( t ) − E r ( t )] / E , where E r ( t ) = (cid:80) Ni =1 m i (cid:2) v i ( t ) · e r,i ( t ) (cid:3) is the kineticenergy due to the motion in radial direction, is also presented. Obviously, inthe exact solution E r ( t ) ≡ E ( t ), so a value of (cid:12)(cid:12) E−E r ( t ) E (cid:12)(cid:12) (cid:28)
10 20 30 40 5000.511.5 r v Figure 9: Same as Fig. 8 in the case of ring loading with method 2 ( t =4 ÷ t E −6 ( E − E r ) / E Figure 10: Time evolution of the total kinetic energy of the ions (blue line) andof the fraction of the perpendicular kinetic energy (red line) obtained with thering method for the same case of Fig. 8 (method 1 for ring loading). Resultsobtained with the ring method are compared with the analytic solutions (blackstars). 14 t E −3 ( E − E r ) / E Figure 11: Same as Fig. 10, using method 2 for ring loading. r z Figure 12: Particle trajectories for the Coulomb explosion of an ion plasmahaving initially a cylindrical shape (the ratio between initial radius R and height H is equal to 0.1) for t = 0 ÷
4. Results obtained with the ring method (bluelines) are compared with those obtained with the PIC method (red dotted lines).15
Figure 13: Angular distribution of E /m for the case of Fig. 12 for t = 0 ÷ . The agreement between the two techniques is excellent. The results presented in the paper and all the tests that have been performedprove the effectiveness of the numerical technique here proposed. The interac-tion between computational particles is not mediated by a grid and, as shown inSects. 2 and 3, the method can be deduced by using a Hamiltonian approach.Consequently, all the physical quantities of interest (e.g., momentum, energyand angular momentum) are conserved exactly by the method, and the only The code makes use of an (
R, z ) uniform grid that is expanding in order to follow themotion of the particles. Moreover, the electrostatic potential is calculated at the border of thecomputational 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. r z Figure 14: Same as Fig. 12, but for a cylinder with
H/R = 1 for t = 0 ÷ Figure 15: Same as Fig. 13, but for a cylinder with
H/R = 1 for t = 0 ÷ A Electrostatic energy of a torus with a (cid:28) R With reference to Figure 16, the electrostatic energy of a torus can be calculatedby dividing the cross section S in a large number of subdomains. Each of themgenerates an electrostatic potential that can be approximated as the one of aring. Indicating by ∆ q i the charge of the i -th subdomain and by ϕ ring ( x i ; x j )the potential in x i due to a unitary charge in x j , the energy of the torus can beapproximated byFigure 16: Cross section of a torus and coordinates employed in the calculation. U (cid:39) (cid:88) i (cid:54) = j ∆ q i ∆ q j ϕ ring ( x i , x j ) . (36)In the limit when the size of the subdomains tends to zero, one obtains U = (cid:90) S d x Q (cid:90) S d x P σ ( x Q ) σ ( x P ) ϕ ring ( x P , x Q ) , (37)where σ ( x ) is the charge density for a unit cross section. If the torus is uniformlycharged and if a (cid:28) R , one can assume σ (cid:39) qπa = Const . (38)18n order to evaluate ϕ ring ( x P , x Q ), the parameters s and ξ , defined in Eq. (29),must be evaluated. One has: ξ = 4( R + x ,P )( R + x ,Q ) s , s = (cid:2) ( R + x ,P + R + x ,Q ) + ( x ,P − x ,Q ) (cid:3) / . (39)It turns out useful to introduce the quantity η = R + x ,P + x ,Q , such that R + x ,P = η + x ,P − x ,Q , R + x ,Q = η − x ,P − x ,Q . In this way, ξ can bewritten as: ξ = 1 − (cid:16) x ,P − x ,Q η (cid:17) (cid:16) x ,P − x ,Q η (cid:17) (cid:39) − (cid:16) r P Q R (cid:17) , (40)with r P Q = ( x P − x Q ) . In fact, η is much larger with respect to | x ,P − x ,Q | ≤ a , so the approximation (cid:15) (cid:39) − (cid:15) can be used; moreover, η can beapproximated by R . Making use of the asymptotic behavior of K [ ξ ] for ξ → K [ ξ ] ∼ ξ → −
12 log(1 − ξ ) + log 4 , (41)and assuming that s (cid:39) R , the following expression for ϕ ring ( x P , x Q ) is ob-tained: ϕ ring ( x P , x Q ) = − πR log (cid:16) r P Q R (cid:17) . (42)Equation (42) can be employed in Eq. (36), which can be rewritten as U = σ (cid:90) S d x Q ϕ torus ( x Q ) , (43)being ϕ torus ( x Q ) = − πR (cid:90) S d x P log (cid:16) r P Q R (cid:17) . (44)For x Q = 0, ϕ torus is readily evaluated: ϕ torus (0) = − πR (cid:90) a πrdr log (cid:16) r R (cid:17) = − a R (cid:20) log (cid:16) a R (cid:17) − (cid:21) . (45)To calculate ϕ torus for a generic x Q ∈ S , one can start by noticing that log( r P Q )is proportional to the Green function for the two-dimensional Poisson’s equation: ∇ Q log r P Q = 2 πδ ( x Q − x P ) . (46)So, by applying the Laplacian operator ∇ Q to Eq. (44), one obtains ∇ Q ϕ torus = − πR (cid:90) S d x P · πδ ( x P − x Q ) = − R . (47)19ue to the symmetry of the problem, ϕ torus is a function of r Q = | x Q | , and theLaplacian operator can be written as ∇ Q = r Q dd r Q r Q dd r Q . Therefore, Eq. (47)can be immediately solved, so obtaining ϕ torus ( r Q ) = ϕ torus (0) − r Q R . (48)Finally, the energy of the torus can be calculated by using Eq. (43): U = q π a · π (cid:90) a r Q d r Q (cid:34) ϕ torus (0) − r Q R (cid:35) = − q πR (cid:20) log (cid:16) a R (cid:17) − (cid:21) . (49)Formula (49) is very accurate for a (cid:28) R . If compared with the value of U obtained from numerical integration of Eq. (36), the relative error is less than0.5% for a/R < .
2. A similar formula (without the term -1/4) has been deducedin a concise, brilliant way in [10] by using the technique of asymptotic matching.20 eferenceseferences