An Accurate P^{3}M Algorithm for Gravitational Lensing Studies in Simulations
DDraft version February 18, 2021
Typeset using L A TEX twocolumn style in AASTeX63
An Accurate P M Algorithm for Gravitational Lensing Studies in Simulations
Kun Xu and Yipeng Jing
1, 2 Department of Astronomy, School of Physics and Astronomy, Shanghai Jiao Tong University, Shanghai, 200240, China Tsung-Dao Lee Institute, and Shanghai Key Laboratory for Particle Physics and Cosmology, Shanghai Jiao Tong University, Shanghai,200240, China (Received xxxx xx, xxxx; Revised xxxx xx, xxxx; Accepted xxxx xx, xxxx)
Submitted to ApJABSTRACTWe present a two-dimensional (2D) Particle-Particle-Particle-Mesh (P M) algorithm with an opti-mized Green function and adaptive softening length for gravitational lensing studies in N-Body sim-ulations. The analytical form of the optimized Green function ˆ G ( k ) is given. The softening schemes( S ) are studied for both the PM and the PP calculations in order for accurate force calculation andsuppression of the particle discreteness effect. Our method is two orders of magnitude more accuratethan the simple PM algorithm with the poor man’s Green function ( ∝ /k ) at a scale of a few meshcells or smaller. The force anisotropy is also much smaller than the conventional PM calculation. Wecan achieve a force accuracy better than 0.1 percent at all scales with our algorithm. When we applythe algorithm to computing lensing quantities in N-Body simulations, the errors are dominated by thePoisson noise caused by particle discreteness. The Poisson noise can be suppressed by smoothing outthe particle distribution, which can be achieved by simply choosing an adaptive softening length in thePP calculation. We have presented a criterion to set the adaptive softening length. Our algorithm isalso applicable to cosmological simulations. We provide a python implementation P3Mlens for thisalgorithm.
Keywords:
Gravitational lensing (670); N-body simulations (1083); Algorithms (1883) INTRODUCTIONNowadays gravitational lensing becomes a powerfultool to investigate the distribution of dark matter, bothfor dark matter halos and for larger scale structures(Treu 2010; Mandelbaum 2018). With the advent ofongoing or planed large imaging and spectroscopic sur-veys such as LSST (Ivezi´c et al. 2019), DESI (Dey et al.2019), HSC-SSP (Aihara et al. 2018) and PFS (Takadaet al. 2014), the statistical uncertainties of observationalmeasurement for gravitational lensing will be very small,but the accuracy will be limited by remaining system-atics in the observation.Accurate theoretical modeling will be required toproperly interpret the lensing observations. A widelyused way for modeling is to perform ray-tracing com-
Corresponding author: Yipeng [email protected] putation on the basis of high resolution N-body simula-tions (Wambsganss et al. 1998; Jain et al. 2000; Takada& Jain 2003a). In nearly all the ray-tracing methods,single or multiple plane approximation is used (Bartel-mann & Weiss 1994; Heymans et al. 2006; Fosalba et al.2008), where the 3D matter distribution around eachplane is projected onto the plane (the 2D distribution).Thus, the main step in the ray-tracing computation is toaccurately calculate the deflection angle of light for thelens plane, which can be reduced to a physical problemof solving a 2D Poisson’s equation.In the literature, 2D Poisson’s equation is usuallysolved using a Particle-Mesh (PM) algorithm. At first,mass of particles is assigned to a grid used for FastFourier Transform (FFT). Then, potential is obtainedby solving Poisson’s equation in Fourier space with aGreen function where ˆ G ( k ) ∝ /k is often used. Differ-ential is approximated by finite-difference to obtain theforce field or the lensing parameters on the grid fromthe potential, and finally the lensing quantities on the a r X i v : . [ a s t r o - ph . C O ] F e b Xu & Jing grid are interpolated to the whole space. Although theGreen function is correct analytically, because the calcu-lation is done on the grid with finite interval in the PM,errors are inevitably generated due to under-sampling(alias) and anisotropies in each step from mass assign-ment to force interpolation. This is usually called the poor man’s
Poisson Solver. Hockney & Eastwood (1981)provides a method to minimize the error in the wholePM process by regarding the Green function as a freeparameter and optimizing it. The force has to be soft-ened on small scale in the PM calculation even withthe optimized Green function. Thus, the PP algorithmis employed to compensate for the force softening in thePM calculation, and the force can be calculated to a veryhigh precision at all scale. The method is usually calledthe Particle-Particle-Particle-Mesh (P M) algorithm.Another source of error in ray-tracing simulations isPoisson noise. Since particle distribution in N-body sim-ulation is just a Monte Carlo sampling of the underlin-ing density distribution, Poisson noise is unavoidable.Smoothing is the main method to reduce the noise andrecover the density distribution, which has been investi-gated in many studies (Bradaˇc et al. 2004; Li et al. 2006;Yi-bin et al. 2009). Kernels of Gaussian or other shapeswith an adaptive smoothing length are usually used inPM algorithm when performing mass assignment. InP M algorithm, as we will show, the smoothing can befully incorporated in the PP calculation by transformingpoint sources to shaped sources. This is another advan-tage of the P M algorithm to predict lensing quanti-ties in N-Body simulations. Furthermore, we can easilyadopt an adaptive smoothing length in the algorithm.In this paper, we present a P M algorithm for lensingstudies in N-Body simulation. We study the accuracyof the computed lensing quantities, and investigate howthe errors are related to the accuracy of the force calcu-lation and to the Poisson noise of the discrete particledistribution. We give an analytical form of the opti-mized Green function and determine the best choice ofthe free parameters in PM. We also give an analyticalform for the force between 2D shaped particles in realspace, which is used both for the PP calculation andfor the smoothing of the particle distribution. Our rec-ommendation for the two softening (smoothing) lengthsis given. Finally, we will briefly introduce our python implementation
P3MLens of the 2D P M algorithm.This paper is organized as follows. In Section 2, webriefly introduce the gravitational lensing theory and thegeneral P M algorithm. In Section 3, our realization ofthe 2D P M algorithm is provided. In Section 4, wetest the accuracy of the algorithm. The python im-plementation is briefly introduced in Section 5 and our conclusions are summarized in Section 6. Throughoutthis paper, Fourier transform of a function f is denotedas ˆ f , and the convention of Fourier transform we use, asexample for one dimension, is f ( x ) = (cid:90) ∞−∞ dk π ˆ f ( k ) e ikx (1)ˆ f ( k ) = (cid:90) ∞−∞ dxf ( x ) e − ikx (2) METHOD2.1.
Gravitational lensing for thin lens
The deflection angle of light by a gravitational lens isthe key quantity calculated in the ray tracing simula-tions. For a thin lens, according to General Relativity(GR), the deflection angle α depends on the gravita-tional potential Φ as follows: α = 2 c (cid:90) ∇ Φ dl = 2 c ∇ (cid:90) Φ dl , (3)where ∇ is the 2D gradient operator in the lens plane.Using Poisson’s equation for 3D gravitational potential,we can easily prove, ∇ (cid:90) Φ dl = (cid:90) ∇ Φ − ∂ Φ ∂l dl = 4 πG (cid:90) ρ dl − ∂ Φ ∂l (cid:12)(cid:12)(cid:12)(cid:12) ∞∞ = 4 πG Σ . (4)In the above equation, Σ is the surface density of thelens. The second term on the rhs vanishes due to thezero boundary condition for the potential. It’s con-venient to define 2D potential ψ = (cid:82) Φ dl and field E = − ∇ ψ . We have α = − c E . Thus, the goalof ray-tracing simulations is to solve the 2D Poisson’sequation from the known matter distribution of an N-body simulation, ∇ ψ = 4 πG Σ (5)It’s easy to derive the 2D force field for a point source: E = − Gm r r , (6)where m is the mass of the point source.For a single thin lens, gravitational lensing parameterslike shear ( γ , γ ), convergence ( κ ) and magnification( µ ) can be calculated if α is known in the lens plane. γ = D l D ls D s ( ∂α x ∂x − ∂α y ∂y ) (7) M lensing γ = D l D ls D s ( ∂α x ∂y + ∂α y ∂x ) (8) γ = γ + iγ (9) κ = D l D ls D s ( ∂α x ∂x + ∂α y ∂y ) (10) µ = 1(1 − κ ) − | γ | (11)where D l , D s and D ls are the angular diameter distancesfrom the observer to the lens, the observer to the sourceand from the lens to the source.2.2. General multi-dimensional P M algorithm P M algorithms are a class of hybrid algorithms de-veloped decades ago to simulate plasma systems usingparticles. They have been used to simulate the evolutionof large-scale structures in the Universe (e.g. Efstathiouet al. (1985); Couchman (1991); Jing et al. (1994)). Theshort-range force on a particle is computed by directlysumming up particle-particle (PP) pair force, and thesmoothly varying long-range force is approximated bythe particle-mesh (PM) force calculation. In this subsec-tion, we briefly introduce the main idea of the P M algo-rithm and we recommend Hockney & Eastwood (1981)for details. The main steps of the P M algorithm in-clude: PM (mass assignment, potential solving, poten-tial difference and force interpolation) and PP (forcesplitting and summation) calculations.2.2.1.
Mass assignment
To take advantage of the Fast Fourier Transform(FFT) algorithm, the density field is assigned to a reg-ular grid. Assignment schemes commonly adopted areNearest Grid Point (NGP), Cloud In Cell (CIC), Trian-gular Shape Cloud (TSC) and Piecewise Cubic Spline(PCS), which corresponds to the zeroth, first, secondand third order piecewise polynomial functions W ( p ) with p = 0,1,2,3 respectively. For each of the schemes,we have a one-dimensional function: W (0) ( x ) = | x | < otherwise (12) W (1) ( x ) = − | x | | x | < otherwise (13) W (2) ( x ) = − x | x | < ( − | x | ) ≤ | x | < otherwise (14) W (3) ( x ) = (4 − x + 3 | x | ) 0 ≤ | x | < (2 − | x | ) ≤ | x | < otherwise (15)The mass assignment weight function is given by W ( x ) = d (cid:81) j W ( p ) ( x j /H ), where d is the dimension and H is the linear size of the grid cell.For a system with the number density of particles: n ( x ) = N p (cid:88) i =1 δ ( x − x i ) (16)where x i is the position of particle i , N P is the numberof particles and δ ( x ) is the Dirac delta function. Themass density after assignment can be described by: ρ ( x ) = mV g (cid:90) n ( x (cid:48) ) W ( x − x (cid:48) ) d x (cid:48) (17)where m is the particle mass and V g is the volume of thegrid. The mass density grid used for FFT can be gen-erated by a sampling process described by the samplingfunction: III( x ) = (cid:88) n δ ( x − n ) (18) ρ (cid:48) ( x ) = III( x H ) ρ ( x ) (19)with n an integer vector for denoting the coordinatesof the grid points. Thus, the whole mass assignmentprocess can be expressed as: ρ (cid:48) ( x ) = mV g III( x H ) W ( x ) ∗ n ( x ) (20)2.2.2. Potential solving, finite difference and forceinterpolation
The Green function G ( x ) describes the potential froma unit mass point source. The solution of Poisson’s equa-tion can be written as: ψ ( x ) = G ( x ) ∗ ρ (cid:48) ( x ) (21)The equation is usually calculated in Fourier space ( k ) inwhich the convolution is transformed to a simple mul-tiplication product. In the literature, there are manychoices for the Green function, e.g., the ”poor man’s”Green function ˆ G ( k ) ∝ /k . However, here we willtreat it as a free parameter which will be optimized af-ter we take into account of all the errors in the PMcalculation. Xu & Jing
The force field is calculated from the potential by E = − ∇ ψ . However, the differential is approximated by adifference in PM: E ( x ) = − D ( x ) ∗ ψ ( x ) (22)The first order approximation for the differential is thetwo-point finite-difference, D ( x ) = δ ( x + H ) − δ ( x − H )2 H (23) D ( x ) = d (cid:88) j =1 D ( x j ) ˆ x j (24)The error of the difference is at the order of O ( H ). Andthe next order approximation is the four-point finite-difference, D ( x ) = β δ ( x + H ) − δ ( x − H )2 H + (1 − β ) δ ( x + 2 H ) − δ ( x − H )4 H (25)Setting β = 4 / O ( H )term and the accuracy of order O ( H ) can be achieved.Higher order approximations can be more accurate, butthey are also more time-consuming in the calculationsince more points are used for the difference. The forcefield grid E (cid:48) calculated by PM is also a sampling of theunderlining force field E : E (cid:48) ( x ) = − III( x H ) D ( x ) ∗ ψ ( x ) . (26)Finally, the force field should be interpolated to thewhole space. We use the same interpolation function forthe force as we did for the mass assignment, F ( x ) = mV g W ( x ) ∗ E (cid:48) ( x ) . (27)When the two interpolation functions are taken thesame, the momentum of the system is conserved.2.2.3. Optimized Green function
Though the Green function ( ˆ G ( k ) ∝ /k ) is accuratein poor man’s Poisson Solver, errors are generated in ev-ery step of the PM calculation. For example, in the massassignment, error called alias is generated during sam-pling if ˆ ρ ( k ) is non-zero beyond the Nyquist frequency:ˆ ρ (cid:48) ( k ) = + ∞ (cid:88) n = −∞ ˆ ρ ( k − n k g ) (28)with k g = 2 π/H . Using a finite-difference to replacedifferential also generates error, and so does the inter-polation. Thus, the Green function is better regarded as a free parameter to be optimized to compensate forall the errors arising in the PM calculation.The Green function is optimized to minimize the av-erage force error between two particles. The force on aunit mass particle at x from another unit mass particleat x is given by: F ( x ) = (cid:90) d k (2 π ) d ˆ F ( k ) e i k · x (29)From equations 20-22 and 26-28, we get ˆ F ( k ) = 1 V g ˆ W ˆ E (cid:48) = − V g ˆ W ˆ D ˆ ψ = − V g ˆ W ˆ D ˆ G ˆ ρ (cid:48) = − V g ˆ W ˆ D ˆ G (cid:88) n ˆ W ( k − n k g ) e − i ( k − n k g ) · x (30)The last step uses the Fourier transform of the δ func-tion that δ ( k ) = e − i k · x . Thus, the force between twoparticles in a periodic system becomes: F ( x ; x ) = 1 V b (cid:88) l ˆ WV g ˆ D ˆ G (cid:88) n ˆ W ( k n ) V g e − i k n · x e − i k · x (31)Where l is the wavenumber of the Fourier series for theperiodic system, V b is the box volume and k n = k − n k g .The displacement-averaged total squared deviation ofthe PM calculated force F ( x = x − x ; x ) from thereference interparticle force R ( x ) is defined as: Q = 1 V g (cid:90) V g d x (cid:90) V b d x | F ( x ; x ) − R ( x ) | (32)Using Fourier transform, Q can be written as: Q = 1 V b (cid:88) l (cid:26) | ˆ D | ˆ G (cid:20) (cid:88) n ˆ U ( k n ) (cid:21) − G ˆ D · (cid:88) n ˆ R ∗ ( k n ) ˆ U + (cid:88) n | ˆ R | (cid:27) (33)where ˆ U = ˆ W /V g . Minimizing Q with respect to ˆ G givesthe optimized Green function:ˆ G ( k ) = ˆ D ( k ) · (cid:80) n ˆ U ( k n ) ˆ R ∗ ( k n ) | ˆ D ( k ) | (cid:20)(cid:80) n ˆ U ( k n ) (cid:21) (34) M lensing
PP force calculation
Due to the very nature of the PM calculation, thePM force is both softened and anisotropic on small scale(smaller than a few H ) compared to the reference forcebetween two point particles. The PP algorithm is usu-ally used to accurately calculate the short-range force, inorder to compensate for the PM softening. To suppressthe error from the anisotropy and to smoothly transitfrom the PM force to the PP force, a shape S ( r ) is as-signed to each particle for the PM, reflected in the soft-ened R : ˆ R = ˆ R p ˆ S (35)where ˆ R p is the Fourier transform of the force betweentwo point sources. The shapes with polynomial formsare usual choices for convenient analytical calculation,such as: S ( r ; a ) = /πa (cid:0) a / − r (cid:1) r < a/ otherwise (36) S ( r ; a ) = /πa (cid:0) a/ − r (cid:1) r < a/ otherwise (37)in the 2D case. The smoothing length is controlled by a = a pm , which also influences the volume to performPP. Larger a pm can better suppress the total error atthe cost of many more particle pairs in PP ( O ( a dpm )).The force used for PP is simply R pp = R tot − R pm ,where R tot and R pm are the reference total force andPM force respectively.
2D P M ALGORITHM FOR GRAVITATIONALLENSINGFirst of all, the analytical form of the optimized Greenfunction should be obtained. The analytic forms of ˆ D and ˆ U are provided in Hockney & Eastwood (1981):ˆ D x ( k ) = iβ sin ( k x H ) H + i (1 − β ) sin (2 k x H )2 H ˆ D y ( k ) = iβ sin ( k y H ) H + i (1 − β ) sin (2 k y H )2 H (38)ˆ U ( k ) = (cid:32) sin ( k x H ) sin ( k y H ) k x H k y H (cid:33) p +1 (39)where p is the order of the piecewise polynomial functionused for mass assignment and force interpolation.Obtaining ˆ R in 2D is more complicated than in 3D.We first calculate ˆ R p and then consider ˆ S for a selectedparticle shape. For unit mass, R p ( r ) = 2 G r /r , the 2D Fourier transform of the x component is:ˆ R px ( k ) = 2 G (cid:90) ∞∞ (cid:90) ∞∞ xx + y e − ( ik x x + ik y y ) dxdy = − πGi (cid:90) ∞∞ sgn ( k x ) e −| k x y |− ik y y dy = − πGi sgn ( k x ) | k x | k x k x + k y = − πG ik x k x + k y (40)Where sgn is the sign function. Similarly we have:ˆ R py ( k ) = − πG ik y k x + k y ˆ R p ( k ) = − πG i k k (41)Two forms of particle shape, Eqs. 36 and 37, are con-sidered and compared here. S is widely used in 3Dsimulations, which is however not preferred in the 2Dcase as will be shown. The Fourier transform of S is,ˆ S ( k ; a ) = (cid:90) π e − ikr cos θ dθ (cid:90) a πa ( a − r ) rdr , (42)where θ is the angle between r and k . Using Ja-cobi–Anger expansion, e iz cos θ = J ( z ) + 2 + ∞ (cid:88) n =1 i n J n ( z ) cos ( nθ ) , (43)where J n is the n -th order Bessel functions of the firstkind, we getˆ S ( k ; a ) = 64 a (cid:90) a ( a r − r ) J ( kr ) dr = 128 k a J ( ka − k a J ( ka S ,ˆ S ( k ; a ) = 12 k a [ J ( ka ( ka − J ( ka ( ka , (45)where H n is the n -th Struve function.The analytic form of ˆ G can be derived with equations38-39,41,44 and 45. The summation (cid:80) n can be donewithin | n | <
2, for very fast decay of ˆ U and ˆ R withrespect to | k | .Free parameters p , a pm and the shape form remainto be determined. The main improvement with higher Xu & Jing a pm / H e rr o r Q + ( S ) Z + ( S ) P + ( S ) Q + ( S ) Z + ( S ) P + ( S ) Figure 1.
The dependence of (cid:112) Q + , √ Z + and √ P + on a pm . The solid lines show the results for particle shape S and dashed ones for S . order mass assignment p is to suppress the anisotropy(Efstathiou et al. 1985). Thus, in the later part of thispaper, we set p = 2 by default with the TSC mass as-signment scheme. As we will show, with this scheme,we can achieve a sufficiently good isotropy with a rea-sonable choice of a pm .The error Q (equation 33) in PM with different choicesof particle shape and a pm can be calculated for the op-timized Green function ˆ G . To investigate the sources ofthe error, Q can be split into two parts, Z = (cid:90) d x |(cid:104) F ( x ; x ) (cid:105) − R ( x ) | = 1 V b (cid:88) l (cid:26) | ˆ D | ˆ G (cid:20) (cid:88) n ˆ U (cid:21) − G ˆ D · (cid:20) (cid:88) n ˆ R ∗ ˆ U (cid:21) + (cid:88) n | ˆ R | (cid:27) , (46)and P = 1 V g (cid:90) V g d x (cid:90) V b d x | F ( x ; x ) − (cid:104) F ( x ; x ) (cid:105)| = 1 V b (cid:88) l | ˆ D | ˆ G (cid:26)(cid:20) (cid:88) n ˆ U (cid:21) − (cid:88) n ˆ U (cid:27) . (47)In the above equations, Z is the total squared de-viation of the displacement-averaged mesh force fromthe reference force, which describes the deviation of themean force at separation r from the reference force, and P is the total squared deviation of the mesh force fromits displacement average, which reflects the anisotropyof the PM force. By definition, we have Q = Z + P . To construct the dimensionless error estimator, we define Q + = Q ( GH ) ( πH ) (48) Z + = Z ( GH ) ( πH ) (49) P + = P ( GH ) ( πH ) . (50)The dependence of (cid:112) Q + , √ Z + and √ P + on a pm isshown in Figure 1. Results for the two different particleshapes are also presented. In general, S performs bet-ter than S . And considering that the calculation of J n is much faster than H n , we prefer S in our 2D lensingcalculation. As expected, with increasing a pm , the erroris decreasing. At small a pm , the error is mainly fromthe deviation of the mean force from the reference force( Z ). Gradually, the anisotropy ( P ) becomes a domi-nant part of the error as a pm increases. The decreaseof the error is much slower once a pm reaches 6 H . Con-sidering the computational cost increases with a pm inthe PP calculation, we recommend a pm ≈ H . More-over, using higher order finite-difference approximationsmay also reduce the error, which we don’t investigatefurther here, as the PM force is accurate enough (betterthan 0.1 percent on all scales, cf. Fig.2) with the aboveparameters.To perform PP on the basis of PM, an analytic form of R pp ( r ), determined by R tot and R pm , is needed. R tot should also be smoothed to recover the underlining mat-ter distribution from the Monte Carlo sampling parti-cles. Hence, real space force R ( r ; a ) between two par-ticles of unit mass with shape S for any a should beknown. Actually, knowing the magnitude of the force R ( r ; a ) at distance r is enough since R doesn’t dependson the direction of r and the direction of R can be eas-ily calculated from the particle positions. The Fouriertransform of the x component of R ( r ; a ) is:ˆ R x ( k x , k y ) = − πG ik x ˆ S ( k ; a ) k . (51)Transforming it to real space, R x ( x, y ) = − πG (2 π ) (cid:90) ∞∞ (cid:90) ∞∞ ik x ˆ S ( k ) k e ik x x + ik y y dk x dk y , (52) M lensing R ( r ; a ) from R x ( x, y ) by setting x = r and y = 0, R ( r ; a ) = R x ( r, − Gπ (cid:90) ∞∞ (cid:90) ∞∞ ik x ˆ S ( k ) k e ik x r dk x dk y = − Gπ (cid:90) ∞ (cid:90) π i cos θ ˆ S ( k ) e ikr cos θ dkdθ . (53)Using Jacobi–Anger expansion, and integrating θ com-ponent, we get (cid:90) π cos θe ikr cos θ dθ = 2 + ∞ (cid:88) n =1 (cid:90) π i n J n ( kr ) cos ( nθ ) cos θdθ . (54)Only the n = 1 component survives in (cid:82) π cos ( nθ ) cos θdθ , thus (cid:90) π cos θe ikr cos θ dθ = 2 (cid:90) π iJ ( kr ) cos θdθ = 2 πiJ ( kr ) . (55)Finally we have R ( r ; a ) = 2 G (cid:90) ∞ J ( kr ) ˆ S ( k ; a ) dk . (56)The integration is hard to calculate analytically, so weprovide an accurate ( < − ) fitting formula for fastcalculation:2 G a (0 . ξ − . ξ − . ξ +2 . ξ − . ≤ ξ < a ( − . ξ + 4 . ξ − . ξ +15 . ξ − . . /ξ ) 1 ≤ ξ < /r < ξ (57)where ξ = 2 r/a . Thus, R pp ( r ) = R ( r ; a pp ) − R ( r ; a pm ) (58)The PP calculation for each point of interest should beperformed within the space range max ( a pm , a pp ) since R pp vanishes beyond, so the computational cost highlydepends on the the choice of these two parameters. Wehave discussed the choice of a pm already, and we willdetermine a pp in next section.The purpose of this paper is to accurately calculatethe 2D potential field E , which is proportional to thedeflection angle α . Although we only discussed aboutthe calculation of R tot , it can be converted to E bysimply dividing it by the particle mass. U H ( ff t h e o / f t h e o ) PM (1/ k )P M (optimized)10 U H f / f Figure 2.
Top: The mean deviation of the computed forcefrom the theoretical prediction at distance r from a pointsource. Bottom: Anisotropy, defined as the rms deviationof the force field ( σ f ) divided by the mean force field ( ¯ f ), atdistance r from a point source. The blue lines show the re-sults for our algorithm and the red ones for the conventionalPM algorithm. 4. TESTIn this section, we test the accuracy of our algorithmand determine the smoothing length for gravitationallensing studies in N-body simulations.4.1.
Point source
We first calculate the 2D force field from a point sourcewith our algorithm, and compare the performance withthat of the poor man’s
PM Poisson solver ( G ( k ) ∝ /k ). We generate particles with different positions in thesame grid cell, and calculate force field for different di-rections at a particular r . Errors of our algorithm andof the poor man’s PM are shown in Figure 2. The testis done with H = 1 kpc/h , a pm = 6 H , a pp = 0, particlemass m = 5 × M (cid:12) /h and box size 10000 kpc/h . Thelarge box is used to reduce the effect of the periodic im-ages. The mean deviation from the theoretical force fieldand the standard deviation in different directions of ouralgorithm are both two orders of magnitude smaller thanthe PM at r smaller than a pm , which can be attributedto our optimized PM and PP calculations. At large r ,the mean deviation of the two algorithms is nearly thesame but the standard deviation is slightly smaller in ourP M algorithm due to optimized ˆ G . Overall, the forcebetween the two particles is calculated with the errorsmaller than 10 − in our P M algorithm at all scales.
Xu & Jing
Halo with NFW profile
As shown in the previous subsection, we can achievea high precision in the force calculation. However, darkmatter distribution is represented by particles in N-bodysimulations, which inevitably leads to a Poisson noise inthe lensing quantities. To suppress the Poisson noise,one needs to smooth out the density distribution, whichcan be done in our algorithm by just choosing a propersoftening length a pp in the PP calculation. The soft-ening length needs to be large enough to smooth thedensity field but be small enough to maintain the un-derlining density distribution. We adopt an adaptivesoftening length a pp for each point of interest as the dis-tance to the N nb -th nearest neighbor particle. We cal-culate lensing quantities κ , γ and µ for NFW (Navarroet al. 1997) halos, and compare the results with the the-oretical prediction. From the comparison, we test theperformance for different N nb and try to give a recom-mendation for the choice of N nb .We test our algorithm with an NFW halo of the virialmass M vir = 10 . M (cid:12) /h and the virial radius R vir =709 kpc/h ( z = 0 . × M (cid:12) /h , which isthe typical mass resolution (10 ∼ M (cid:12) /h ) for cos-mological simulations. The calculation is done with H = 1 kpc/h , a pm = 6 H and box size 6 R vir . Red-shift is set to be 1.0 and 0.5 for the source and lensplanes respectively. In Figure 3, we show the rela-tive deviation of the mean values from the theoreticalpredictions and the standard deviation normalized bythe mean value for κ , γ and δµ = µ − r . The former (the mean deviation) shows the system-atic deviation caused by the smoothing, and the latter(the standard deviation) reflects the shot noise fluctua-tion left by the smoothing. As expected, the standarddeviation is decreasing with increasing N nb for all lens-ing parameters at nearly all radius. If only the Poissonnoise were considered, N nb should be as large as possi-ble. However, the mean quantities will deviate from thepredictions due to the smoothing. We take the adaptivesoftening length a pp at the center of the halo as a min which is also the smallest softening length for the halo.As shown in Figure 3, the mean deviation at the cen-tral region r < a min increases with the increase of N nb .Thus, a compromise for the choice of N nb can be reachedin order to make the mean deviation and the standarddeviation to be comparable. From the test, we find that N nb = 400 ∼
800 are recommended if the very inner part r < a min is not considered. Under this choice, therelative standard deviation is smaller than 7% for κ , 20%for γ and 8% for δµ , and the mean deviation is smallerthan 4% for κ , 12% for γ and 5% for δµ at r > a min .To apply our code to cosmological simulations, perfor-mance on smaller halos should also be investigated. InFigure 4, we show the results for halos with mass from10 . M (cid:12) /h to 10 . M (cid:12) /h . The concentration parame-ter is taken again from Dutton & Macci`o (2014). N nb isset to 400. For the smaller halos, our code still performswell for their outer parts ( r > a min ), and our resultsdoes not depend on the halo mass.Let us consider a halo simulated with different res-olutions, and investigate the choice of N nb . We stilluse the NFW halo of the virial mass M vir = 10 . M/h .The higher resolution one has a particle mass 10 M (cid:12) /h ,and the lower resolution one has the same resolution asbefore. The performance with different N nb is shownin Figure 5. The standard deviation can be reducedroughly as N − / nb for all the three quantities at the ra-dius r > a min . It is also interesting to note the standarddeviation remains the same for the different resolutionswhen N nb is taken the same. In contrast, the systematicdeviations are only slightly reduced the higher resolutionsimulations even with the increase of N nb . Because theshear γ is a long-range interaction quantity, the smooth-ing at the central region of a halo can result in a signifi-cant systematic deviation at r smaller than a few a min ,which in turn leads to the error of the magnification.From the above tests, we can conclude that adaptivesmoothing with N nb ≈
400 should work for most lens-ing studies. If one needs a smaller statistical deviation(noise), one may increase N nb but at the expense of los-ing resolution at high density regions. The halo coresof r < a min cannot be resolved, and the shear in theinner regions could be significantly underestimated bythe smoothing. P3MLENS : A PYTHON IMPLEMENTATION OFTHE 2D P M ALGORITHMWe provide a python script
P3MLens that imple-ments the 2D P M algorithm for gravitational lensing.
P3MLens constructs a lens plane from input particles andreturn the quantities like E , α , κ , γ and µ for positionsof interest. We have taken full use of the numerical pack-ages numpy , scipy and numba , so the computation isdone in an efficient way in python . https://numpy.org http://numba.pydata.org M lensing | t h e o |/ t h e o | t h e o |/ t h e o U U s | t h e o |/ t h e o / / N=100N=200N=400N=800 U U s / Figure 3.
The relative deviation of the mean value from the theoretical prediction (left) and the standard deviation normalizedby the mean value (right) for κ (top), γ (middle) and δµ (bottom) as function of radius. The results for different smoothinglengths, denoted by N nb , are shown in different colors. For each N nb , a vertical dashed line is drawn to indicate the smallestsmoothing length a min around the halo center. | t h e o |/ t h e o | t h e o |/ t h e o M / h M / h M / h U U s | t h e o |/ t h e o / / U U s / Figure 4.
The same as Figure 3 but for halos with a different mass. N nb is set to 400 Xu & Jing | t h e o |/ t h e o | t h e o |/ t h e o U U s | t h e o |/ t h e o / / M / h , 40010 M / h , 40010 M / h , 400 5010 M / h , 20000 U U s / Figure 5.
Comparison of the performance for simulations with different mass resolution for a NFW halo of the virial mass M vir = 10 . M (cid:12) /h . The low resolution halo has a particle mass 5 × M (cid:12) /h , and the high resolution one has a particle mass10 M (cid:12) /h . The results are presented for the high resolution halo with different N nb = 400 , √ , N nb = 400.6. CONCLUSIONIn the paper, we introduce a 2D P M algorithm withoptimized Green function and adaptive softening lengthfor gravitational lensing studies in cosmological simula-tions. The algorithm yields a precise calculation for the2D force field between particles, which is two orders ofmagnitude more accurate than a simple PM algorithmat a small r , and has a smaller anisotropy at large r .Comparing the total error ( Q ), mean deviation ( Z ) andanisotropy ( P ) for different smoothing schemes in thePM, we prefer a pm ≈ H and particle shape S . Fromour tests of computing κ , γ and δµ for a halo of the NFWprofile, we found an adaptive smoothing with N nb = 400in PP can largely suppress the Poisson noise in the cos-mological simulation, though the high density cores of a mass corresponding to N nb are smoothed out. Onemay reduce the Poisson noise further by increasing N nb but at the expenses of losing resolution at high densityregions.We will release the python implementation P3MLens of this algorithm once the paper is published.ACKNOWLEDGMENTSThe work is supported by NSFC (11890691, 11621303,11533006) and by 111 project No. B20019. Xu thanksWei Tian for useful discussion about analytical calcula-tion of integration involving Bessel functions. This workmade use of the Gravity Supercomputer at the Depart-ment of Astronomy, Shanghai Jiao Tong University.REFERENCES
Aihara, H., Arimoto, N., Armstrong, R., et al. 2018, PASJ,70, S4, doi: 10.1093/pasj/psx066Bartelmann, M., & Weiss, A. 1994, A&A, 287, 1.https://arxiv.org/abs/astro-ph/9311074Bradaˇc, M., Schneider, P., Lombardi, M., et al. 2004, A&A,423, 797, doi: 10.1051/0004-6361:20040168Couchman, H. M. P. 1991, ApJL, 368, L23,doi: 10.1086/185939 Dey, A., Schlegel, D. J., Lang, D., et al. 2019, AJ, 157, 168,doi: 10.3847/1538-3881/ab089dDutton, A. A., & Macci`o, A. V. 2014, MNRAS, 441, 3359,doi: 10.1093/mnras/stu742Efstathiou, G., Davis, M., White, S. D. M., & Frenk, C. S.1985, ApJS, 57, 241, doi: 10.1086/191003Fosalba, P., Gazta˜naga, E., Castander, F. J., & Manera, M.2008, MNRAS, 391, 435,doi: 10.1111/j.1365-2966.2008.13910.x M lensing11