The Mathematics of Smoothed Particle Hydrodynamics (SPH) Consistency

Since its inception Smoothed Particle Hydrodynamics (SPH) has been widely employed as a numerical tool in different areas of science, engineering, and more recently in the animation of fluids for computer graphics applications. Although SPH is still in the process of experiencing continual theoretical and technical developments, the method has been improved over the years to overcome some shortcomings and deficiencies. Its widespread success is due to its simplicity, ease of implementation, and robustness in modeling complex systems. However, despite recent progress in consolidating its theoretical foundations, a long-standing key aspect of SPH is related to the loss of particle consistency, which affects its accuracy and convergence properties. In this paper, an overview of the mathematical aspects of the SPH consistency is presented with a focus on the most recent developments.


INTRODUCTION
The basic form of Smoothed Particle Hydrodynamics (SPH) was introduced more than 40 years ago by Gingold and Monaghan [1] and Lucy [2] for solving the equations of fluid dynamics in the context of astrophysical flows. Since then, SPH has become a popular numerical method due to its countless applications to the study of a great variety of astrophysical processes along with its ever-increasing applications to fluid and solid mechanics, magnetohydrodynamics, cosmology, geophysics, engineering, health sciences, and computer-generated animations, among others.
SPH is a Lagrangian scheme that is based on particle interpolation to compute smooth field variables. Such particles 1 act as control masses and carry all physical properties of the system to be simulated. Because of its Lagrangian nature, SPH has clear advantages over traditional meshdependent Eulerian methods. For example, it does not suffer from mesh distortions that affect the numerical accuracy in simulations of large material deformations. Most importantly, advection is performed exactly and therefore material history information can be tracked free of numerical diffusion. While these advantages cast SPH as a potential method for solving many kinds of problems, the method in its classical flavor suffers from slow numerical convergence due to loss of particle consistency, which has motivated in the last 20 years different modifications and corrections to the original formulation [3][4][5][6][7][8][9][10][11][12][13]. Consistency is a mathematical concept that is related to how well the discrete equations approximate the exact differential equations [14]. In other words, it is a measure of the local truncation error carried by the numerical discretization of the governing equations. In SPH, consistency is always approached under the assumption that the numerically approximated function is sufficiently smooth between particles. This is a quite acceptable assumption because information of the function is not available at scales smaller than the interparticle separation distances. Unlike traditional grid-based methods, the SPH discretization is performed in two separate steps: the kernel approximation and the particle approximation. Therefore, any error analysis will involve a two step process. In particular, the loss of particle consistency in standard SPH arises from the particle discretization procedure and the discrepancy that emerges between the kernel and the particle approximations is referred to as particle inconsistency [6]. In general three different sources of particle inconsistency have been identified, namely: truncation of the kernel support due to the presence of a physical boundary, irregular distributions of particles, and spatially varying smoothing lengths in adaptive calculations [6]. In particular, the problem of particle inconsistency in the presence of twodimensional irregular boundaries was revised more recently by Fourtakas et al. [15], where C 0 consistency is restored by discretizing the wall by means of a set of virtual particles.
Early corrective strategies to restore particle consistency were first advanced by Li and Liu [16] and Liu et al. [17], where the kernel function is modified to ensure that polynomial functions are exactly interpolated up to a given degree. Based on a variational formulation of SPH, Bonet and Lok [18] proposed kernel gradient corrections to allow exact evaluation of the gradient of a linear function, while Liu et al. [19] proposed a reconstruction approach of the kernel function that obeys the consistency conditions of SPH in continuous form (i.e., the conditions that the moments of the kernel must satisfy to guarantee completeness). However, this latter approach was quickly put aside because the reconstruction procedure may result in partially negative, non-symmetric, and nonmonotonically decreasing smoothing functions, compromising the stability properties of SPH. More stable schemes for restoring particle consistency rely on Taylor series expansions of the kernel approximation of a function and its spatial derivatives. For example, if m derivatives are retained in the series expansion, then the resulting kernel and particle approximations will have (m + 1)th-order accuracy, or C m consistency. Based on this observation, Chen et al. [3,4] proposed the Corrective Smoothed Particle Method (CSPM), which is obtained by ignoring all derivatives in the expansion of the function and retaining only first-order terms in the expansions of the derivatives of the function. This method restores C 1 kernel consistency everywhere except at the boundaries, where the kernel consistency drops to C 0 . However, C 1 particle consistency is restored only for uniformly distributed particles, while at the boundaries and for irregular particle distributions the particle approximation is only C 0 -consistent. The CSPM approximations obtained this way have the same form of a Shepard's interpolation for the function and its gradient [20]. Improvements to the CSPM approach were further introduced by Liu and Liu [6]. Their FPM scheme was designed to restore C 1 particle consistency everywhere by retaining first-order terms in the Taylor series expansions for the function and its derivatives and allowing the simultaneous solution of the resulting set of linear equations. An even better approach was proposed by Zhang and Batra [5] (their MSPH scheme) by retaining up to second-order terms in the Taylor series expansions. This method was claimed to restore C 2 particle consistency (i.e., third-order convergence rates) for the function and C 1 particle consistency for the first-order derivatives. However, a drawback of this approach is that when adding higher-order derivatives to the Taylor series expansions, the number of algebraic linear equations that must be solved simultaneously increases. Since the solution involves a matrix inversion, this not only affects the computational efficiency when large matrices are tried, but also compromises the stability of the overall scheme in situations where the matrix becomes ill-conditioned. A FPM-like scheme, which is free of kernel gradients and results in a symmetric corrective matrix was introduced by Huang et al. [10].
Alternative formulations to methods relying on Taylor series expansions were also reported in the literature. For example, Ferrand et al. [21,22] proposed an approach based on a Shepard renormalization of the boundary integrals in the kernel approximation of the spatial derivatives in order to restore C 0 consistency close to physical wall boundaries in weakly compressible SPH calculations, while Macià et al. [23] extended the scheme for incompressible SPH formulations, where a Poisson equation is solved for the pressure field, and performed a consistency analysis. A consistent Shepard interpolation for both the kernel and the kernel gradient was recently introduced by Reinhardt et al. [24] for applications to fluid animation. These authors noted that in the classical Shepard interpolation scheme the correction factors and the fluid quantities are calculated using different kernels, i.e., the uncorrected kernel is used for the correction factors, while the corrected kernel is employed in the calculation of all fluid quantities. This results in an inconsistent scheme since new errors are introduced, which may distort the density field even more than in simulations without a kernel correction. In this formulation C 0 consistency is restored by correcting the volumes of source and neighboring particles simultaneously. A further recent contribution to SPH consistency from the computer graphics and animation community consisted in the formulation of a physically consistent implicit viscosity solver for incompressible and highly-viscous flows [25]. This method produces realistic animations of very complex flows as, for example, the buckling and rope coiling effects observed in nature when a thin stream of honey or caramel flows from a spoon. Such improvements are not only useful for animation purposes, but also for other science and engineering applications. Improved simulations of highly-viscous free-surface flows were recently reported by Kondo et al. [26]. A completely different scheme based on a piecewise high-order Moving-Least-Squares WENO reconstruction and the use of Riemann solvers was presented by Avesani et al. [27]. In particular, this latter approach was found to improve the accuracy and stability of SPH in the presence of sharp discontinuities.
As was first realized by Rasio [28] and then analyzed by Zhu et al. [13], a new source of inconsistency is due to the use of small numbers of neighbors per particle in simulations with compactly supported kernels. A formal derivation of the error carried by the SPH interpolation by Sigalotti et al. [29] has shown that the inconsistency is removed by increasing the number of neighbors. This was thereafter tested in the astrophysical context with reasonably good results [30,31]. In this approach no corrections to the SPH discretization are required and full particle consistency is restored in the asymptotic limit N → ∞, h → 0, and n neigh → ∞ with n neigh /N → 0, where N is the total number of particles, h is the smoothing length, and n neigh is the number of neighbors within the kernel support. Unlike most modern applications of SPH, Zhu et al. [13] demonstrated that C 0 particle consistency can be achieved when working with a large number of neighbors and correspondingly small values of h. This result is consistent with the analysis performed by Read et al. [32], who found that working with a small number of neighbors, as it is common practice in SPH simulations, particle consistency is lost due to irreducible zerothorder errors that will appear even though N → ∞ and h → 0. In spite of these improvements, the issue of convergence in a rigorous mathematical sense and its relation to the conservation properties have not been well understood. However, encouraging preliminary steps in this line were first put forward by Di Lisio et al. [33], Moussa and Vila [34], Quinlan et al. [35], Vignjevic et al. [36], Vaughan et al. [37], Fatehi and Manzari [38], and more recently by Sigalotti et al. [39], who developed an error analysis in n-dimensional space using the Poisson summation formula and derived for the first time the functional dependence of the error bounds on the SPH interpolation parameters, namely N, h, and n neigh . In the following, an overview of the mathematical aspects of the SPH consistency will be given by focusing on the most recent advances. The paper is organized as follows. The standard SPH formulation along with the problem of particle inconsistency are briefly described in Section 2, while Section 3 deals with an overview of the most popular corrective schemes. Section 4 presents and discusses new insights and developments that have been made to improve the consistency of standard SPH. A brief discussion on the strategies that SPH practitioners must follow to achieve consistency in their applications is given in Section 5, while Section 6 contains the conclusions.

