Rapid Computation of the Plasma Dispersion Function: Rational and Multi-pole Approximation, and Improved Accuracy

The plasma dispersion function $Z(s)$ is a fundamental complex special integral function widely used in the field of plasma physics. The simplest and most rapid, yet accurate, approach to calculating it is through rational or equivalent multi-pole expansions. In this work, we summarize the numerical coefficients that are practically useful to the community. Besides the Pade approximation to obtain coefficients, which are accurate for both small and large arguments, we also employ optimization methods to enhance the accuracy of the approximation for the intermediate range. The best coefficients provided here for calculating $Z(s)$ can deliver twelve significant decimal digits. This work serves as a foundational database for the community for further applications.

The plasma dispersion function Z(s) is a fundamental complex special integral function widely used in the field of plasma physics.The simplest and most rapid, yet accurate, approach to calculating it is through rational or equivalent multi-pole expansions.In this work, we summarize the numerical coefficients that are practically useful to the community.Besides the Padé approximation to obtain coefficients, which are accurate for both small and large arguments, we also employ optimization methods to enhance the accuracy of the approximation for the intermediate range.The best coefficients provided here for calculating Z(s) can deliver twelve significant decimal digits.This work serves as a foundational database for the community for further applications.

I. INTRODUCTION AND MOTIVATION
The plasma dispersion function [1][2][3] is defined as with s = x + iy, which is valid for y > 0. For y ≤ 0, the function is analytically continued from the above upper plane form to the lower plane.The function is relevant to the Faddeev function w(s) and error function erf(s) by = i √ πe −s 2 − 2F (s), (5) where i = √ −1, erf(s) = 2 √ π s 0 e −t 2 dt and erfi(s) = −i • erf(is).Here, F (x) is Dawson integral F (x) = e −x 2 x 0 e t 2 dt = √ π 2 e −x 2 • erfi(x).Note also dZ/ds = −2[1 + sZ(s)].The below symmetry properties (the asterisk denotes complex conjugation) hold And the two-side Taylor expansion approximation and Γ is Euler's Gamma function.A further expansion is 8) would not match well of the imag part for the range y < √ πx 2 e −x 2 when x ≫ 1.Since e −x 2 < 10 −40 for x > 10, the term actually mainly affects the intermediate range, e.g., 1 < x < 10.
It is also surprising that the rational and multi-pole approximation of the plasma dispersion function inspire two unexpected applications: (1) The development of Landau fluid models to mimic the kinetic Landau damping effects [18][19][20]; (2) The first solver to obtain all the kinetic dispersion relation solutions without the requirement of an initial guess [15,16].These two applications are possible only when the approximation of Z(s) keeps the following features: (a) The approximation should be rational functions; (b) One formula for all the interesting regions; (c) To maintain the same accuracy with as few terms as possible.
Hence, segmentation calculations (cf., [8]) and nonrational expansion are not our choices.We find our only choice is the rational and multi-pole approximation.The standard approach to obtaining the rational coefficients is to use Pade approximation to match Eq.( 8), which is accurate for small (s → 0) and large (s → ∞) arguments.The coefficients can be calculated rigorously via matrix inverse.The Pade approximation is less accurate at the intermediate range s ≃ 1.We then use optimization methods to reduce the error of the approximation at the intermediate range, which can improve one order of magnitude.Weideman's [9] method is also in rational form; however, it will require series of high-order terms.
We found rational and multi-pole coefficients are not provided systematically in literature.The methods used in this work are probably not new, and some of the coefficients have been provided in different literature, for example, the Pade rational form with also analytic coefficients up to J = 2 − 8 [4,18], multi-pole for small J ≤ 8 [2,4,15,17] and extended to J = 24 [16].The coefficients with optimization for J = 5, 6, 7 can also be found in rapid calculation of Z, such as Refs.[5][6][7].A systematic summary of the rational coefficients and corresponding multi-pole coefficients could be useful to the community, especially for beginners.
The purpose of the present work is to provide comprehensive numerical Pade coefficients from small order to high accuracy, typically J = 2 to J = 24, which can have an error less than 10 −13 , and could be used for rapid numerical calculation of Z(s) to high accuracy.The coefficients of improved fitting for small J ≤ 8 are also provided, which can be used to save computation time or improve accuracy with fewer terms in Landau fluid model and kinetic dispersion relation solver.
In Section II, we describe the approaches we used and the results we obtained.In Section III, summary and conclusion are given.

