Smoothed Particle Hydrodynamics (SPH) is a spatial discretisation method to solve partial differential equations. It is a mesh-free, Lagrangian method in which the system is represented by a finite set of particles. SPH was originally developed to solve astrophysical problems   , but many advantages of the method compared to conventional methods for specific types of problems made it attract attention in other areas, such as fluid and solid mechanics. For instance, with SPH, free surfaces are very easy to deal with  . Also, SPH is quite straightforwardly applicable to multiphase flows    , flows with fluid- structure interaction   , and so on.
As opposed to the astrophysical problems that SPH was initially applied to, most fluid and solid mechanics problems have solid/physical domain boundaries. As a consequence, the support domain of particles close to these boundaries overlaps with the empty area behind it, which has a major negative influence on the accuracy of the SPH approximations. Consistency is also lost if particles are not uniformly distributed, which is the case in most simulations. Restoring the consistency for functions and first-derivative approximations is fairly straight- forward and can be achieved by using correction terms as in, e.g.  . An extensive comparison between correction methods can be found in  . Correcting second-derivative approximations, however, is more complicated.
Second derivatives (or the Laplacian) appear in viscous terms, conduction equations and in the pressure Poisson equation that is used in incompressible SPH. Flebbe et al.  straightforwardly computed the second-derivative terms by first computing the gradient of the unknown variable and then computing the divergence of the result. This sometimes led to non-physical oscillations in the solution  . Originally, approximations for second derivatives included second derivatives of the kernel function, which are sensitive to particle disorder   . Therefore, several ways to restore the consistency have been proposed in the literature. Many researchers have suggested to use expressions based on first derivatives of the kernel function rather than the second derivatives, because they are less sensitive to the particle distribution   . However, these expressions do not solve the consistency problem near boundaries. Some re- searchers have therefore proposed to include boundary terms in the approximation   . Also, there is quite some literature available on Reproducing Kernel Particle Methods (RKPM) and similar approaches    , which are designed to give approximations up to a desired order of consistency. They can lead, however, to partially negative, non-symmetric and non-monotonically decreasing smoothing functions and are therefore not preferred for hydro-dy- namic simulations  . A different approach was followed by Chen and Beraun  , who present a method in which an estimate for the leading error term is subtracted from the approximation. Other researchers use methods in which a function and all of its desired derivatives are solved simultaneously    . These methods have higher accuracy, but are also computationally expensive.
In this work we only consider methods that use conventional smoothing functions. This excludes the previously mentioned RKPM. We describe the original SPH method to approximate second derivatives, as well as the methods by Chen and Beraun  and Zhang and Batra  . We also propose an improvement to the estimate by Chen and Beraun. Our method is computationally only slightly more expensive than the method by Chen and Beraun, but its accuracy is similar to that of the method by Zhang and Batra. We compare numerical results in one and two dimensions for both uniformly and non-uniformly distributed particles, but are especially interested in the latter, because it most accurately resembles particle distributions from actual simulations. Initial results have been presented in   .
2. Smoothed Particle Hydrodynamics
There are two approaches to derive the SPH equations for fluid flows in the literature. The first one is primarily concerned with choosing a density estimate. Substituting this estimate into the Lagrangian then leads to a system of equations conserving both linear and angular momentum. This approach is explained in detail by Price  . The second approach considers SPH to be a spatial discretisation scheme which can be applied to any set of equations, e.g. the Navier-Stokes equations. In this paper we adopt the latter view on SPH. A short derivation of the SPH equations will follow. For more details we refer to, e.g.  .
2.1. The Kernel Approximation
Consider a function , defined on a domain . Then the value of at an arbitrary point can be written as:
Here, is the Dirac delta distribution:
In SPH, is replaced by a continuous function , with a smoothing parameter. is called the smoothing or kernel function and should converge to the Dirac delta distribution as goes to zero. Preferably it is also radially symmetric, has compact support and satisfies the unity condition for all :
Not satisfying the unity condition―to be more precise, its discrete equiva- lent―is the main reason for inaccurate approximations. This will be explained in more detail later. Replacing the Dirac delta function by leads to the following approximation:
which is called the kernel approximation of . Note that it is a weighted average over the values of at all points , including itself.
2.2. The Particle Approximation
To get a numerically useful expression, a quadrature rule is applied to (4). This results in a partitioning of the (computational) domain by a finite number of so-called particles. The summation for a particular particle is limited to those particles that are within the support domain of the particle. This is illustrated in Figure 1. Here, is the support radius. Defining as the set of particles within the support domain of particle , (4) is approximated by:
where is the volume of particle , denotes and . The approximation in (5) is called the particle approxima-
Figure 1. Support domain of a particle illustrated in two dimensions.
tion of . When applied to a physical problem, the volume of a particle is usually replaced by the ratio of its mass and density, , but since this is not necessary for our derivations we will stick to the notation in (5).
2.3. Derivatives in SPH
Approximations for derivatives are found by replacing by its derivatives in (4). For instance, substituting leads to:
Integrating by parts gives:
where is a function depending on and and denotes the boundary of . The first term on the right-hand side of (7) is a boundary term that is usually neglected, since it is zero for particles sufficiently far away from the boundaries. Hence, only the second term on the right-hand side remains and, after discretising, we find:
Similarly, we can derive an approximation for the Laplacian by substituting for in (4). Integrating by parts twice, neglecting boundary terms and discretising leads to:
Note that in both the approximation for the gradient and the Laplacian the derivative operator has switched from to . Since is known, these expressions can be computed.
Substituting a Taylor series expansion for around in (9) reveals that the leading error term of this estimate can be subtracted to find the more accurate estimate:
Here, the subscript “SPH” was added to indicate that this is the standard SPH equation used in Section 5. In practice, (10) is rarely used. The reason is that for most kernels the second derivatives have very steep gradients, making the approximation very sensitive to irregularities in the particle distribution   . An alternative expression is given by Brookshaw  :
where . This approximation uses first derivatives of the kernel function instead of second derivatives, which makes it less sensitive to changes in the particle distribution  . However, it has accuracy issues that are explained in the next section.
2.4. Error Analysis
The simple expressions in (10) and (11) have a downside: their accuracy decreases for particles close to boundaries or when the particle distribution is irregular. This can be shown by substituting a Taylor series expansion for around . For (10), for instance, this gives:
where denotes the Hessian matrix of . The term that we are approxi- mating, the Laplacian, is contained in the second term on the right-hand side of (12).
The first term on the right-hand side of (12) is an -error term that vanishes in the ideal conditions of a uniform particle distribution and no boundaries, but it does not if a particle is close to a boundary or is surrounded by irregularly located particles. A similar -error term can be found for (11).
Furthermore, it is impossible to distill the exact Laplacian from the second term, because the separate second-derivative terms have different coefficients. This leads to an additional -error, which also holds for (11). In the next section we discuss several methods that were designed to give more accurate approximations.
3. Second Derivative Approximations with Higher Accuracy
This chapter describes two methods that approximate the Laplacian with higher accuracy than the conventional estimate in (11). We discuss the Corrective Smoothed Particle Method (CSPM) and the Modified Smoothed Particle Hydrodynamics (MSPH) method.
3.1. Corrective Smoothed Particle Method
This method, described in various papers by Chen and his colleagues      starts from the Taylor series expansion of around :
In all derivations to come, we assume is a two-dimensional domain and and refer to the two independent spatial directions. We stress, however, that it is possible to extend everything in Sections 3 and 4 to three dimensions. To compute second derivatives, the first term on the right-hand side of (13) is first subtracted from both sides of the equation. The result is multiplied by the vector , with defined as:
Multiplying (13) by and summing over all the particles leads to:
In (15) and in the rest of this paper, the derivatives in are taken with respect to . The first term on the right-hand side of (15) is a -term. We wish to subtract it from both sides of the equation, but since is unknown we are forced to use an approximation for that as well. A derivation similar to the one that led to (15), but with instead of , gives:
where is the normalisation matrix defined by:
Multiplying (16) by the inverse of leads to a first-order accurate approximation for the gradient:
The gradient approximation in (19) is substituted for in (15). Subtracting the first term on the right-hand side gives:
It can be shown that the first term on the right-hand side of (20) can be written as , where is the normalisation matrix defined by:
with the following vector:
Multiplying the left-hand side of (20) by the inverse of gives the final CSPM approximation for second derivatives:
By using an approximation for in (15), an error was introduced, as indicated by the approximation symbol in (20). As a result, normalising the first term on the right-hand side of (20) does not lead to the desired first-order accurate approximation. Instead, the first-order error of the gradient estimate makes that the approximation in (23) is zeroth-order accurate, i.e.:
In Section 4 we show how this approximation can be improved to be first- order accurate.
3.2. Modified Smoothed Particle Hydrodynamics
Zhang and Batra   describe a different approach to compute second derivatives. In their method, a vector is computed that consists of an approximation of itself, all its first-order derivatives and all its second-order derivatives:
This is achieved by multiplying (13) by , and , leading to six equations for six unknowns (in two dimensions) for each particle. Because all unknowns are put in one single vector, this method requires the inversion of - matrices. This is computationally expensive, but it leads to more accurate results than with the method described in Section 3.1. In fact, if we isolate from , we find:
It is possible to achieve similar accuracy without using these larger matrices. For that we need only one adjustment to CSPM. This is described in Section 4.
4. Improved CSPM: A First-Order Accurate Estimate for Second Derivatives
4.1. Improving the Normalisation Matrix
In Section 3.1 we showed that the CSPM approximation for second derivatives is zeroth-order consistent. The crucial step in improving the accuracy is keeping the first-order terms in (16) separate from the second and higher-order terms. This gives:
so that, instead of (19), we find:
Accounting for the extra term in (28) leads to the adjustment of (20) to:
with and as defined in (21) and (22), respectively. Note that the nor- malisation matrix is the normalisation matrix used in CSPM. We propose to use the following approximation for second derivatives instead of (23):
with the normalisation matrix:
which we designate as improved CSPM (ICSPM). This method is slightly more expensive than CSPM, but it requires the same effort regarding matrix inversions. Yet, it follows directly from the definition of and (29) that this approximation is first-order accurate:
In Section 5 we explore the behaviour of ICSPM and compare it with the standard SPH, CSPM and MSPH approximations described in Sections 3.1, 3.2 and 2.3.
4.2. Accuracy versus Computational Cost
The methods described in Sections 3.1, 3.2 and 4 are more accurate, but also computationally more expensive than (11). Therefore, if the available computation time is limited, (11) is the approximation to go for. Also, if the domain has no boundaries, particles are distributed regularly, or if accuracy is simply not of great importance, (11) is a perfectly good estimate for the Laplacian. However, if accuracy is important and a somewhat higher computational cost is not an issue, CSPM, MSPH and ICSPM are attractive alternatives. Moreover, these methods give approximations for all second derivative-terms instead of just the Laplacian.
The MSPH method by Zhang and Batra   requires the inversion of 6 ´ 6- matrices (in two dimensions) and 10 ´ 10-matrices (in three dimensions). This makes the method computationally more expensive than CSPM, which splits the problem into smaller sets of equations, so that it only requires the inversion of 2 ´ 2 and 3 ´ 3-matrices (in two dimensions) and 3 ´ 3 and 6 ´ 6-matrices (in three dimensions). MSPH, however, is more accurate, as it is first-order accurate (26), whereas CSPM is only zeroth-order accurate (24).
With the improved normalisation step in Section 4, we have found a method that is computationally similar to CSPM―it only additionally requires the explicit computation of the sums in (31)―but which has the same order of accuracy as MSPH. We will verify this with numerical examples in Section 5.
In principle it is possible to increase the accuracy to any order. This might not be obvious with CSPM, but it is quite straightforward with MSPH. However, every increment in order introduces a number of extra equations to solve (simultaneously), so that it is not worth the effort for any order higher than one. Hence, if accurate approximations for second derivatives are required, we recommend ICSPM (30). This method is more accurate than CSPM, while the extra computational effort is negligible. MSPH has a similar accuracy, but is computationally more expensive.
5. Numerical Experiments
This section describes three numerical case studies to compare the various methods that were discussed previously. The first one yields the computation of a one-dimensional second derivative. In the second and third case study we consider a time-dependent problem on a two-dimensional domain. In all three experiments we use the Wendland kernel, which, in one-dimension, is given by:
where and is the radius of the support domain. As suggested by Dehnen and Aly  , we choose the smoothing length equal to twice the standard deviation, so that . In two dimensions, the Wendland kernel has a different form:
with . In our test cases we use these kernels with the ratio between particle distance and .
5.1. Second Derivative on a One-Dimensional Domain
We start with a one-dimensional domain , on which we compute the second derivative of the given function . This test case was also performed by Zhang and Batra  . The analytical solution is available and is equal to . Therefore, we can calculate the exact error, for which we use the infinity norm:
where is the total number of particles. We use , with the (uniform) inter-particle distance, equal to that of Zhang and Batra  . The two left panels of Figure 2 show the results with particles. The particles are distributed both uniformly (top) and pseudo-randomly (bottom), where in the latter case the boundary particles are kept on the boundary. As we can see in the figure, the particle distribution has no great effect on the results, except for the standard SPH discretisation.
(a) (b)(c) (d)
Figure 2. Approximating a second derivative with uniformly and pseudo-randomly distributed particles. A comparison between the analytical solution and four estimates. (a) Estimates (uniform); (b) Errors (uniform); (c) Estimates (random); (d) Errors (random).
We also calculate the errors according to (35) for several values of . These are shown in the two right panels of Figure 2. They clearly show the - convergence of SPH, the -convergence of CSPM and the -conver- gence of MSPH and ICSPM. Note that the error plots for MSPH and ICSPM coincide.
5.2. Heat Conduction on a Two-Dimensional Domain
Next, we compute the time-dependent temperature distribution on a two-di- mensional domain , a test case by Cleary and Monaghan  . The equation governing the conduction process is:
with parameter choices and . Initially, we have the following temperature distribution:
which is shown in Figure 3. We assume four isothermal walls with . The
Figure 3. The initial temperature distribution.
analytical solution of this problem is known and reads:
where . Since our main interest is the comparison between the various spatial discretisations, it is sufficient to use the explicit Euler time stepping scheme. The time step used depends on the spatial discretisation distance and is chosen as , where , with the reference volume, equal to the volume of an interior particle.
In this case study we use two expressions for the error. The first one is similar to (35), but looks at the temperature itself rather than the Laplacian:
The second one considers the relative errors at the positions of the particles instead of the absolute errors:
We choose the smoothing length to be , equal to that of Cleary and Monaghan  . We use pseudo-random particle distributions, so as to compare the methods with distributions that resemble realistic particle configurations. Starting with a uniform distribution, every particle is randomly shifted in both the horizontal and vertical direction, either to the left or to the right and up or down, with a maximum displacement of in either direction. Boundary particles are only shifted in one direction, such that they stay on the boundaries, while corner particles are entirely excluded from the shifting process.
The results in the left panel of Figure 4 show the temperature distributions at for the various methods, which we generated with a 17 ´ 17 particle distribution. Clearly, the standard SPH scheme performs worse. Yet, approximations with similar behaviour are often used in the literature. The CSPM, MSPH and ICSPM solutions are very close and in this figure cannot by distinguished by the naked eye. The left panel in Figure 5 shows, however, that there is a difference in behaviour between CSPM on the one hand, and MSPH and ICSPM on the other hand. Although initially the three methods have comparable errors, for
Figure 4. The temperature profile at for the particles for which and the particle distribution showing the volumes of particles, which form a Voronoi tessellation. (a) Temperature profile; (b) Particle distribution.
Figure 5. The absolute errors (39) and relative errors (40) for the heat conduction problem at for several values of . (a) Absolute errors; (b) Relative errors.
larger the CSPM method is less accurate than both MSPH and ICSPM.
The results in the right panel in Figure 5, for which we used the error in (40), show even more clearly that there is a significant difference in behaviour between CSPM and ICSPM (and MSPH). Clearly, the zeroth-order consistency of CSPM results in a zeroth-order convergence for the temperature as well, in contrast to ICSPM and MSPH, which are second-order convergent. MSPH is slightly more stable and accurate than ICSPM.
5.3. Heat Conduction in Slabs with Discontinuous Parameters
In the final test case we study the effect of discontinuous parameters. Cleary and Monaghan  describe the heat conduction in a slab of unit width which is periodic in the -direction. The slab consists of two different materials, each with their own set of parameters, touching each other at the interface at . Initially, the left half of the slab has zero temperature, , while the right half has temperature . Since is no longer constant on the whole domain, in the approximations it is replaced with:
When the temperature variation of the outermost points is small, the analytical solution can be approximated by   :
We start with a discontinuity only in the conductivity. For the left half it is set to , while for the right it is . The densities are and for both halves of the slab. We uniformly distribute 40 particles in the - direction and 20 in the -direction and we apply an explicit Euler time stepping scheme with time step size . Figure 6 shows the temperature distribution along the -direction for the ICSPM derived in this work. The left panel shows the temperature at , while at the right . As time progresses, the errors become smaller. This is confirmed by the errors in Table 1, which shows the maximum relative errors (40) for the Brookshaw method used
Figure 6. Heat distribution along the -direction. The conductivity of the left half of the slab is , while for the right it is . (a) Temperature at ; (b) Temperature at .
Table 1. Maximum relative errors for the particles for which , for various number of particles, when and . (a) Errors at ; (b) Errors at .
by Cleary and Monaghan and the ICSPM. To avoid dividing by small numbers we only consider particles for which .
Next, we consider the case where both the conductivity and the specific heat are discontinuous. We choose , while . The temperature distribution at is shown in the left panel of Figure 7. Again, we find that the temperature distribution is captured very well. The maximum relative errors are shown on the left in Table 2. The Brookshaw solution decreases only marginally when going from to . For the case , the error in the left table is one order lower for the ICSPM for these particular parameter choices. This behaviour seems to be coincidental, however, since the Brookshaw approximation and the ICSPM usually perform similarly in the infinite slab test case.
Finally, we also take the density to be discontinuous at the interface. We choose and , while , and . The temperature distribution at is shown in the right panel in Figure 7. The relative errors are shown on the right in Table 2. Just like in the previous cases, the temperature distribution is captured very well and the errors become smaller for increasing numbers of particles. Since the particles are uniformly distributed and there are no boundaries causing problems, both methods behave similarly. The test cases in this section show, however, that ICSPM can also handle discontinuous parameters.
Figure 7. Heat distribution along the -direction (a) at when and (b) at when and . (a) Temperature; (b) Temperature.
Table 2. Maximum relative errors for particles (a) at when and (b) at when and for various number of particles. (a) Errors; (b) Errors.
We studied approximations of second-derivative terms (or the Laplacian) in Smoothed Particle Hydrodynamics (SPH). Second derivatives appear in viscosity and conduction terms, and in Poisson equations. Traditionally, second derivatives have been approximated with the second derivatives of the kernel function, but these showed to be sensitive to particle disorder. Therefore various other methods and improvements have been suggested in the literature.
We proposed an improvement to the Corrective Smoothed Particle Method (CSPM) that increased the consistency of the estimates from zeroth-order to first-order. This method, called improved CSPM (ICSPM), was―together with several other methods―applied to one and two-dimensional test cases. The one- dimensional test case clearly showed the high accuracy of ICSPM compared to CSPM. In the first two-dimensional test case we solved the heat equation and compared errors of the temperature itself, rather than the Laplacian. With the absolute error, the difference between CSPM and ICSPM was already clearly visible, but when we considered the relative error the difference between the two methods was even more obvious. In the third test case we computed the temperature distribution for an infinite slab consisting of two different materials. The results showed that our improved approximation is perfectly capable of handling discontinuous parameters.
ICSPM is computationally only slightly more expensive than CSPM. Yet, its accuracy is similar to that of more expensive methods, such as the Finite Particle Method (FMP) and Modified Smoothed Particle Hydrodynamics (MSPH), although MSPH proved to be a bit more stable in specific cases.
The first author kindly acknowledges support from the Nederlandse Organisatie voor Wetenschappelijk Onderzoek (NWO) VICI Grant 639.033.008. We also like to thank the reviewers for their helpful comments.
 Gingold, R.A. and Monaghan, J.J. (1977) Smoothed Particle Hydrodynamics: Theory and Application to Nonspherical Stars. Monthly Notices of the Royal Astronomical Society, 181, 375-389.
 Colagrossi, A. and Landrini, M. (2003) Numerical Simulation of Interfacial Flows by Smoothed Particle Hydrodynamics. Journal of Computational Physics, 191, 448-475.
 Liu, X., Xu, H., Shao, S. and Lin, P. (2013) An Improved Incompressible SPH Model for Simulation of Wavestructure Interaction. Computers and Fluids, 71, 113-123.
 Randles, P.W. and Libersky, L.D. (1996) Smoothed Particle Hydrodynamics: Some Recent Improvements and Applications. Computer Methods in Applied Mechanics and Engineering, 139, 375-408.
 Belytschko, T., Krongauz, Y., Dolbow, J. and Gerlach, C. (1998) On the Completeness of Meshfree Particle Methods. International Journal for Numerical Methods in Engineering, 43, 785-819.
 Flebbe, O., Munzel, S., Herold, H., Riffert, H. and Ruder, H. (1994) Smoothed Particle Hydrodynamics: Physical Viscosity and the Simulation of Accretion Disks. The Astrophysical Journal, 431, 754-760. https://doi.org/10.1086/174526
 Fatehi, R. and Manzari, M.T. (2011) Error Estimation in Smoothed Particle Hydrodynamics and a New Scheme for Second Derivatives. Computers and Mathematics with Applications, 61, 482-498. https://doi.org/10.1016/j.camwa.2010.11.028
 Schwaiger, H.F. (2008) An Implicit Corrected SPH Formulation for the Thermal diffusion with Linear Free Surface Boundary Conditions. International Journal for Numerical Methods in Engineering, 75, 647-671. https://doi.org/10.1002/nme.2266
 Liu, W.K., Chen, Y., Jun, S., Chen, J.S., Belytschko, T., Pan, C., Uras, R.A. and Chang, C.T. (1996) Overview and Applications of the Reproducing Kernel Particle Methods. Archives of Computational Methods in Engineering, 3, 3-80.
 Liu, M.B., Liu, G.R. and Lam, K.Y. (2003) Constructing Smoothing Functions in Smoothed Particle Hydrodynamics with Applications. Journal of Computational and Applied Mathematics, 155, 263-284.
 Bonet, J. and Kulasegaram, S. (2002) A Simplified Approach to Enhance the Performance of Smooth Particle Hydrodynamics Methods. Applied Mathematics and Computation, 126, 133-155.
 Chen, J.K. and Beraun, J.E. (2000) A Generalized Smoothed Particle Hydrodynamics Method for Nonlinear Dynamic Problems. Computer Methods in Applied Mechanics and Engineering, 190, 225-239.
 Zhang, G.M. and Batra, R.C. (2004) Modified Smoothed Particle Hydrodynamics Method and Its Application to Transient Problems. Computational Mechanics, 34, 137-146. https://doi.org/10.1007/s00466-004-0561-5
 Korzilius, S.P., Schilders, W.H.A. and Anthonissen, M.J.H. (2013) An Improved Corrective Smoothed Particle Method Approximation for Second-Order Derivatives. Proceedings of the 8th International SPHERIC Workshop, Trondheim, 4-6 June 2013, 38-43.
 Korzilius, S.P., Schilders, W.H.A. and Anthonissen, M.J.H. (2015) An Improved CSPM Approximation for Multidimensional Second-Order Derivatives. Proceedings of the 10th International SPHERIC Workshop, Parma, 16-18 June 2015, 39-44.
 Liu, M.B. and Liu, G.R. (2010) Smoothed Particle Hydrodynamics (SPH): An Overview and Recent Developments. Archives of Computational Methods in Engineering, 17, 25-76.
 Chen, J.K., Beraun, J.E. and Carney, T.C. (1999) A Corrective Smoothed Particle Method for Boundary Value Problems in Heat Conduction. International Journal for Numerical Methods in Engineering, 46, 231-252.
 Chen, J.K., Beraun, J.E. and Jih, C.J. (1999) An Improvement for Tensile Instability in Smoothed Particle Hydrodynamics. Computational Mechanics, 23, 279-287.
 Chen, J.K., Beraun, J.E. and Jih, C.J. (1999) Completeness of Corrective Smoothed Particle Method for Linear Elastodynamics. Computational Mechanics, 24, 273-285.
 Chen, J.K., Beraun, J.E. and Jih, C.J. (2001) A Corrective Smoothed Particle Method for Transient Elastoplastic Dynamics. Computational Mechanics, 27, 177-187.
 Zhang, G.M. and Batra, R.C. (2007) Wave Propagation in Functionally Graded Materials by Modified Smoothed Particle Hydrodynamics (MSPH) Method. Journal of Computational Physics, 222, 374-390.
 Dehnen, W. and Aly, H. (2012) Improving Convergence in Smoothed Particle Hydrodynamics Simulations without Pairing Instability. Monthly Notices of the Royal Astronomical Society, 425, 1068-1082.