CONVENTIONAL SPH FORMULATION
The foundation of SPH is interpolation theory. The starting point of the method lies on the Dirac sifting property where f(x) is some smooth function defined in the domain Ω ⊆ R n and evaluated at the generic point x ∈ Ω, δ(x − x′) is the Dirac delta distribution, and n denotes the spatial dimension (n 1, 2, or 3).

Kernel Approximation
The kernel approximation is obtained by substituting the Dirac delta distribution by an interpolation function such that where is the interpolation or smoothing kernel function, and h is the smoothing length. In general, W is a bell-shaped function defined such that: (1) it must satisfy the normalization condition 2) it must equal the Dirac delta distribution in the limit h → 0, 3) it must be positive definite, 4) it must be an even function such that W( x − x′ , h) W( x′ − x , h), and 5) it must be a monotonically decreasing function. In addition, modern applications of SPH employ kernels of compact support for which W 0 except within a sphere of radius kh from the center x, where k is some integer that depends on the kernel. Hence, W 0 for x − x′ ≥ kh. If f(x′) in the integrand of Eq. 2 is expanded in Taylor series about x, the kernel approximation becomes [39].
where ∇ (l) is the product of the ∇ operator l times with respect to coordinate x, the symbol ::. . .: is used to denote the lth-order inner product, and (x′ −x) l is a tensor of rank l. It is clear from Eq. 4 that 〈f(x)〉 → f(x) only if the family of consistency relations are satisfied for the moments of the kernel for l 1, 2, . . . , where 0 (1) (0, 0, 0) is the null vector and 0 (l) is the null tensor of rank l. C 0 kernel consistency is always guaranteed by virtue of the normalization condition Eq. 3, while C 1 kernel consistency is automatically ensured because for l 1 Eq. 5 is always satisfied due to the symmetry of the kernel. The same is also true for all odd values of l. Therefore, only for l even the above integrals do not necessarily vanish and for l 2 a second-order error follows for the estimate of the function. Therefore, C 2 kernel consistency is not achieved by the kernel approximation unless a higher-order kernel is used for which M 2 0 (2) [40]. Similar expressions to Eqs 2, 5 for the estimate of the gradient and its moments can be derived as and for l 2, 3, . . . , where I is the unit tensor of second rank. Therefore, exact approximation to the gradient requires that the kernel consistency relations given by Eq. 7 are exactly fulfilled.

Particle Approximation
The SPH discretization makes reference to a set of Lagrangian particles which may, in general, be disordered as a consequence of their being convected with the fluid. The discretization is shown schematically in Figure 1, where the domain Ω is divided into N subdomains, Ω k , within each of which lies a Lagrangian particle at point x k chosen such that x k ∈ Ω k . It is common practice to define the boundaries of the Lagrangian subdomains, i.e., zΩ k , such that their masses remain constant. The kernel support is termed Ω s , which may be a circle in two-dimensions or a sphere in threedimensions, and is centered at the observation point x a . Thus, the numerical integration is performed over the intersection set between Ω s and the finite model domain Ω. The elements of the intersection set are the neighbors around particle a (or subdomain Ω a ). Using the mean value theorem so that x b ∈ Ω b for each neighbor b of the observation particle a and denoting by V b the volume of Ω b , Eq. 2 in discretized form becomes where the last term in Eq. 8 is called the summation interpolant. From now on we will simplify the notation by setting W ab W( x a − x b , h) and n a n neigh for the number of neighbors of particle a, which should not be confused with the same letter used for the dimension. Here the subscript a is used to denote a generic observation particle and the subscript b is used to denote its neighbors. Therefore, in simplified form the particle estimate of function f(x) at the observation point x a becomes as it is traditionally written in the SPH literature [41][42][43][44]. Here N a {1, 2, . . . , n a } is the set of all neighbors of particle a and the second summation stresses the fact that subscript b takes only values pertaining to this set. In general, all particles have the same number of neighbors so that the number of elements of N a is always the same (i.e., n n a ). In actual SPH applications it is customary to replace the particle volume V b by the ratio m b /ρ b , where m b and ρ b are the mass and density of particle b, respectively. For quasi-ordered particles following a lowdiscrepancy sequence, Fulk [45] proved that the SPH approximation (9) is consistent with the identity operator, is the exact value of the function evaluated at point x a , while Kf 〈f(x a )〉 and Sf f a are, respectively, the kernel and particle estimates of the function at x a defined by Eqs 2, 9. In the limit of vanishing inter-particle distances Sf − Kf ∞ → 0, whereas Kf − If ∞ → 0 as h → 0. Therefore, it has always been common practice to write Eq. 9 as Following similar steps, the particle approximation of Eq. 6 can be written as A better representation of the gradient, which vanishes exactly when the function is a constant, is given by [46] ∇f a b∈N a Other alternative representations can be found in the literature guided by momentum conservation requirements [41,42].

