Numerically stable form factor of any polygon and polyhedron

Coordinate-free expressions for the form factors of arbitrary polygons and polyhedra are derived using the divergence theorem and Stokes’s theorem. Series expansions are used to ensure numeric precision close to apparent singularities.


Overview
The term 'form factor' has different meanings in science and in engineering. Here, we are concerned with the form factor of a geometric figure as defined in the physical sciences, namely the Fourier transform of the figure's indicator function, also called the 'shape transform' of the figure.
This form factor has important applications in the emission, detection and scattering of radiation. Two-dimensional shape transforms are used in the theory of reflector antennas (Lee & Mittra, 1983). The three-dimensional form factors of the sphere and the cylinder go back to Lord Rayleigh (1881). Shapes of three-dimensional nanoparticles are investigated by neutron and X-ray small-angle scattering (Hammouda, 2010). Particles grown on a substrate (Henry, 2005) develop many different shapes, especially polyhedral ones, as observed by grazing-incidence neutron and X-ray small-angle scattering (GISAS, GISANS, GISAXS) (Renaud et al., 2009). Large collections of particle shape transforms have therefore been derived for and implemented in GISAS software (Lazzari, 2006;Pospelov et al., 2020); another GISAS software package uses surface triangulation for computing approximative form factors (Chourou et al., 2013). For attempts at direct reconstruction of polyhedral shapes from scattering patterns see the article by Engel & Laasch (2020) and literature cited therein.
In this paper, we derive a numerically stable algorithm for computing the form factor of any polygon or polyhedron, as implemented in the GISAS software BornAgain (Pospelov et al., 2020). Originally, this algorithm was documented in a terse mathematical note (Wuttke, 2017). In the present paper, derivations and results have been simplified, the material has been completely reorganized for better readability, and additional literature is taken into account.

Different ways to compute form factors
The form factor of a three-dimensional solid body Å & R 3 is Fðq; ÅÞ :¼ RR Å R d 3 r expðiqrÞ: ð1Þ ISSN 1600-5767 In most applications, the wavevectors q are real. In GISAS, however, the incident and scattered radiation may undergo substantial absorption, which can be modeled by an imaginary part of q. Therefore, we admit complex wavevectors q 2 C 3 . For any polyhedron, (1) can be evaluated analytically by successive integration in the three coordinates. This is straightforward for a cuboid with edges along the coordinate axes. In most other cases, the algebra is cumbersome, and the resulting expressions are complicated and unattractive in that they do not reflect symmetries of the underlying shape. Striking examples are the form factors of the Platonic solids worked out in a tour de force by Li et al. (2011).
It is therefore preferable to derive a coordinate-free solution of (1) that expresses the form factor of a generic polyhedron in terms of its topology and vertex coordinates. This has been undertaken in different ways by Senesi & Lee (2015), Croset (2017) and Wuttke (2017). Senesi & Lee (2015) decomposed the polyhedron into pyramids and wrote the polyhedral form factor as the sum of the pyramidal form factor evaluated at different rotated q. Croset (2017) decomposed the polyhedron into simplexes, as explored previously by Lien & Kajiya (1984) and most recently by Li & Xie (2020), for the integration of multinomials. Following Wuttke (2017), we here present a different derivation that is based on use of the divergence theorem and Stokes's theorem to reduce the volume integral to integrals over polygonal faces, and further reduce these surface integrals to line integrals over straight edges. This approach has been previously used for the integration of polynomials (Cattani & Paoluzzi, 1990;Bernadini, 1991) and for the computation of inertia moments (Mirtich, 1996).
Applications to nanoparticle assemblies typically require some averaging over particle sizes or/and orientations. How to compute these averages efficiently and with sufficient accuracy is an interesting and important question, which however is beyond the scope of the present work.