II. APPROACHES AND RESULTS
The rational approximation and multi-pole expansion of Z(s) is (typos in Ref. [15] are fixed here) with q 0 = 1.This form is valid for the upper plane to high accuracy and is analytical (i.e., automatically be analytically continued) and thus would also have a good approximation at the real axis and lower plane in case the y is not far from the real axis.If one needs a more accurate value in the lower plane, Eq.( 7) can be used to use the value from the upper plane and symmetry property.Hence, the entire plane can be calculated to high accuracy.

A. Pade method
The Pade expansion is to match Eq. (10) to Eq.( 8), i.e., with terms to terms match To obtain the coefficients, the system of equations to be solved are where I + K = 2J, and p j = 0 for j > J − 1 and j < 0, and q j = 0 for j > J and j < 0. The 2J equations determine 2J coefficients p j and q j .Here I means keeping I equations for s → 0, and K means keeping K equations for s → ∞.The above equation can be solved via matrix inverse, both analytically for small J or numerically for arbitrary J.
The term iσ √ πe −s 2 for s → ∞ in Eq.( 8) is omitted, and can be explicitly rewritten as And Eq.( 12) is rewritten as For a given J and I, the coefficients p l and q k can be solved easily via matrix inverse.Solving for the multipole coefficients b j and c j is also straightforward (e.g., using the residue() function in Matlab), i.e., c j are the roots of the equation q 0 + J k=1 q k s k = 0.There are also symmetric features to ensure that b j and c j occur in pairs: b j = b * J+1−j and c j = −c * J+1−j .For multi-pole expansion, we have Comparing Eq.( 15) with Eq.( 13), we have √ π, and j b j = −1, j b j c j = 0, and j b j c 2 j = −1/2.For kinetic dispersion relation solver [15,16], j b j c 2 j = −1/2 is used, which means that we need to keep to O(1/s 3 ) when calculating the multi-pole coefficient.Hence, we should have K ≥ 3.This also implies that −q J = p J−1 , −q J−1 = p J−2 , and −q J−2 − 1 2 q J = p J−3 hold, and the R(s) = 1 + sZ(s) coefficients in the Landau fluid model can be obtained straightforwardly.For small J = 2, 3, 4, we may only keep to O(1/s) or O(1/s 2 ).The coefficients taken by Ronnmark [17] are J = 8, I = 10.
It also appears that a slightly larger I [4] than K can provide a better overall approximation.The notation Z IK is used to describe different orders of approximation.

B. Search minimum method
It appears that p 2j and q 2j+1 are pure imaginary numbers, and p 2j+1 and q 2j are pure real numbers, for j = 0, 1, 2, • • • .
We minimize both the absolute and relative errors by with and Z A (s) is the approximation value, and Z(s) is the accurate value.In this work, we set w = 0.5.We take the accurate Z(s) to be the one using the Dawson function in Matlab (though it may not be accurate for all ranges).We use some of the constraint Eq.( 14) to make small s to O(s 2 ) and large s to O(1/s 3 ).There are some standard approaches to obtain the optimized p l and q k .Here, we use the Matlab function fminsearch() to perform the calculations and use the Pade coefficients as initial guesses.The results could be sensitive to the initial guess and the iterative convergence criteria.Hence, the final results may not be determined.We choose our best obtained results to list here.The optimization is performed for s = x + iy, with x = [−50, 50] and y = 0.1.