Particle Inconsistency
The SPH approximation given by Eq. 11 introduces two types of error. The first one comes from the kernel approximation which is ∝ h m , where m is a small integer that depends on the order of completeness that is enforced. This error is independent of the FIGURE 1 | Schematic of the SPH discretization with a compactly supported kernel function.
Frontiers in Applied Mathematics and Statistics | www.frontiersin.org December 2021 | Volume 7 | Article 797455 particle distribution and for most commonly used kernels it contributes with a second-order term (i.e., m 2). The second error originates from the discretization procedure and is a function of the number of neighbors, ϕ ϕ(n), and depends on how they are actually distributed within Ω s . It has long been conjectured that ϕ(n) ∝ log n/n when the neighbors are quasiordered [46], while if they are randomly distributed then ϕ(n) ∝ n −1/2 [13]. The dependence of the discretization error on n can be parameterized as ϕ(n) ∝ n −p , where 1/2 ≤ p ≤ 1. Zhu et al. [13] noted that the reduction of the overall error would demand making h smaller and, at the same time, n larger so that the combined limit N → ∞, h → 0, and n → ∞, with n/N → 0, would be necessary to achieve full consistency. This limit was first noted by Rasio [28] using a simple analysis based on sound waves. From the matching between the smoothing and particle discretization errors, Zhu et al. [13] proposed the following power-law dependences of h and n on the total number of particles (N) where β − 3 + m/p. Thus, for m 2 and a random distribution (p 1/2), β − 7 and n ∝ N 0.57 , while for quasi-ordered distributions (p 1), β − 5 and n ∝ N 0.4 . For an intermediate value of β(− 6), h ∝ N −1/6 and n ∝ N 1/2 . If higher-order kernels (m > 2) are chosen, a stronger scaling of n with N would be required. From a truncation error analysis of the SPH approximation of spatial derivatives in one dimension, Quinlan et al. [35] suggested that the SPH errors will decay as h 2 if h is reduced while maintaining Δ/h constant, where Δ is the inter-particle distance. If, on the other hand, Δ/h is reduced while maintaining h fixed, the error will also decrease, which is equivalent to increase the number of neighbors. Using the Euler equations, Read et al. [32] realized that the loss of particle consistency is due to the presence of h-independent, zeroth-order error terms that are irreducible even when N → ∞, h → 0, and n is maintained constant at the low values typically used in most SPH applications. They suggested as a convergence criterion that n must increase to ensure that ϕ(n) shrinks faster than the smoothing error (∝ h m ). A similar result was reported by Vaughan et al. [37], who found that if C 0 consistency is not restored the error does not converge with h, whereas if C 1 consistency is achieved then the error goes as −h 2 . The claim of particle consistency based solely on the equivalence V b → m b /ρ b is not sufficient because in general the conditions on the kernel moments given by Eqs 3, 5 are not fulfilled in discrete form, i.e., and similarly for the particle form of the integral conditions in Eq. 7 for the kernel gradient. The loss of C 0 particle consistency is a clear manifestation of the violation of the discrete form of the probability function given by Eq. 15 even for regularly ordered particles because of noise errors that scale with n as −n −1 . In practical applications, as the particles are convected with the fluid, node disorder results in amplified noise, leading to loss of C 0 consistency because of more severe deviations of the particle normalization condition from unity. Loss of C 0 consistency also occurs when Ω s is truncated by a physical boundary or when spatially and temporally adaptive kernels are employed [6].

Conservation Principles
The loss of particle consistency in fluid simulations may affect the fundamental principles of physics, such as the principles of conservation of mass, energy, linear momentum, and angular momentum. For most applications, the conservation of mass and linear momentum represents irreducible requirements. Restoring particle consistency is therefore important to guarantee local conservation because if a condition is satisfied locally everywhere within a domain, then it will be also satisfied globally.
Local and global mass conservation can be satisfied exactly by using Eq. 11 with f → ρ to evaluate the density provided that all particles with Lagrangian boundaries preserve their mass. As was demonstrated by Vignjevic et al. [36], if C 0 particle consistency is achieved, then the homogeneity of space is not affected by the process of spatial discretization, which in turn has, as a consequence, the conservation of linear momentum. In other words, under C 0 particle consistency the SPH interpolation becomes independent of rigid-body translations of the coordinates. To put this in mathematical form, let us first approximate the position vector x a using Eq. 11 The transformed coordinates x′ x + Δx will then be estimated according to where . . , n ′ n a ′ }, and n′ is the number of neighbors in the transformed support domain Ω s ′ of radius kh′. Preservation of space homogeneity under uniform translations demands that where we have made m b ′ m b , W ab ′ W ab , and n′ n because under solid-body translations the coordinates of a point are independent of the translation of the coordinate axes. Therefore, Eq. 19 expresses that homogeneity of the discrete space is satisfied by the SPH interpolation only when the condition M 0,a 1, i.e., only when C 0 particle consistency is restored 2 . Vignjevic et al. [36] also demonstrated that angular momentum conservation is exactly preserved only if the isotropy of space is not affected by the SPH interpolation. This is equivalent to demand that the SPH interpolation must be invariant under solid-body rotation of the coordinate axes. In order to simplify the mathematical operations, let us consider only small rotations so that the coordinates are transformed according to where dw is the differential rotation vector. Under solid-body rotation, the coordinates of a point are invariant with respect to a rotation of the coordinate axes and therefore [∇(dw ×x)] a ∇(dw ×x), or using the particle approximation given by Eq. 12 which implies that the first moment of the gradient given by Eq. 7 in discrete form must satisfy the condition in order to preserve the isotropy of the discrete space and conserve angular momentum. In summary, the particle approximation does not affect the conservation of linear and angular momentum provided that C 0 and C 1 consistencies are restored for the function and its gradient. When thermal effects are important, energy conservation is ensured by the first law of thermodynamics from which an equation for either the internal energy or the enthalpy is derived. Rotational invariance using the more accurate representation given by Eq. 13 requires that, in addition to the condition represented by Eq. 22, the discrete zeroth moment of the kernel gradient, M 0 ′ , must vanish.