Singularities and asymptotes
All analytical expressions for polyhedral form factors, derived by whatever method, contain denominators that vanish at q = 0. Croset (2017) suggested, and we will confirm, that the degree of this singularity is closely related to the asymptotic envelope of F(q, Å) for large q, which goes as q À1 , q À2 or q À3 depending on whether q is perpendicular to a face or an edge or points in an off-symmetric direction.
However, there is nothing fundamental about the singularities at q = 0: From the definition (1) in conjunction with Leibniz's integral rule of differentiation we see that F(q, Å) is infinitely many times differentiable for all q 2 C 3 ; therefore F is a holomorphic function of each of the Cartesian components of q; therefore any apparent singularity is removable. Croset (2018) rederived the asymptotic envelopes by classifying the endpoint singularities of the section normal to q as function of height. In Section 3.6, we obtain them directly from our form factors.
The main purpose of this paper is to overcome numeric instabilities for small q and q k . The latter is the wavevector component in the plane of a polygonal face. We will explain how rounding errors can grossly distort form factors when q or q k is of the order of /a, where is the machine precision and a is a typical edge length.
At this point the reader may wonder whether wavevectors with extremely small, but nonzero, q or q k have any practical importance. If wavevectors were drawn at random from an entire Brillouin zone, then the chance of ever hitting numerically problematic values would indeed be negligible. Often, however, q is chosen along a face normal. Roundoff errors then easily yield a tiny nonzero q k , which causes huge, and symmetry breaking, errors in the form factor. Actually, this entire study started from the unexpected discovery of such artifacts in conventionally computed form factors.

Notation
A flat polygon À, embedded in three-dimensional space, shall be specified by its J vertices V j (j = 1, . . . , J). Vertex indices shall be understood modulo J so that V 0 V J . With this convention, the vertex sequence forms the closed loop @À. Edge j of the polygon is a straight line from V jÀ1 to V j . In most of this work it is more advantageously specified through its position and mid-to-end vector The normal vectorn n of the plane spanned by V j shall be oriented such that @À has the winding number +1 (fulfills the right-hand rule with respect ton n). The oriented plane characterized byn n induces a decomposition of any vector v 2 C 3 into a component perpendicular to the plane, v ? :¼ ðvn nÞn n; and an in-plane component, This decomposition will be applied to position vectors r and to wavevectors q. The oriented plane is fully specified by its normal vectorn n and its distance from the origin, r ? . Complex conjugation is denoted by a superscript asterisk. The absolute value of a complex vector is written q : Note that the in-plane unit vector differs from the in-plane componentq q k of the unit vectorq q. In this work, we shall only use b q k q k andq q, notq q k . The triple product is denoted research papers with the standard operators dot (Á) for the scalar product and cross (Â) for the vector product. Between adjacent vector symbols, as in the parentheses in (4), we omit the dot. The cardinal sine function sincðzÞ :¼ sinðzÞ=z has the analytic continuation sincð0Þ ¼ 1. The numeric implementation for |z| ! 0 is unproblematic: as sinðzÞ has full floatingpoint accuracy, so has sinðzÞ=z.

Form factor
We define the form factor of a flat figure À, embedded in three-dimensional space, as Proposition. The form factor (8) of a flat J-gon À is for q k 6 ¼ 0, with notations from Section 2.1, and with an arbitrary constant c that can be chosen for computational convenience. The value at q = 0 is the area of À, Values at q 6 ¼ 0, q k = 0 can be obtained from Proof. For any vector field G, we have Stokes's theorem: To prove (9), we choose GðrÞ :¼n n Â q Ã ½expðiqrÞ À c. The lefthand side of (12) is The right-hand side of (12) is Each edge can be written as a parametric curve r j ðÞ : With q k 6 ¼ 0, we obtain (9). Equation (11) follows directly from (8) and the fact that À is a flat figure with constant r ? . To prove (10), we use Stokes's theorem (12) with G ¼n n Â r=2. &

Remarks and example
A closed expression for the form factor of the polygon has long since been known (Lee & Mittra, 1983, equation 6). A more symmetric expression was obtained by Croset (2017) (equation 4, where i þ 1 and i À 1 should be swapped). In our notation, it reads The equivalence with our equation (9) is proven in Appendix A. Equation (15) is esthetically more pleasing than (9), but (15) is problematic for computer implementation and ill suited for the theoretical study of singularities, because for each j there are two q planes for which the denominator vanishes. Equation (10) is well known as the 'surveyor's formula'. The standard proof uses triangular tessellation (Braden, 1986). See Appendix B for a derivation of (10) from the q k ! 0 limit of (9). Fig. 1 shows for an equilateral triangle how strongly the form factor, plotted as | f(q)| versus q, varies with the wavevector directionq q.