C. Coefficients tables
After the methods described in the subsections above, we performed calculations to obtain the coefficients p l , q k , b j , c j , and corresponding errors δ a and δ r in a table.For Pade coefficients, we calculated from J = 2 to 24, and from J = 2 to 8 for optimized coefficients.The error data are taken for the lower plane with y = −0.1 and x = [−50, 50], ensuring that the approximation can also accurately capture the Landau damping effect even without the 2i √ πe −s 2 term.Data for J > 24 are not listed here, as the double precision is not adequate and not necessary for most applications.
Figures 4, 5, and 6 show the coefficients, while Figures 1, 2, and 3 show the results of the errors.'Pade best' in the figures means that I is chosen such that the error is minimized according to the tables.
It is possible to optimize for J ≥ 9, but it is less useful, and thus we do not provide it here.The reason is that J = 8 is sufficient for most purposes, and those requiring higher accuracy can use slightly higher J Pade coefficients; for example, the J = 10 Pade may have better accuracy than the optimized J = 8.For J > 20, random error becomes the major issue for the approximations, and reducing the error to less than 10 −13 becomes difficult.Therefore, for double precision usage, J = 20 is sufficient.The J = 24 coefficients listed here are for reference, for those who hope to develop more accurate approximations.Due to the random off error, calculating using p l and q l would yield better accuracy for large z than using b j and c j .
It is claimed [8] that if the term iσ √ πe −s 2 is omitted, Eq.( 8) would not match well for the range y < √ πx 2 e −x 2 when x ≫ 1.However, in practical tests, this only slightly affects the imag part of intermediate x, and the Weideman method also holds well for y ≃ 0.
In Appendix B, we provide Matlab and Python code examples of how to use these coefficients to calculate Z(s).In Appendix A, we also provide the coefficients using the Weideman method.

D. Validation
The accuracy of Z is not just about the function itself, but rather whether it can accurately capture the physical phenomenon.We employ the one-dimensional electrostatic Landau damping problem to demonstrate the accuracy of the data provided in the table.The problem can be simplified as follows [21] with z = ω/( √ 2k).For a given k, we aim to solve for ω, focusing on the least damping branch.
Results are depicted in Fig. 7 and Fig. 8.The J = 2 method cannot be used, as the large error leads to incorrect roots.J = 3, 4 could be viable options for low accuracy calculations.We observed that for the Pade method, larger J generally results in better accuracy across all arguments.However, for the same J, optimized coefficients may only control total errors, with local errors potentially increasing, especially for small s.Thus, optimization coefficients may not always outperform unoptimized coefficients.This suggests that if only accurate Z(s) is needed (i.e., not aiming to reduce the expansion order), using coefficients from the rational expansion with higher orders (i.e., larger J) is preferable over optimizing the coefficients.Hence, we recommend using J ≥ 8 from our table, instead of the older lower J data [5][6][7].The coefficients provided here for calculating Z(s) can deliver up to   twelve significant decimal digits, limited by double precision data rather than the approach itself.Therefore, the Pade expansion method can compete with Weideman's method and can even be simpler.
Breakdown region is for k ≤ 0.15, where y < 10 −7 ≪ 1, incorrect positive numerical solutions of ω i may occur.However, since ω i /ω r ≪ 10 −8 is small, it is less critical to achieve high accuracy for most applications.Humlicek [7] reported that any rational approximation suffers inevitable failure near the real axis.For methods to address this issue, one can refer to [7,8].

III. SUMMARY AND CONCLUSION
We revisit the issue of rapid calculation of the plasma dispersion function and provide comprehensive rational and multi-pole coefficients for reference.A practical application of this work is to accelerate the computation of the PDRK/BO [15,16] code.For instance, an optimized J = 6 can achieve the same maximum error as the former Pade J = 8 method, with a complexity of O(J 2.7 ), resulting in a speedup of approximately 2.1 times.The optimized J = 8 method also reduces the maximum er-