CORRECTIVE SPH SCHEMES
As was outlined above, SPH particle inconsistency is a manifestation of the discrepancy between the spatially discretized equations and their corresponding kernel approximation in continuous form. Improvements on this drawback is therefore mandatory for SPH to become a robust and reliable simulation tool. Earlier approaches to restore particle consistency were based on kernel reconstruction techniques [19,47,48]. Although these methods can reproduce a polynomial function and its derivatives to any desired order, they can sometimes lead to reconstructed kernels that may be non-positive, not symmetric, and not monotonically decreasing, yielding inaccurate results. In addition, the computational effort for constructing higherorder reproducing kernel functions could be quite intensive.
More robust and stable approaches rely on a Taylor series expansion of the function and its derivatives in the kernel estimate. In what follows the most popular methods based on Taylor series expansions will be briefly reviewed along with some alternative recent schemes.

The Corrective Smoothed Particle Method (CSPM)
The CSPM scheme was originally developed by Chen et al. [3,4] as a corrective scheme for approximating a function and its derivatives. In this formulation, the conventional kernel estimate of a function is expressed in terms of its Taylor series expansion about an arbitrary point of its domain of definition.
Expansion of the function f(x) about the point x a up to secondorder terms gives Multiplying the above expression by the kernel function W a W( x a − x , h) and integrating over the domain volume V, the result is which after solving for f(x a ) and neglecting all terms containing derivatives, the corrective kernel estimate of the function becomes which is identical to the Shepard kernel interpolation [20]. From now on, the angle brackets used to denote the kernel estimate of the function will be dropped for simplicity. The approximations for the derivatives are constructed the same way by multiplying Eq. 23 by ∇W a for the first derivatives and by ∇∇W a for the second derivatives and then integrating both expressions over the whole volume of the domain. Retaining terms up to first-order, the expansion for the first derivatives becomes Similarly, retaining terms up to the second-order in the expansion for the second derivatives, the result is In three-space dimensions, Eqs 26, 27 represent two independent sets of three and six coupled equations, respectively. Due to the symmetry of W, the second term on the right-hand side of Eq. 24 vanishes identically for points that are far from the boundaries of the domain, implying that the estimate given by Eq. 25 is of second-order for unbounded domains. However, for points close to a boundary Eq. 25 is only a first-order approximation because the second term on the right-hand side of Eq. 23 no longer vanishes. The same holds for Frontiers in Applied Mathematics and Statistics | www.frontiersin.org December 2021 | Volume 7 | Article 797455 6 the kernel approximation of the derivatives. Therefore, the whole scheme achieves C 1 kernel consistency away from the boundaries and C 0 kernel consistency for points near or on the boundary. For numerical work, the particle approximation of Eq. 25 is just while the particle approximation of Eq. 26 for the corrective firstorder derivatives can be written in terms of the linear matrix equation where A is a 3 × 3 matrix, F′ is the vector solution whose elements are the three derivative components, and b is a column vector. Using index notation, these quantities are given by where in three dimensions l, m 1, 2, and 3. Since the integral on the left-hand side of Eq. 27 is a tensor of rank four and the other two integrals on the right-hand side are of second rank, the six components of the second-order derivatives can be obtained by solving the linear equation where in index notation for l, k, s, m 1, 2, and 3 in three dimensions. The solution of Eqs 29, 31 uniquely determine the particle estimate of the three components of the first-order derivatives and the six components of the secondorder derivatives at the position of particle a.

The Finite Particle Method (FPM)
The FPM approach was proposed by Liu and Liu [6] as an improved extension of the CSPM method for the kernel and particle approximations of the function and its first derivatives. The main differences with the CSPM scheme is that terms up to second-order are retained in the expansion given by Eq. 24 for the kernel and particle approximations of both the function and its first-order derivatives and that the particle estimates of the function and its derivatives are solved simultaneously. Using Cartesian coordinates the problem reduces to solve a linear system of the form where A is a 4 × 4 matrix, which in three dimensions looks like In discrete form, the solution vector containing the particle estimate of the function and its first derivatives is given by where and x ba x b − x a , y ba y b − y a , and z ba z b − z a are the components of vector x ba x b − x a in a Cartesian coordinate system. It is often claimed that this correction scheme restores C 1 kernel and particle consistency for both the interior and boundary points, implying full second-order accuracy.
Although it improves the C 0 consistency of CSPM for the particle approximation of the function, it can be shown that the method achieves only C 0 consistency in its discrete form given by Eqs 35, 36 for both the function and its first derivatives [29].

The Modified Smoothed Particle Hydrodynamics (MSPH) Scheme
Although this method was reported by Zhang and Batra [5] in 2004 before the FPM paper of Liu and Liu [6] was published, MSPH is just a straightforward extension of FPM by including second-order terms in the expansion given by Eq. 24 for the function and its first and second derivatives. The coupling of the particle estimates of the function and its first and second derivatives results in a system of 10 equations for the 10 unknowns f a , zf a /zx l a , and z 2 f a /zx l a zx k a at point x a , where l, k 1, 2, and 3 in three dimensions. The details of the method are similar to those displayed for the FPM scheme and therefore they will not be repeated here. Although C 2 consistency is achieved by the kernel estimates of the function and its first and second derivatives, the particle approximation is actually C 1consistent for the estimate of the function and C 0 -consistent for the first-order derivatives. Retention of third-order derivatives in the Taylor series expansions will result in a system of 20 linear equations that must be solved simultaneously for the particle approximation.
This method can be straightforwardly generalized to ensure any order of consistency. A generalization to kth-order was reported by Di Blasi et al. [8], where the estimates of the function and all its higher-order derivatives are obtained with kth-order consistency. Thus, increasing the order of the Taylor series expansions will also increase the computational cost because this will require inverting large corrective matrices. In addition, a serious drawback of this type of schemes is that the corrective matrix may be close to singularity especially when the kernel gradients do not exist or when they are either close or equal to zero. In these cases, the matrix will not admit inversion and the method will fail. More recently, a modification to the MSPH scheme was presented by Huang et al. [10] in which case the kernel gradient is not necessary in the whole calculation. This method which was dubbed "a kernel gradient free (KGF) SPH method" presents some advantages over the conventional FPM and MSPH corrective approaches in that kernels that are not differentiable or not sufficiently smooth can all be used in the KGF-SPH formulation. Also, due to the symmetry of the correction matrix, KGF-SPH is insensitive to the particle distribution, free of ill-conditioning of the correction matrix, and performs better than FPM and MSPH at low spatial resolution. A reformulation of the FPM and MSPH schemes was suggested by Sibilla [12], where only the first-, second-, and third-order derivatives are retained as unknowns. This scheme enforces C 1 particle consistency for the first-order derivatives and C 0 consistency for the second-order derivatives, both inside the domain and close to the domain boundaries. Like the KGF-SPH scheme, it is not affected by irregular particle distributions.