Removable singularity
The closed expression (9) for the polygonal form factor f(q, À) has a singular factor q À1 k . As discussed in Section 1.3, the definition of f guarantees its analyticity. So we know that the apparent singularity at q k ! 0 is removable. We also know the value at q k = 0, given by (10) and (11). Nonetheless, the presence of a divergent factor may cause numeric instabilities Modulus of the form factor of an equilateral triangle as function of wavenumber q for three different wavevector directionsq q. The triangle lies in the xy plane and has a symmetry axis along x. The center of gravity is at the origin. The edge length is L = 1.
for small values of q k . To investigate this more closely, let us write (9) as with the function The constant c can be chosen differently for different q. At large q k , c ¼ 0 prevents roundoff errors in the numerator of (17). For small q k , c ¼ 1 is the better choice. To see this, we expand j as The leading, j-independent term in the expansion, with the apparent singularity (iq k ) À1 , contributes nothing to the sum in (16), whatever the value of c, because P j E j ¼ 0. This, however, holds only in exact arithmetics; in a floating-point implementation, roundoff errors can make the sum nonzero. For q k ! 0, any such error outgrows all other terms in the expansion (Fig. 2). Therefore, in the small-q k regime, the only sensible choice is c = 1, which lets the leading term vanish.
Unfortunately though, the subtraction of c ¼ 1 can cause roundoff errors in the numerator of (17). As a straightforward remedy, we compute j for small q k from its series expansion As a consistency check, we note that the limit j ðq ! 0Þ, at constant in-plane directionq q, is the n ¼ 1 termq qR j . Plugging this value into (16), we recover the surveyor's formula (10). The algebra is quite lengthy and therefore relegated to Appendix B. Fig. 3 compares the series expansion with the closed expression (9). The series expansion works well even beyond the first minima in f(|q|). In practice, the series expansion is only needed for qL ( 1, and therefore only a few expansion orders are needed to keep errors close to machine precision.

Polygon with inversion center
Computations can be simplified, and the numeric instability at q k ! 0 avoided, if a polygon has a perpendicular twofold symmetry axis (Schoenflies group S 2 ) or, equivalently, an inversion center at q. The form factor has the symmetry f (q ? + q k , À) = f (q ? À q k , 2q À À). As the number of vertices is even, we can write J = 2J 0 and use R jþJ 0 À q ¼ ÀðR j À qÞ and E jþJ 0 ¼ ÀE j . For q k 6 ¼ 0, the form factor is In contrast to (9), the summand in (20) has no constant contribution but is linear in q. There is no cancellation for q k ! 0 and no need to use a series expansion for the accurate computation of f.

Notation and parameterization
An orientable polyhedron Å shall be specified by its K polygonal faces À k (k = 1, . . . , K). For each face À k , the normaln n k :¼n nðÀ k Þ shall point to the outside of Å; this then determines the order of the vertices in the sequence V k1 ; . . . ; V kJ k .
In a computer implementation, the topology and geometry of a polyhedron can be specified through two arrays: Array C holds one coordinate triple V for each of the polyhedron's vertices. Array T holds one array k for each of the polyhedron's faces À k ; k holds the global indices kj of the vertices that belong to face À k , such that V kj ¼ V kj . In short, array C holds the coordinates and array T holds the topology of the polyhedron. For the latter, Schlegel diagrams (Fig. 4) provide a helpful visualization. In physical simulations, C is typically generated by a parametric function, whereas T is static. An Same form factor modulus as in Fig. 1, forq q ¼ ð0:6; 0:8; 0Þ, now in a double logarithmic plot, computed in double-precision arithmetic according to the closed expression (9) with c = 0. For some small q, the results are totally wrong owing to imperfect cancellation in the leadingorder term that ought to vanish.

Figure 3
Modulus of the form factor of the equilateral triangle of Fig. 1, as function of wavenumber q for wavevector directionq q ¼ ð0; 1; 0Þ. The black chain is computed using the analytical expression (9). The colored curves are computed using the series expansion (19) up to the indicated order n.
assertion in the computer code should ensure that all faces are planar for any geometry parameters.
Additionally, it is advantageous to foresee boolean parameters to indicate the presence or absence of inversion centers. One needs one such parameter for the entire polyhedron and one for each of its polygonal faces.