A Boundary Integral SPH Formulation
A different strategy to restore consistency near wall boundaries was first introduced by Ferrand et al. [21,22] in the context of a weakly compressible SPH formulation by just considering the contribution of the boundary integral terms in standard SPH, while Macià et al. [23] extended the scheme to incompressible SPH formulations. Apart from the treatment of the Poisson equation for the pressure, both schemes use the same principle. Therefore, for simplicity the method is described here in terms of the onedimensional analysis reported by Macià et al. [23]. The generalization to multi-dimensions is straightforward. The method starts by considering a symmetric kernel function such that which is the one-dimensional form of Eq. 3. The kernel is assumed to have compact support and the kernel estimate of the function within the interval (A, B) ⊂ R is given by the Shepard interpolation where f(x) is a scalar function on R, W a W(|x a − x|, h), and the angle brackets are re-used to distinguish the kernel estimate of the function from its exact value evaluated at x x a ∈ (A, B). The kernel approximation defined via Eq. 38 has C 1 consistency for all interior particles in the interval (A, B) and C 0 consistency for particles at distances ≤ kh from a domain boundary. For any x a ∈ (A, B) far from a boundary Γ(x a ) 1 in Eq. 38 and the kernel estimate becomes By expanding f(x) in Eq. 39 in Taylor series about x x a and subtracting the result from the exact value of the function evaluated at x x a , it is easy to show that 〈f(x a )〉 − f(x a ) ∝ O(h 2 ), while the same operation on Eq. 38 produces the result because the first moment of the kernel no longer vanishes. Replacing f(x) by df(x)/dx in Eq. 38 produces the following expression for the kernel estimate of the first-order derivative Subtracting this from the exact derivative evaluated at x x a , it can be demonstrated that the error carried by Eq. 41 far from a boundary is O(h 2 ), while in the proximity of a boundary it is O(h). In a completely analogous form, the expression for the kernel estimate of the second-order derivative reproduces the exact derivative to second-order accuracy for unbounded domains and for points far from a boundary, while in the proximity of a domain boundary the approximation is only first-order accurate. In particular, when a gradient is calculated near a domain boundary, C 0 particle consistency is restored for the estimate of the gradient if the boundary terms are included. The same is true for the particle estimate of the second-order derivatives, which also restores C 0 consistency close to a domain boundary only when boundary terms along with the normalization factor Γ(x a ) are included. Mathematically, this can be seen by integrating by parts the kernel estimates given by Eqs 41,42. After straighforward integration by parts Eq. 41 becomes Similarly, after integration by parts of Eq. 42, the kernel estimate of the second derivative becomes where one-sided finite differences have been used to evaluate the first-order derivatives of the function at both extremes of the interval (A, B). Macià et al. [23] suggested the following particle discretization of Eqs 43, 44 and respectively, where the notation f 1 f(x 1 ) f(A) and f N f(x N ) f(B) has been used. Also, note that the index b employed for the neighbors of particle a only runs for interior particles and excludes the lower and upper ends of the interval (A, B), i.e., b ≠ 1 and b ≠ N.

An Integral Approach to First Derivatives
An interesting formalism based on a matrix approach to the SPH Euler equations that improves the interpolation of physical magnitudes was introduced by García-Senz et al. [49]. In this approach first derivatives are approximated from an integral expression that leads to a tensor, instead of the usual vector, estimation of the gradients. The scheme reduces to the standard formulation in the continuum limit and results into a less noisy estimate of the gradients. Assuming that f(x): R 3 → R is a smooth and differentiable function, the method starts with the integral expression Expanding the term f(x′) − f(x) in Taylor series about x′ x and retaining terms up to first-order, the above integral becomes Using Cartesian coordinates and writing where x 1 x, x 2 y, and x 3 z, Eq. 48 can be written as a linear matrix equation of the form that must be solved simultaneously for the three components of the gradient, where for i, j 1, 2, and 3, and for k 1, 2, and 3. For spherically symmetric kernels, A 11 A 22 A 33 and all off-diagonal terms vanish, i.e., A ij A ji 0 for i ≠ j. When the integrals in Eq. 51 are converted into finite sums, most of the symmetry properties of matrix A are lost because the diagonal elements may have different values and the off-diagonal elements may not vanish. However, the formulation of the gradient according to Eq. 50 preserves the property that the derivative of a linear function is always exact by construction. This can be easily seen in one-space dimension, where Eq. 50 in particle form reduces to where Although this method leads to improved results compared to standard SPH, it has only been tested for two-dimensional problems and the authors have commented on some difficulties that will probably come up in three-dimensional simulations. These difficulties include: 1) a much higher level of numerical noise as a result of an increased sensitivity to particle disorder compared to conventional SPH and 2) an inherent difficulty in handling sharp boundaries which is typical of matrix methods.

RECENT INSIGHTS AND DEVELOPMENTS
Other alternative attempts have been recently proposed to improve the convergence properties of standard SPH. Advances in this direction were presented by Adami et al. [9], who proposed a transport-velocity formulation for solution of the fluid-dynamics equations that combines the Frontiers in Applied Mathematics and Statistics | www.frontiersin.org December 2021 | Volume 7 | Article 797455 homogenization of the particle configuration by a background pressure and reduces artificial numerical dissipation. A wide number of numerical tests have demonstrated the accuracy and stability of the scheme. In a further work, Litvinov et al. [11] have outlined the importance of partition of unity as a condition under which the conventional SPH approximation can achieve consistency and convergence. In particular, they argued that the inconsistency of standard SPH approximations for a random particle distribution may be due to the SPH definition of a particle volume, V a m a /ρ a , which does not exactly partition the total volume, as is indeed required by the Quasi Montecarlo method whose error goes as (log N) n /N under randomly perturbed particles, where now n 1, 2, or 3 denotes the dimension. Exact partition of unity requires that that is, the sum of all particle volumes must equal the total volume of the system. In general, this condition is not exactly fulfilled by standard SPH because of either gaps between or overlaps of the particle volumes. Litvinov et al. [11] found that standard SPH will achieve better error properties and improve on the partition of unity by relaxing a particle distribution under a constant pressure field and invariant particle volume. The use of the principle of least action to derive the equations of motion with a regularized density field in the context of measures has been illustrated by Evers et al. [50]. Depending on the order in which regularization and the principle of least action are applied two different equations can be derived, whose discrete counterparts coincide with the traditional scheme used in SPH. Measure theory provides a framework to study convergence as the number of particles goes to infinity. In particular, both the particle approximation and the limiting continuum setting can be formulated in terms of measures. For example, the Wasserstein distance on the space of probability measures provides a natural tool to characterize convergence. Convergence of measure-valued solutions was first studied by Di Lisio et al. [33], who proved the convergence of the SPH method with the use of measures in combination with the Wasserstein distance. The measure-valued formulation of Evers et al. [50] not only generalizes Di Lisio et al. [33] work in that the SPH-particle approach is just a special case because it can be applied to a class of approximating measures that is much broader than just sums of Dirac delta distributions. An important concluding remark from this study is that a favorable approach is to have h depending on the total number of particles N in such a way that h − O(N −1/n ) (where n denotes dimension) as N → ∞, and hence the joint limit N → ∞ and h → 0 must hold simultaneously. This entails a close relationship with the fact that the smoothing kernel asymptotically approaches a Dirac delta distribution when h → 0 and this can be achieved in the continuum limit when N → ∞ simultaneously.

Consistency and Self-Consistency Conditions
The particle density estimate as is frequently used in SPH can be derived from Eq. 11 by replacing f by ρ such that where in the original SPH formulation the sum is taken over all particles of the system. However, in modern applications of SPH a smoothing function of compact support is employed to limit the volume to a locally finite extent in order to improve the computational efficiency. This means that the sum will then operate only on a finite number (n) of particles around particle a. As was pointed out by Zhu et al. [13], consistency of this step requires that n → ∞ in order to approach the continuous limit. Therefore, there is a lack of self-consistency when working with finite values of N and n, where the condition is not satisfied with the density estimate given by Eq. 56 when the sum is made to operate on a finite number of neighbors. In order words, the violation of Eq. 57 implies that partition of unity is not exactly fulfilled by standard SPH as was claimed by Litvinov et al. [11]. Therefore, the joint limit N → ∞, h → 0, and n → ∞, with n/N → 0, is a necessary requirement to achieve consistency and self-consistency [13,28]. This complies with the suggestion advanced by Evers et al. [50] that h → 0 as N → ∞ if the sum interpolant is taken over all particles of the system, as in Eq. 56. A previous error analysis of the SPH representation of the density and momentum equations have shown that for finite values of n, the errors carried by the SPH discretization will not decay even when N → ∞ and h → 0 because of an irreducible zeroth-order error that is independent of h and that will depend only on n [32]. Based on the above discussions, restoring particle consistency in SPH applications is impeded by the loss of selfconsistency when using smoothing functions of compact support with a small number of neighbors. As was already outlined in Section 2.3, Zhu et al. [13] derived power-law dependences of n and h on N complying with the joint limit N → ∞, h → 0, and n → ∞ by matching the SPH error from the kernel approximation (smoothing procedure), ∝ h m with m 2 in most cases, with that resulting from the particle approximation, ϕ(n) −n −p , where 1/2 ≤ p ≤ 1. By balancing both errors, they suggested to use the scaling relations given by Eq. 14 as approximate requirements for SPH convergence. However, increasing the number of neighbors with resolution will demand using a particular class of kernel functions. It is wellknown that conventional kernels, like the widely used cubic B-spline kernel [51], suffer from a pairing instability when supporting large numbers of neighbors. The instability manifests by particles coming into so close pairs that resolution is halved and, as a consequence, they become less sensitive to small perturbations within the volume of the kernel support [52][53][54]. One route to overcome this difficulty has been to use the Wendland functions [52,55]. Wendland functions are a family of spherically symmetric interpolation kernels that are smooth everywhere and have positive Fourier transforms under compact support, thereby allowing the use of large number of neighbors without suffering from particle clumping. Also, the Frontiers in Applied Mathematics and Statistics | www.frontiersin.org December 2021 | Volume 7 | Article 797455 exact particle distribution depends on the actual flow dynamics and on the kernel function that is employed to perform particle interpolation. For example, in contrast to most conventional kernels, Wendland functions do not allow particle motions on a sub-resolution scale (i.e., within the kernel support) and therefore they maintain quasi-ordered particle distributions even in highly dynamical test problems [56]. A Wendland function that has been adopted in many recent tests and simulations [13,30] is the Wendland C 4 kernel function given by for q ≤ 1 and 0 otherwise, where q x − x′ /h. The magnitude of self-inconsistency due to violation of the partition of unity was further examined by Zhu et al. [13] and later on by Gabbasov et al. [30] in three-dimensional simulations of the protostellar collapse, where the full set of gravitohydrodynamics equations were solved. In particular, Zhu et al. [13] found through a simple test that for a random particle distribution the density distribution approaches a Dirac-δ distribution as n is increased, leading to approximate volume partitioning. Also, the standard deviation measured in the density estimate as a function of varying n showed the expected inferred trend of σ(ρ) ∝ n −1/2 . When the particles were relaxed into a quasi-ordered distribution, the density distribution as a function of varying n was much narrower at all resolutions compared to the random case and σ(ρ) ∝ n −1 with increasing n as expected. These results were further confirmed through full hydrodynamical simulations by Gabbasov et al. [30], where particle consistencies close to C 1 were confirmed for values of N above 9 million particles and n higher than 24 thousand neighbors.