Form factor
Proposition. If q 6 ¼ 0, then the form factor (1) of a K-hedron Å is with an arbitrary constant C that can be chosen for computational convenience. Otherwise, for q = 0, the form factor is just the volume of Å, Proof. For a polyhedron, the divergence theorem takes the form RR With the choice H :¼ q Ã ½expðiqrÞ À C, this yields iq 2 Fðq; ÅÞ ¼ P K k¼1 q Ãn n k RR À k d 2 r expðiqrÞ À C ½ : With the notation (8), this proves (21). With the choice H :¼ r=3, we obtain the volume formula (22). & Similar to c in Section 2, the constant C can and should be chosen differently for different q domains. At large q, the best choice is C = 0. The small-q case is discussed in Section 3.3.
To see the equivalence of (12) with equation 15 of Croset (2017), we let C = 0, take f(q, À) from (15) and use the fact that E jÀ1 Â E j is colinear withn n.

Removable singularity
The closed expression (21) for the polyhedral form factor F(q, Å) contains two removable singularities: the explicit factor q À1 , and the factor q À1 k contained in the polygonal form factors f(q, À k ). For the case that only q k , but not q, is close to 0, we rely on the numerically stable computation of f derived in Section 2.4. Here we address the case q ! 0.
In analogy to Section 2.4, it is sufficient to invoke analyticity to convince ourselves that the singularity of F is removable. The value of F in the limit q ! 0 is just the volume of Å. The expansion of (21) starts with The leading, apparently singular term is identically zero because P kn n k ArðÀ k Þ ¼ 0. This, however, holds only in exact arithmetics; in a floating-point implementation, roundoff errors can make the sum nonzero. For q ! 0, any such error outgrows all other terms (Fig. 5). Therefore, in the small-q regime, the only sensible choice is C = 1, which lets the leading term vanish.
Unfortunately though, this can lead to roundoff errors in the difference f(q) À f(0) in (21). As a remedy, we compute the form factor from a series expansion as follows: Combine (21) and (16) to write the form factor as Modulus of the form factor of a truncated tetrahedron (trigonal pyramidal frustum). The base is an equilateral triangle in the xy plane, oriented so that an edge points in the y direction, with edge length L = 1; the dihedral angle is 72 ; the height H = L/2. The plot shows ReF versus q for the off-symmetric directionq q ¼ ð1; 2; 3Þ=14 1=2 . Blue spots are computed using the analytical expressions (9) and (21). For qL < 10 À6 , roundoff errors dominate. The orange line is computed according to (26) with summation (28) up to n = 19.

Figure 4
Schlegel diagram of a facetted cube. This representation of polyhedral topology can help to assign vertex indices in a systematic way, which then facilitates the coding of C and T . The topology array T has elements (0, 1, 2, 3, 4, 5, 6, 7), (0, 8, 1) etc. The coordinate array C, parameterized on lengths a and b, has elements (a, b, a), (b, a, a), (Àb, a, a) etc.
The constant (0, 1) neutralizes the first term in the expansion (19) of (q, 1) so that Fig . 5 shows that there is good overlap between the domains of the closed expression and the series expansion.

Polyhedron with inversion center
If a polyhedron has an inversion center at q (Schoenflies group C i ), then the form factor has the symmetry F(q, Å) = F(Àq, 2q À Å). As the number of faces is even, we can write K = 2K 0 . We require that faces numbered k and k + K 0 be opposite to each other. We usen n kþK 0 ¼ Àn n k to write the form factor as Fðq; ÅÞ ¼ expðiq qÞ 1 iq wheref f ðq; ÀÞ :¼ f ðq; ÀÞ À f ðq; ÀÀÞ ð 30Þ is the form factor of a pair of opposite faces. The symmetry f(q, ÀÀ) = f(Àq, À) allows some economy in computing F from the generic closed expression.
In the small-q case, the expansion (26) is symmetrized as and in consequence in (28) the terms with odd n cancel.

Prism
For a prism Å = {r | r k 2 À k , |r ? | < h/2} a much simpler solution is available. We return to the definition (1). Applying Fubini's theorem to factorize the triple integral from the onset into an integral (8) over the base À k of the prism and an integral along the normal directionn n, we obtain the form factor for all q. Thanks to the sinc function in (32), there is no singularity in q ? and therefore no series expansion is needed for q ? ! 0.