Consistent Shepard Interpolation
Since the introduction of SPH to the fascinating field of computer graphics and animation in 1996 by Desbrun and Gascuel [57], the method has been extensively revised and improved to produce realistic animations of deformable bodies and complex flow scenarios. The ever-increasing demand for good visual quality and realism for computer-generated animations has imposed new challenges to many SPH practitioners. Following this line of research, Reinhardt et al. [24] have recently implemented a generalization of the classical Shepard interpolation to achieve constant completeness (or true C 0 consistency) for both interior particles and particles near a rigid boundary. The essence of the new approach consisted in realizing that the classical particle Shepard correction of the smoothing kernel defined as introduces new errors because the densities of neighboring particles entering in the normalization factor c sh a are calculated with the uncorrected kernel, while any fluid attribute, saỹ is calculated using the corrected kernel. This approach is inconsistent and when applied to solve the hydrodynamics equations, the density calculation is inaccurate and produces erroneous pressure forces that may induce instabilities in the course of the simulation. In order to restore consistency with the Shepard interpolation, Reinhardt et al. [24] have proposed to use the corrected kernel,W ab , to compute the correction factors, c a , in Eq. 59. To do so, it is necessary to correct the densities of all neighboring particles, which is equivalent to correct their volumes according to the prescription V b m b /(c b ρ b ) such that the correction factors become which must be solved implicitly. They found that use of the power method provides an unconditionally stable and fast converging algorithm for the calculation of the correction factors. To this end, they proposed working with the equivalent form of Eq. 61 where a ab m b W ab /ρ b are the elements of matrix A (a ab ) (a,b)∈N×N , which satisfy the condition 0 ≤ a ab ≤ 1. Solving Eq. 62 is equivalent to solve the linear system A ·c c, which implies finding an eigenvectorc to the eigenvalue λ 1 of matrix A, where the elements of vectorc are the correction factorsc a for a 1, 2, . . . , N. Because of the condition 0 ≤ a ab ≤ 1, matrix A is a column stochastic matrix and so the eigenvalue λ 1 exists, while for all other eigenvalues |λ| < 1. This allows the use of the power method to find the solution vectorc through the recurrence relatioñ For details of the mathematical derivation and convergence proof of the method, the interested reader is referred to Reinhardt et al. [24] paper. Once the correction factors are determined, the corrected kernel follows asW which can be used to perform the SPH simulation step. In order to restore full C 0 consistency, the kernel gradients must be corrected as well. From Eq. 64 it follows that where from Eq. 12, the gradient of the correction factors has the form ∇c a b∈N a Although the use of Eqs 65, 66 increases the accuracy, they do not ensure full C 0 consistency because in general the discrete form of the zeroth moment of the kernel gradient, M 0,a ′ ≠ 0 (1) , implying that the gradient of constant functions is not exactly reproduced.
In fact, reproduction of the consistency condition M 0,a ′ 0 (1) would require that ∇c a , as given by Eq. 66, vanishes identically. However, in general this is not true. To overcome this trouble, Reinhardt et al. [24] have proposed to adjust the gradient of a fluid quantity A a according to where, for instance, the form is suggested when animating the fluid with a weakly compressible SPH formulation. For particles near a solid boundary, C 0 consistency is restored by adjusting the density according to where N fluid a is the set of neighbors of particle a pertaining to the fluid domain, N boundary a is the set of boundary neighbors of particle a, ρ 0 is a reference density, and Ψ b k (ρ 0 ) is a pseudomass given by the expression A consistent approach for particles near the boundary is maintained by including the boundary particles into the calculation of the correction factors, which is accomplished by replacing the uncorrected kernels in Eqs 69, 70 by the corrected ones such that the density of a boundary particle, namely ρ ba , is given by and With the use of the above prescriptions smooth density distributions can be obtained at the rigid boundaries of a body with only very small deviations from the reference density ρ 0 .

The Poisson Summation
The Poisson summation formula was first used by Monaghan [41] as a tool to estimate the errors carried by the SPH summation interpolant for a linear function, f(x) α + βx, in one-space dimension using equidistant particles and a Gaussian kernel However, no clear conclusion about the SPH particle consistency resulted from this analysis. The difficulty was problably due to the fact that it was not until 2016 that it was discovered that the Poisson summation formula is not unique, but belongs to a wider class allowing the construction of crystalline measures with locally finite support and spectrum [58]. This brought out the Poisson summation formula as a powerful tool for the error analysis of particle methods involving the evaluation of quadratures. As was recently demonstrated by Sigalotti et al. [39], the use of this formula will enable the simultaneous treatment of both the kernel and the particle approximation errors for arbitrary particle distributions. In particular, the method allowed to derive for the first time the explicit dependence of the error bounds on the SPH interpolation parameters, namely N, h, and n, from which the consistency scaling relations introduced by Zhu et al. [13] arise naturally. In Sigalotti et al. [39] the analysis is first performed in one dimension and then its generalization to n dimensions is fully developed. Here we shall describe the salient features of the analysis and its results in three-space dimensions, which is the case of interest in most real life applications. As a starting point, let Λ ⊂ R 3 be a crystalline lattice and Φ(x): R 3 → R a smooth function of locally finite support Γ belonging to space D(R 3 ) of all smooth functions with compact support in R 3 . In passing, we note that every function of D(R 3 ) belongs to the Schwartz space, S(R 3 ), of all infinitely continuous functions on R 3 with fast decay at infinity along with all derivatives. The distributional Fourier transform of Φ, namelyΦ, in the dual lattice Λ p is given by the Poisson formula [58]. b1,b2,b3∈Λ where the triplet of integers (b 1 , b 2 , b 3 ), with b i ∈ Z(i 1, 2, 3), denotes the projections of the lattice node (or particle) labels b ∈ Z on the axes of a 3-dimensional Cartesian coordinate system and j (j 1 , j 2 , j 3 ) belongs to the dual of the space labels b. If we set the summation on the left-hand side of Eq. 74 equal to the righthand side of Eq. 8 for the particle estimate of a funtion f(x): R 3 → R at particle position x a , the Poisson summation formula becomes where J x b is the Jacobian of the transformation whose determinant is given by where V 4πk 3 h 3 /3 is the spherical volume of the kernel support and kh is its radius. Note that the first equality in Eq. 76 only holds in the SPH sense. Since it is common practice in SPH to assume the mass of particles to be constant, the particle volumes will change according to the value of the local density. If the local density increases, the particle volume shrinks. If, on the other However, solving for the ratio m a /ρ a in Eq. 80 yields the expression m a ρ a 4π 3 where the volume of the kernel support is V 4πk 3 h 3 /3. Hence, Eq. 83 becomes If the ratio m a /ρ a is associated with the volume of particle a, V a , it follows from Eq. 84 that the sum over all particles lying within the kernel support of particle a is given by where it has been assumed for simplicity that all particles have the same mass. Equation 85 implies that for sufficiently small (but finite) sizes of the kernel support, and independently of whether h is fixed or variable, partition of unity can be achieved only in an approximate sense. Only in the continuous limit when h → 0, exact partition of unity can be fulfilled provided that the joint limit n → ∞, m → 0, and N → ∞ is satisfied. Thus, Eq. 85 differs from Eq. 55 in that for compact kernel supports exact partition of unity demands that the volumes of all neighboring particles exactly fit the volume of the kernel support. If this requirement is fulfilled, the conditions M 0,a 1 and M 0,a ′ 0 will be exactly satisfied. For example, Zhu et al. [13]  numerically the deviation from an exact partition of unity by evaluating the standard deviation measured in the distributions of M 0,a as a function of n for a compactly supported Wendland kernel and a low-discrepancy set of particles, finding that σ(M 0,a ) − n −1 as is indeed predicted by the error bound given by Eq. 81. Therefore, the inconsistency in the volume estimate declines as the number of neighbors increases because in this limit h must also decrease. As this limit is observed, it was demonstrated numerically that the particle approximations of the derivatives converge essentially at the same rate as the particle estimate of the function, regardless of the particle distribution [29].

SOME STRATEGIES FOR SPH PRACTITIONERS
The issue of consistency restoring in SPH has advanced to the point that the SPH practitioners can choose from a number of possible strategies to ensure convergence and consistency of their model simulations. When dealing with incompressible fluid flows, C 0 particle consistency may be enough to produce acceptably good results at a relatively low computational cost. Such schemes include the FPM approach of Liu and Liu [6] and the boundary integral SPH formulation introduced by Ferrand et al. [21] and Macià et al. [23]. However, in the light of the most recent developments, the consistent Shepard interpolation formulation advanced by Reinhardt et al. [24] seems to be a better choice for restoring C 0 particle consistency and improving the volume preservation of the fluid. For example, the consistent Shepard interpolation scheme can be combined with the implicit viscosity solver introduced by Weiler et al. [25] to model complex physical effects in highly viscous fluids, such as rope coiling and melting. Very recent progress in improving the physical consistency of highly viscous, free surface flow calculations has been reported by Kondo et al. [26]. Unlike incompressible fluids, with the exception perhaps of turbulent flows, compressible fluid flows are likely to develop small-scale structures. Such flows are commonplace in many astrophysical phenomena, in the transport of gas through pipelines, in the motion of airplanes and similar highvelocity vehicles, in shock wave formation, and in commercial and industrial applications such as abrasive blasting and supersonic wind tunnels, among many others. For example, formation processes in astrophysics involve flows over many orders of magnitude increase in density and temperature by simultaneously spanning a wide range of spatial and time scales. As a good example, the star formation process involves a continuous gas flow from galactic scales (measured in kiloparsecs (kpc), where 1 kpc 3.086 × 10 16 km) down to stellar scales (measured in astronomical units (AU), where 1 AU 1.496 × 10 8 km). In this case and in cases where compressions, expansions, and shock waves form in small-scale systems, it is important to solve the smallest features that compose the structure of the flow. Resolution requirements become an important factor in this kind of simulations. However, for such flows particle consistency demands using a large number of neighbors and correspondingly small values of the smoothing length according to the power-law scalings advanced by Zhu et al. [13], which define the number of neighbors and the size of the kernel support in terms of the total number of particles filling the computational domain. As the number of neighbors is increased, the minimum resolvable mass is reduced which is a key aspect for solving small local structures in compressible flows at all scales. In these simulations consistency and volume partitioning is improved as the number of neighbors per particle is increased and the smoothing length is reduced with increasing total number of particles. This has, as a consequence, increased accuracy and a reduction of the minimum resolvable mass, leading in most cases to better than C 0 completeness [30].

CONCLUSION
Smoothed Particle Hydrodynamics (SPH) is a Lagrangian particlebased method that has emerged in recent years as a popular numerical technique for the simulation of a large variety of problems in computational fluid and solid mechanics and related areas. On the surface SPH looks like deceptively simple and robust, but when we delve into its convergence properties, in spite of its widespread applications and recent progress in consolidating its theoretical foundations, the method still has unknown properties that need to be investigated. Undoubtedly a long-standing and non-trivial problem that has concerned many practitioners is the lack of mathematical consistency that the SPH approximation typically experiences. The lack of consistency is a key feature of SPH because it not only affects its accuracy and convergence properties, but also the conservation of fundamental principles of physics, such as the principles of conservation of linear momentum, angular momentum, and energy. In the light of these shortcomings and deficiencies, the issue of SPH consistency has become a very hot and important topic of research.
Several corrective methods and strategies have been proposed over the years to restore SPH consistency and misunderstandings regarding its concept have been clarified. In this overview, we have provided a comprehensive survey of the most successful corrective methods and strategies that have been recently developed to address the problem of SPH consistency. In particular, methods based on Taylor series expansions of the kernel approximations of a function and its derivatives can be formulated to achieve any degree of consistency. However, the improved accuracy of these corrective schemes comes at the price of involving the inversion of large matrices that apart from implying an increased computational cost for time-evolving simulations, may also involve a loss of numerical stability due to matrix illconditioning for some specific problems. Other methods consisted in generalizing the kernel approximation of a function and its derivatives by including the boundary integrals and in specialized representations of the kernel and particle estimates of the derivatives. It was not until recently that the Poisson summation formula was employed to study the convergence properties of SPH, allowing for the derivation of the functional dependence of the error bounds carried by the particle interpolation formula on the SPH parameters, namely the total number of particles, N, the smoothing length, h, and the total number of neighbors within the kernel support, n. This dependence indicates that the error of the particle discretization decays as n is increased in compliance with previous findings that have pointed to the joint limit N → ∞, h → 0, and n → ∞ as the only possible way to guarantee full SPH consistency. Whereas this route seems quite promising, it is not free of difficulties. For example, regardless of how large N and n may be and of how small the size of h may be consistency will be restored only in an approximate sense. Although, the use of large values of n along with small sizes of the kernel support may result in drastic improvements of the SPH accuracy and convergence, they may incur high computational costs which, on the other hand, can discourage many practitioners. However, at the present state, the SPH practitioner can choose from a number of possible strategies to mitigate the lack of consistency in their model applications. Unquestionably the improvement of SPH certainly rests on a more accurate volume estimate, or in other words, on the elimination of the zeroth-order error so that the resulting scheme will be confidently second-order accurate because of the symmetry of the smoothing function.

AUTHOR CONTRIBUTIONS
LS wrote the paper, while JK and MG. worked on improving the manuscript and adding some sections. All authors agree to be accountable for the content of the present work.

ACKNOWLEDGMENTS
One of us (LS) acknowledges support from the Department of Basic Sciences of the Autonomous Metropolitan University (UAM-A) of Mexico. We are grateful to the reviewers who have raised a number of suggestions and comments that have greatly improved the style and content of the manuscript. In particular, we thank one of the reviewers for having directed our attention to the newest research on consistency in the field of computer graphics and animation.