Asymptotic envelopes
We now come back to the asymptotic power-law envelopes for large q discussed in Section 1.3. A cube Å with side lengths L, centered at the origin and oriented along the coordinate axes, has the form factor Fðq; ÅÞ ¼ L 3 sincðq x L=2Þ sincðq y L=2Þ sincðq z L=2Þ: For large q it has the asymptotic envelope |F | 8/|q x q y q z |, which goes as q À3 for fixed directionq q, provided none of the three componentsq q x ,q q y ,q q z is zero. Ifq q is perpendicular to one of the edges of the cube, then one of the three sinc functions has the fixed argument 0 and value 1. And ifq q is perpendicular to one of the faces of the cube, then (33) has two constant factors sincð0Þ ¼ 1. As Croset (2017) has worked out, these observations can be generalized to any polygon. Within our present formalism, this can be confirmed as follows.
For q 6 ¼ 0, the form factor (21) of any K-hedron Å is limited by For q k 6 ¼ 0, the form factor (9) of any J-gon À is limited by For qE j 6 ¼ 0, the sinc function in (35) is limited by So if all the above conditions are fulfilled, then the polyhedron form factor F has an asymptotic envelope $q À3 . If there is any edge perpendicular to q, then (36) is not applicable, and max jsincðqE j Þj takes the q-independent value 1, so that the envelope of F goes with q À2 . If there is any face perpendicular to q, then (35) is not applicable, and max j f j in (34) takes the q-independent value ArðÀÞ; so the envelope of F goes with q À1 .

Implementation
Code for computing the form factor of any polygon or polyhedron, based on all the above, has been implemented as part of the open-source GISAS simulation package Born-Again (Pospelov et al., 2020). All floating-point numbers, internal and external, have double precision. A summary of the algorithm is given in Appendix C.
The code underwent extensive tests for internal consistency and for compatibility with conventional form factor formulae. Checks of BornAgain against the reference code IsGISAXS (Renaud et al., 2009;Lazzari, 2006) have been documented by Pospelov et al. (2020). In the following, we describe form factor consistency checks that have been permanently added to the BornAgain unit tests.

Tests
The internal consistency tests address symmetry, specialization and continuity. Symmetry tests are performed for particle shapes that are invariant under some rotation or reflection R. For a suite of wavevectors q, it is checked that the relative deviation of form factors F(q) and F(Rq) stays below a given bound.
The specialization tests address pairs of figures Å 1 , Å 2 with different topologies that coincide for certain geometry parameters. For instance, if Å 1 is a box with side lengths a, b, c, and Å 2 a truncated cube with side length a and truncation length t, then the choices a = b = c and t = 0 reduce Å 1 and Å 2 to the research papers same cube. For a suite of wavevectors q, it is checked that the relative deviation of form factors F(q, Å 1 ) and F(q, Å 2 ) stays below a given bound.
The continuity tests search for possible discontinuities due to a change in the computational method. They need special instrumentation of the code, activated through a CMake option and a precompiler macro. Under this option, additional variables tell us whether the analytical expression or the series expansion has been used in the latest form factor computation, and, if applicable, at which expansion order the summation was terminated. For a given directionq q, bisection is used to determine wavevectors where one of these variables changes. Then, the form factor F is computed for wavevectors slightly before and slightly after the transition, and it is checked that the relative step in F remains below a given bound.
All these tests are performed for a suite of particle shapes, for different wavevector directionsq q with different degrees of symmetry, for a logarithmically wide range of magnitudes q and for a range of complex phases.

Crossover metaparameters
For large q, the polyhedral form factor is computed from (21) with C = 0. For small q, we use (26) with the expansion (28). Therefore, we need a heuristic metaparameterq q that determines which algorithm to use. For large q, it still can happen that q k is small. Therefore, a second metaparameterq q k is needed to determine whether face form factors are computed from the closed expression (9) or from (16) with the expansion (19). Asq q andq q k are dimensionless, the choice of algorithm is based on qr and q k r, where r is the radius of the circumscribing sphere of figure Å. Under a multitude of tests, we obtained the best results withq q ¼q q k ' 10 À2 .

Accuracy
Currently, the bounds for maximum relative form factor discrepancies are 10 À11 in symmetry tests, 6 Â 10 À12 in specialization tests and 6 Â 10 À9 in continuity tests. Discrepancies reaching the order of magnitude of these bounds are only observed for a few out of hundreds of thousands of test cases. Most often, errors are smaller than 10 À15 , i.e. a small multiple of the machine precision. Some of the larger discrepancies are compiler or processor dependent.
The cases of relatively large discrepancy that we have investigated so far all involve special wavevectors that make the integral (1) more symmetric than the underlying figure Å. Appendix D presents one such case: a pyramid that acquires the inversion symmetry of a bipyramid if q lies in the base plane.
It remains to be seen whether such cases warrant closer attention and improved code. So far, we have not encountered a single q, Å combination where symmetry, specialization or continuity tests revealed numeric errors larger than singleprecision machine error.

APPENDIX A Polygon form factor in the literature
As a complement to Section 2.3, we demonstrate the equivalence of our form factor (9) with the result (15) of Lee & Mittra (1983) and Croset (2017). We start from (9), choose c = 0 and expand the unit in-plane vector b q k q k : Expanding the sinc function, and using (2) and (3) to simplify R j AE E j , we obtain Usingn nE ¼ 0, and hence qE = q k E, Usingn nq k ¼ 0, we obtain (15).

APPENDIX B Polygon area as lowest expansion coefficient
As a complement to Section 2.4, we demonstrate the equivalence of two different expressions for the area of a polygon: the q k ! 0 limit of the generic form factor (9), and the surveyor's formula (10). For any constant directionq q with q k 6 ¼ 0, we have from (16) and (19) Splitting q = q ? + q k , drawing the constant q ? r ? in front of the sum and using P E j ¼ 0 we obtain Inserting the definitions of E j and R j , Multiplying out and shuffling indices ðjÀ1Þ ! j for some terms under the sum gives Usingn nq k ¼ 0, we recover (10).

APPENDIX C Algorithm summary
In this appendix, we summarize the algorithm for the computation of the form factor F(q; Å) as derived in this work. For clarity and brevity, we only consider a polyhedron Å without any inversion symmetry. For polyhedra with inversion symmetry, and for other details omitted here, see the actual implementation in the open-source code BornAgain.
Read the wavevector q, the topology T , the vertex coordinates C and the symmetry flags (Section 3.1, ignored here). Discard faces with zero or negligible area. For each face, merge adjacent vertices with zero or negligible distance. Assert that all remaining faces are planar. Compute the circumscribing radius r of Å (Section 4.3). For each face À k , compute the circumscribing r k .
If qr < 10 À2 (Section 4.3), then compute F according to (26) and (28), with changed summation order so that the outermost summation runs over the expansion power n. Terminate if the relative contribution of two subsequent expansion terms is less than = 2 Â 10 À16 .
Otherwise (qr !q q), compute F according to (21) with C = 0. This is a weighted sum of polygonal face form factors f(q, À k ). They are computed as follows: If q k r k < 10 À2 , then compute f according to (16) and (19) with changed summation order so that the outermost summation runs over the expansion power n. Terminate if the relative contribution of two subsequent expansion terms is less than = 2 Â 10 À16 .

APPENDIX D Additional symmetry at special wavevectors
As a complement to Section 4.4, we give one example of how special wavevectors cause extra symmetries in the integral (1) that lead to cancellation and roundoff errors. Consider a pyramid Å with a base À 0 that has at least a twofold rotation axis, as considered in Section 2.5. Letn n be the normal of À 0 , pointing towards the outside of Å. With (4) and (5) it induces the decomposition r = r k + r ? . Now consider a wavevector q in the plane of À 0 so that q ? = 0. To see how this causes an extra symmetry in the integral (1), consider the bipyramid Å 2 that is the union of Å and the mirror image of Å under reflection about the À 0 plane, Å 2 ¼ fr j r k þ r ? 2 Å _ r k À r ? 2 Åg: Obviously, for the special in-plane q under consideration, Thereby, the form factor integral (1) for F(q, Å) has the extra symmetries of the bipyramid Å 2 . The combination of the S 2 symmetry assumed for À 0 and the mirror symmetry that constitutes the bipyramid implies an inversion symmetry around the center of À 0 . As discussed in Section 3.4 this leads to cancellation of some terms. Unless these terms are removed from the implemented formulae, severe roundoff errors must be expected.